limdyn.F90 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. MODULE limdyn
  2. !!======================================================================
  3. !! *** MODULE limdyn ***
  4. !! Sea-Ice dynamics :
  5. !!======================================================================
  6. !! history : 1.0 ! 2002-08 (C. Ethe, G. Madec) original VP code
  7. !! 3.0 ! 2007-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle) LIM3: EVP-Cgrid
  8. !! 3.5 ! 2011-02 (G. Madec) dynamical allocation
  9. !!----------------------------------------------------------------------
  10. #if defined key_lim3
  11. !!----------------------------------------------------------------------
  12. !! 'key_lim3' : LIM3 sea-ice model
  13. !!----------------------------------------------------------------------
  14. !! lim_dyn : computes ice velocities
  15. !! lim_dyn_init : initialization and namelist read
  16. !!----------------------------------------------------------------------
  17. USE phycst ! physical constants
  18. USE dom_oce ! ocean space and time domain
  19. USE sbc_oce ! Surface boundary condition: ocean fields
  20. USE sbc_ice ! Surface boundary condition: ice fields
  21. USE ice ! LIM-3 variables
  22. USE dom_ice ! LIM-3 domain
  23. USE limrhg ! LIM-3 rheology
  24. USE lbclnk ! lateral boundary conditions - MPP exchanges
  25. USE lib_mpp ! MPP library
  26. USE wrk_nemo ! work arrays
  27. USE in_out_manager ! I/O manager
  28. USE prtctl ! Print control
  29. USE lib_fortran ! glob_sum
  30. USE timing ! Timing
  31. USE limcons ! conservation tests
  32. USE limvar
  33. IMPLICIT NONE
  34. PRIVATE
  35. PUBLIC lim_dyn ! routine called by ice_step
  36. !! * Substitutions
  37. # include "vectopt_loop_substitute.h90"
  38. !!----------------------------------------------------------------------
  39. !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
  40. !! $Id: limdyn.F90 7597 2017-01-20 18:40:06Z clem $
  41. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  42. !!----------------------------------------------------------------------
  43. CONTAINS
  44. SUBROUTINE lim_dyn( kt )
  45. !!-------------------------------------------------------------------
  46. !! *** ROUTINE lim_dyn ***
  47. !!
  48. !! ** Purpose : compute ice velocity and ocean-ice stress
  49. !!
  50. !! ** Method :
  51. !!
  52. !! ** Action : - Initialisation
  53. !! - Call of the dynamic routine for each hemisphere
  54. !! - computation of the stress at the ocean surface
  55. !! - treatment of the case if no ice dynamic
  56. !!------------------------------------------------------------------------------------
  57. INTEGER, INTENT(in) :: kt ! number of iteration
  58. !!
  59. INTEGER :: ji, jj, jl, ja ! dummy loop indices
  60. INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology
  61. REAL(wp) :: zcoef ! local scalar
  62. REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice
  63. REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask
  64. REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity
  65. !
  66. REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
  67. !!---------------------------------------------------------------------
  68. IF( nn_timing == 1 ) CALL timing_start('limdyn')
  69. CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
  70. CALL wrk_alloc( jpj, zswitch, zmsk )
  71. CALL lim_var_agg(1) ! aggregate ice categories
  72. IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only)
  73. IF( ln_limdyn ) THEN
  74. !
  75. ! conservation test
  76. IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
  77. u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
  78. v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
  79. ! Rheology (ice dynamics)
  80. ! ========
  81. ! Define the j-limits where ice rheology is computed
  82. ! ---------------------------------------------------
  83. IF( lk_mpp .OR. lk_mpp_rep ) THEN ! mpp: compute over the whole domain
  84. i_j1 = 1
  85. i_jpj = jpj
  86. IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
  87. CALL lim_rhg( i_j1, i_jpj )
  88. ELSE ! optimization of the computational area
  89. !
  90. DO jj = 1, jpj
  91. zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line
  92. zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line
  93. END DO
  94. IF( l_jeq ) THEN ! local domain include both hemisphere
  95. ! ! Rheology is computed in each hemisphere
  96. ! ! only over the ice cover latitude strip
  97. ! Northern hemisphere
  98. i_j1 = njeq
  99. i_jpj = jpj
  100. DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
  101. i_j1 = i_j1 + 1
  102. END DO
  103. i_j1 = MAX( 1, i_j1-2 )
  104. IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : NH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
  105. CALL lim_rhg( i_j1, i_jpj )
  106. !
  107. ! Southern hemisphere
  108. i_j1 = 1
  109. i_jpj = njeq
  110. DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
  111. i_jpj = i_jpj - 1
  112. END DO
  113. i_jpj = MIN( jpj, i_jpj+1 )
  114. IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : SH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
  115. !
  116. CALL lim_rhg( i_j1, i_jpj )
  117. !
  118. ELSE ! local domain extends over one hemisphere only
  119. ! ! Rheology is computed only over the ice cover
  120. ! ! latitude strip
  121. i_j1 = 1
  122. DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
  123. i_j1 = i_j1 + 1
  124. END DO
  125. i_j1 = MAX( 1, i_j1-2 )
  126. i_jpj = jpj
  127. DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
  128. i_jpj = i_jpj - 1
  129. END DO
  130. i_jpj = MIN( jpj, i_jpj+1)
  131. !
  132. IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : one hemisphere: i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
  133. !
  134. CALL lim_rhg( i_j1, i_jpj )
  135. !
  136. ENDIF
  137. !
  138. ENDIF
  139. ! computation of friction velocity
  140. ! --------------------------------
  141. ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
  142. zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
  143. zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
  144. ! frictional velocity at T-point
  145. zcoef = 0.5_wp * rn_cio
  146. DO jj = 2, jpjm1
  147. DO ji = fs_2, fs_jpim1 ! vector opt.
  148. ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
  149. & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1)
  150. END DO
  151. END DO
  152. !
  153. ! conservation test
  154. IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
  155. !
  156. ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean
  157. !
  158. zcoef = SQRT( 0.5_wp ) * r1_rau0
  159. DO jj = 2, jpjm1
  160. DO ji = fs_2, fs_jpim1 ! vector opt.
  161. ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
  162. & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1)
  163. END DO
  164. END DO
  165. !
  166. ENDIF
  167. CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point
  168. IF(ln_ctl) THEN ! Control print
  169. CALL prt_ctl_info(' ')
  170. CALL prt_ctl_info(' - Cell values : ')
  171. CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
  172. CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :')
  173. CALL prt_ctl(tab2d_1=divu_i , clinfo1=' lim_dyn : divu_i :')
  174. CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :')
  175. CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :')
  176. CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_dyn : cell area :')
  177. CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :')
  178. CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :')
  179. CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_dyn : vt_s :')
  180. CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn : stress1_i :')
  181. CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn : stress2_i :')
  182. CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn : stress12_i:')
  183. DO jl = 1, jpl
  184. CALL prt_ctl_info(' ')
  185. CALL prt_ctl_info(' - Category : ', ivar1=jl)
  186. CALL prt_ctl_info(' ~~~~~~~~~~')
  187. CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_dyn : a_i : ')
  188. CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_dyn : ht_i : ')
  189. CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_dyn : ht_s : ')
  190. CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_dyn : v_i : ')
  191. CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_dyn : v_s : ')
  192. CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_dyn : e_s : ')
  193. CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_dyn : t_su : ')
  194. CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_dyn : t_snow : ')
  195. CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_dyn : sm_i : ')
  196. CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_dyn : smv_i : ')
  197. DO ja = 1, nlay_i
  198. CALL prt_ctl_info(' ')
  199. CALL prt_ctl_info(' - Layer : ', ivar1=ja)
  200. CALL prt_ctl_info(' ~~~~~~~')
  201. CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn : t_i : ')
  202. CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn : e_i : ')
  203. END DO
  204. END DO
  205. ENDIF
  206. !
  207. CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
  208. CALL wrk_dealloc( jpj, zswitch, zmsk )
  209. !
  210. IF( nn_timing == 1 ) CALL timing_stop('limdyn')
  211. END SUBROUTINE lim_dyn
  212. SUBROUTINE lim_dyn_init
  213. !!-------------------------------------------------------------------
  214. !! *** ROUTINE lim_dyn_init ***
  215. !!
  216. !! ** Purpose : Physical constants and parameters linked to the ice
  217. !! dynamics
  218. !!
  219. !! ** Method : Read the namicedyn namelist and check the ice-dynamic
  220. !! parameter values called at the first timestep (nit000)
  221. !!
  222. !! ** input : Namelist namicedyn
  223. !!-------------------------------------------------------------------
  224. INTEGER :: ios ! Local integer output status for namelist read
  225. NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, &
  226. & nn_nevp, rn_relast
  227. !!-------------------------------------------------------------------
  228. REWIND( numnam_ice_ref ) ! Namelist namicedyn in reference namelist : Ice dynamics
  229. READ ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
  230. 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
  231. REWIND( numnam_ice_cfg ) ! Namelist namicedyn in configuration namelist : Ice dynamics
  232. READ ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
  233. 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
  234. IF(lwm) WRITE ( numoni, namicedyn )
  235. IF(lwp) THEN ! control print
  236. WRITE(numout,*)
  237. WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
  238. WRITE(numout,*) '~~~~~~~~~~~~'
  239. WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr
  240. WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf
  241. WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg
  242. WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio
  243. WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar
  244. WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg
  245. WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl
  246. WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc
  247. WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp
  248. WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast
  249. ENDIF
  250. !
  251. usecc2 = 1._wp / ( rn_ecc * rn_ecc )
  252. rhoco = rau0 * rn_cio
  253. !
  254. END SUBROUTINE lim_dyn_init
  255. #else
  256. !!----------------------------------------------------------------------
  257. !! Default option Empty module NO LIM sea-ice model
  258. !!----------------------------------------------------------------------
  259. CONTAINS
  260. SUBROUTINE lim_dyn ! Empty routine
  261. END SUBROUTINE lim_dyn
  262. #endif
  263. !!======================================================================
  264. END MODULE limdyn