limhdf.F90 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317
  1. MODULE limhdf
  2. !!======================================================================
  3. !! *** MODULE limhdf ***
  4. !! LIM ice model : horizontal diffusion of sea-ice quantities
  5. !!======================================================================
  6. !! History : LIM ! 2000-01 (LIM) Original code
  7. !! - ! 2001-05 (G. Madec, R. Hordoir) opa norm
  8. !! 1.0 ! 2002-08 (C. Ethe) F90, free form
  9. !! 3.0 ! 2015-08 (O. Tintó and M. Castrillo) added lim_hdf (multiple)
  10. !!----------------------------------------------------------------------
  11. #if defined key_lim3
  12. !!----------------------------------------------------------------------
  13. !! 'key_lim3' LIM3 sea-ice model
  14. !!----------------------------------------------------------------------
  15. !! lim_hdf : diffusion trend on sea-ice variable
  16. !! lim_hdf_init : initialisation of diffusion trend on sea-ice variable
  17. !!----------------------------------------------------------------------
  18. USE dom_oce ! ocean domain
  19. USE ice ! LIM-3: ice variables
  20. USE lbclnk ! lateral boundary condition - MPP exchanges
  21. USE lib_mpp ! MPP library
  22. USE wrk_nemo ! work arrays
  23. USE prtctl ! Print control
  24. USE in_out_manager ! I/O manager
  25. USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
  26. IMPLICIT NONE
  27. PRIVATE
  28. PUBLIC lim_hdf ! called by lim_trp
  29. PUBLIC lim_hdf_init ! called by sbc_lim_init
  30. LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call)
  31. INTEGER :: nn_convfrq !: convergence check frequency of the Crant-Nicholson scheme
  32. REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: efact ! metric coefficient
  33. !! * Substitution
  34. # include "vectopt_loop_substitute.h90"
  35. !!----------------------------------------------------------------------
  36. !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
  37. !! $Id: limhdf.F90 4990 2014-12-15 16:42:49Z timgraham $
  38. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  39. !!----------------------------------------------------------------------
  40. CONTAINS
  41. SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i )
  42. !!-------------------------------------------------------------------
  43. !! *** ROUTINE lim_hdf ***
  44. !!
  45. !! ** purpose : Compute and add the diffusive trend on sea-ice variables
  46. !!
  47. !! ** method : Second order diffusive operator evaluated using a
  48. !! Cranck-Nicholson time Scheme.
  49. !!
  50. !! ** Action : update ptab with the diffusive contribution
  51. !!-------------------------------------------------------------------
  52. INTEGER :: jpl, nlay_i, isize, ihdf_vars
  53. REAL(wp), DIMENSION(:,:,:), INTENT( inout ),TARGET :: ptab ! Field on which the diffusion is applied
  54. !
  55. INTEGER :: ji, jj, jk, jl , jm ! dummy loop indices
  56. INTEGER :: iter, ierr ! local integers
  57. REAL(wp) :: zrlxint ! local scalars
  58. REAL(wp), POINTER , DIMENSION ( : ) :: zconv ! local scalars
  59. REAL(wp), POINTER , DIMENSION(:,:,:) :: zrlx,zdiv0, ztab0
  60. REAL(wp), POINTER , DIMENSION(:,:) :: zflu, zflv, zdiv
  61. CHARACTER(lc) :: charout ! local character
  62. REAL(wp), PARAMETER :: zrelax = 0.5_wp ! relaxation constant for iterative procedure
  63. REAL(wp), PARAMETER :: zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
  64. INTEGER , PARAMETER :: its = 100 ! Maximum number of iteration
  65. !!-------------------------------------------------------------------
  66. TYPE(arrayptr) , ALLOCATABLE, DIMENSION(:) :: pt2d_array, zrlx_array
  67. CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) :: type_array ! define the nature of ptab array grid-points
  68. ! ! = T , U , V , F , W and I points
  69. REAL(wp) , ALLOCATABLE, DIMENSION(:) :: psgn_array ! =-1 the sign change across the north fold boundary
  70. !!---------------------------------------------------------------------
  71. ! !== Initialisation ==!
  72. ! +1 open water diffusion
  73. isize = jpl*(ihdf_vars+nlay_i)+1
  74. ALLOCATE( zconv (isize) )
  75. ALLOCATE( pt2d_array(isize) , zrlx_array(isize) )
  76. ALLOCATE( type_array(isize) )
  77. ALLOCATE( psgn_array(isize) )
  78. CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
  79. CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv )
  80. DO jk= 1 , isize
  81. pt2d_array(jk)%pt2d=>ptab(:,:,jk)
  82. zrlx_array(jk)%pt2d=>zrlx(:,:,jk)
  83. type_array(jk)='T'
  84. psgn_array(jk)=1.
  85. END DO
  86. !
  87. IF( linit ) THEN ! Metric coefficient (compute at the first call and saved in efact)
  88. ALLOCATE( efact(jpi,jpj) , STAT=ierr )
  89. IF( lk_mpp ) CALL mpp_sum( ierr )
  90. IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
  91. DO jj = 2, jpjm1
  92. DO ji = fs_2 , fs_jpim1 ! vector opt.
  93. efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
  94. END DO
  95. END DO
  96. linit = .FALSE.
  97. ENDIF
  98. ! ! Time integration parameters
  99. !
  100. zflu (jpi,: ) = 0._wp
  101. zflv (jpi,: ) = 0._wp
  102. DO jk=1 , isize
  103. ztab0(:, : , jk ) = ptab(:,:,jk) ! Arrays initialization
  104. zdiv0(:, 1 , jk ) = 0._wp
  105. zdiv0(:,jpj, jk ) = 0._wp
  106. zdiv0(1, :, jk ) = 0._wp
  107. zdiv0(jpi,:, jk ) = 0._wp
  108. END DO
  109. zconv = 1._wp !== horizontal diffusion using a Crant-Nicholson scheme ==!
  110. iter = 0
  111. !
  112. DO WHILE( MAXVAL(zconv(:)) > ( 2._wp * 1.e-04 ) .AND. iter <= its ) ! Sub-time step loop
  113. !
  114. iter = iter + 1 ! incrementation of the sub-time step number
  115. !
  116. DO jk = 1 , isize
  117. jl = (jk-1) /( ihdf_vars+nlay_i)+1
  118. IF (zconv(jk) > ( 2._wp * 1.e-04 )) THEN
  119. DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction
  120. DO ji = 1 , fs_jpim1 ! vector opt.
  121. zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
  122. zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
  123. END DO
  124. END DO
  125. !
  126. DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes
  127. DO ji = fs_2 , fs_jpim1 ! vector opt.
  128. zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
  129. END DO
  130. END DO
  131. !
  132. IF( iter == 1 ) zdiv0(:,:,jk) = zdiv(:,:) ! save the 1st evaluation of the diffusive trend in zdiv0
  133. !
  134. DO jj = 2, jpjm1 ! iterative evaluation
  135. DO ji = fs_2 , fs_jpim1 ! vector opt.
  136. zrlxint = ( ztab0(ji,jj,jk) &
  137. & + rdt_ice * ( zalfa * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj,jk) ) &
  138. & + ( 1.0 - zalfa ) * zdiv0(ji,jj,jk) ) &
  139. & ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
  140. zrlx(ji,jj,jk) = ptab(ji,jj,jk) + zrelax * ( zrlxint - ptab(ji,jj,jk) )
  141. END DO
  142. END DO
  143. END IF
  144. END DO
  145. CALL lbc_lnk_multi( zrlx_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
  146. !
  147. IF ( MOD( iter-1 , nn_convfrq ) == 0 ) THEN !Convergence test every nn_convfrq iterations (perf. optimization )
  148. DO jk=1,isize
  149. zconv(jk) = 0._wp ! convergence test
  150. DO jj = 2, jpjm1
  151. DO ji = fs_2, fs_jpim1
  152. zconv(jk) = MAX( zconv(jk), ABS( zrlx(ji,jj,jk) - ptab(ji,jj,jk) ) )
  153. END DO
  154. END DO
  155. END DO
  156. IF( lk_mpp ) CALL mpp_max_multiple( zconv , isize ) ! max over the global domain for all the variables
  157. ENDIF
  158. !
  159. DO jk=1,isize
  160. ptab(:,:,jk) = zrlx(:,:,jk)
  161. END DO
  162. !
  163. END DO ! end of sub-time step loop
  164. ! -----------------------
  165. !!! final step (clem) !!!
  166. DO jk = 1, isize
  167. jl = (jk-1) /( ihdf_vars+nlay_i)+1
  168. DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction
  169. DO ji = 1 , fs_jpim1 ! vector opt.
  170. zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
  171. zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
  172. END DO
  173. END DO
  174. !
  175. DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes
  176. DO ji = fs_2 , fs_jpim1 ! vector opt.
  177. zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
  178. ptab(ji,jj,jk) = ztab0(ji,jj,jk) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj,jk) )
  179. END DO
  180. END DO
  181. END DO
  182. CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
  183. !!! final step (clem) !!!
  184. ! -----------------------
  185. ! IF(ln_ctl) THEN
  186. ! DO jk = 1 , isize
  187. ! zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk)
  188. ! WRITE(charout,FMT="('lim_hdf : zconv =',D23.16, ' iter =',I4)") zconv, iter
  189. ! CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout )
  190. ! END DO
  191. ! ENDIF
  192. !
  193. CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
  194. CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv )
  195. DEALLOCATE( zconv )
  196. DEALLOCATE( pt2d_array , zrlx_array )
  197. DEALLOCATE( type_array )
  198. DEALLOCATE( psgn_array )
  199. !
  200. END SUBROUTINE lim_hdf
  201. SUBROUTINE lim_hdf_init
  202. !!-------------------------------------------------------------------
  203. !! *** ROUTINE lim_hdf_init ***
  204. !!
  205. !! ** Purpose : Initialisation of horizontal diffusion of sea-ice
  206. !!
  207. !! ** Method : Read the namicehdf namelist
  208. !!
  209. !! ** input : Namelist namicehdf
  210. !!-------------------------------------------------------------------
  211. INTEGER :: ios ! Local integer output status for namelist read
  212. NAMELIST/namicehdf/ nn_ahi0, rn_ahi0_ref, nn_convfrq
  213. INTEGER :: ji, jj
  214. REAL(wp) :: za00, zd_max
  215. !!-------------------------------------------------------------------
  216. !
  217. REWIND( numnam_ice_ref ) ! Namelist namicehdf in reference namelist : Ice horizontal diffusion
  218. READ ( numnam_ice_ref, namicehdf, IOSTAT = ios, ERR = 901)
  219. 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in reference namelist', lwp )
  220. REWIND( numnam_ice_cfg ) ! Namelist namicehdf in configuration namelist : Ice horizontal diffusion
  221. READ ( numnam_ice_cfg, namicehdf, IOSTAT = ios, ERR = 902 )
  222. 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in configuration namelist', lwp )
  223. IF(lwm) WRITE ( numoni, namicehdf )
  224. !
  225. IF(lwp) THEN ! control print
  226. WRITE(numout,*)
  227. WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion'
  228. WRITE(numout,*) '~~~~~~~~~~~'
  229. WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0
  230. WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref
  231. WRITE(numout,*) ' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq
  232. ENDIF
  233. !
  234. ! Diffusion coefficients
  235. SELECT CASE( nn_ahi0 )
  236. CASE( -1 )
  237. ahiu(:,:) = 0._wp
  238. ahiv(:,:) = 0._wp
  239. IF(lwp) WRITE(numout,*) ''
  240. IF(lwp) WRITE(numout,*) ' No sea-ice diffusion applied'
  241. CASE( 0 )
  242. ahiu(:,:) = rn_ahi0_ref
  243. ahiv(:,:) = rn_ahi0_ref
  244. IF(lwp) WRITE(numout,*) ''
  245. IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref'
  246. CASE( 1 )
  247. zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
  248. IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
  249. ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2
  250. ! (60deg = min latitude for ice cover)
  251. ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
  252. IF(lwp) WRITE(numout,*) ''
  253. IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
  254. IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp
  255. CASE( 2 )
  256. zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
  257. IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
  258. za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2
  259. ! (60deg = min latitude for ice cover)
  260. DO jj = 1, jpj
  261. DO ji = 1, jpi
  262. ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
  263. ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
  264. END DO
  265. END DO
  266. !
  267. IF(lwp) WRITE(numout,*) ''
  268. IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1'
  269. IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
  270. END SELECT
  271. !
  272. END SUBROUTINE lim_hdf_init
  273. #else
  274. !!----------------------------------------------------------------------
  275. !! Default option Dummy module NO LIM sea-ice model
  276. !!----------------------------------------------------------------------
  277. #endif
  278. !!======================================================================
  279. END MODULE limhdf