Changeset 14834 for NEMO/trunk/src/OCE/ZDF
- Timestamp:
- 2021-05-11T11:24:44+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/ZDF
- Files:
-
- 12 edited
-
zdfddm.F90 (modified) (5 diffs)
-
zdfdrg.F90 (modified) (3 diffs)
-
zdfevd.F90 (modified) (6 diffs)
-
zdfgls.F90 (modified) (26 diffs)
-
zdfiwm.F90 (modified) (16 diffs)
-
zdfmfc.F90 (modified) (7 diffs)
-
zdfmxl.F90 (modified) (6 diffs)
-
zdfphy.F90 (modified) (8 diffs)
-
zdfric.F90 (modified) (2 diffs)
-
zdfsh2.F90 (modified) (3 diffs)
-
zdfswm.F90 (modified) (1 diff)
-
zdftke.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdfddm.F90
r14053 r14834 83 83 REAL(dp) :: zavfs ! - - 84 84 REAL(wp) :: zavdt, zavds ! - - 85 REAL(wp), DIMENSION( jpi,jpj) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd385 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 86 86 !!---------------------------------------------------------------------- 87 87 ! … … 95 95 !!gm and many acces in memory 96 96 97 DO_2D( 1, 1, 1,1 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==!97 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==! 98 98 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 99 99 !!gm please, use e3w at Kmm below … … 111 111 END_2D 112 112 113 DO_2D( 1, 1, 1,1 ) !== indicators ==!113 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== indicators ==! 114 114 ! stability indicator: msks=1 if rn2>0; 0 elsewhere 115 115 IF( rn2(ji,jj,jk) + 1.e-12 <= 0. ) THEN ; zmsks(ji,jj) = 0._wp … … 135 135 END_2D 136 136 ! mask zmsk in order to have avt and avs masked 137 zmsks(:,:) = zmsks(:,:) * wmask( :,:,jk)137 zmsks(:,:) = zmsks(:,:) * wmask(A2D(nn_hls),jk) 138 138 139 139 … … 141 141 ! ------------------ 142 142 ! Constant eddy coefficient: reset to the background value 143 DO_2D ( 1, 1, 1,1 )143 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 144 144 zinr = 1._wp / zrau(ji,jj) 145 145 ! salt fingering -
NEMO/trunk/src/OCE/ZDF/zdfdrg.F90
r13558 r14834 117 117 ! 118 118 IF( l_log_not_linssh ) THEN !== "log layer" ==! compute Cd and -Cd*|U| 119 DO_2D ( 0, 0, 0, 0)119 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 120 120 imk = k_mk(ji,jj) ! ocean bottom level at t-points 121 121 zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm) ! 2 x velocity at t-point … … 129 129 END_2D 130 130 ELSE !== standard Cd ==! 131 DO_2D ( 0, 0, 0, 0)131 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 132 132 imk = k_mk(ji,jj) ! ocean bottom level at t-points 133 133 zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm) ! 2 x velocity at t-point … … 432 432 l_log_not_linssh = .FALSE. !- don't update Cd at each time step 433 433 ! 434 DO_2D( 1, 1, 1, 1) ! pCd0 = mask (and boosted) logarithmic drag coef.434 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! pCd0 = mask (and boosted) logarithmic drag coef. 435 435 zzz = 0.5_wp * e3t_0(ji,jj,k_mk(ji,jj)) 436 436 zcd = ( vkarmn / LOG( zzz / rn_z0 ) )**2 -
NEMO/trunk/src/OCE/ZDF/zdfevd.F90
r13295 r14834 62 62 ! 63 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zavt_evd, zavm_evd 64 ! NOTE: [tiling] use a SAVE array to store diagnostics, then send after all tiles are finished. This is necessary because p_avt/p_avm are modified on adjacent tiles when using nn_hls > 1. zavt_evd/zavm_evd are then zero on some points when subsequently calculated for these tiles. 65 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: zavt_evd, zavm_evd 65 66 !!---------------------------------------------------------------------- 66 67 ! 67 IF( kt == nit000 ) THEN 68 IF(lwp) WRITE(numout,*) 69 IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)' 70 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 71 IF(lwp) WRITE(numout,*) 68 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 69 IF( kt == nit000 ) THEN 70 IF(lwp) WRITE(numout,*) 71 IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)' 72 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 73 IF(lwp) WRITE(numout,*) 74 ENDIF 75 76 ALLOCATE( zavt_evd(jpi,jpj,jpk) ) 77 IF( nn_evdm == 1 ) ALLOCATE( zavm_evd(jpi,jpj,jpk) ) 72 78 ENDIF 73 79 ! 74 80 ! 75 zavt_evd(:,:,:) = p_avt(:,:,:) ! set avt prior to evd application 81 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 82 zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) ! set avt prior to evd application 83 END_3D 76 84 ! 77 85 SELECT CASE ( nn_evdm ) … … 79 87 CASE ( 1 ) !== enhance tracer & momentum Kz ==! (if rn2<-1.e-12) 80 88 ! 81 zavm_evd(:,:,:) = p_avm(:,:,:) ! set avm prior to evd application 89 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 90 zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) ! set avm prior to evd application 91 END_3D 82 92 ! 83 93 !! change last digits results … … 87 97 ! END WHERE 88 98 ! 89 DO_3D ( 0, 0, 0, 0, 1, jpkm1 )99 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 90 100 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 91 101 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) … … 94 104 END_3D 95 105 ! 96 zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:) ! change in avm due to evd 97 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 106 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 107 zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) - zavm_evd(ji,jj,jk) ! change in avm due to evd 108 END_3D 109 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 110 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 111 DEALLOCATE( zavm_evd ) 112 ENDIF 98 113 ! 99 114 CASE DEFAULT !== enhance tracer Kz ==! (if rn2<-1.e-12) … … 103 118 ! END WHERE 104 119 105 DO_3D ( 0, 0, 0, 0, 1, jpkm1 )120 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 106 121 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 107 122 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) … … 110 125 END SELECT 111 126 ! 112 zavt_evd(:,:,:) = p_avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd 113 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 127 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 128 zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) - zavt_evd(ji,jj,jk) ! change in avt due to evd 129 END_3D 130 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 131 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 132 DEALLOCATE( zavt_evd ) 133 ENDIF 114 134 IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_evd, zavt_evd ) 115 135 ! -
NEMO/trunk/src/OCE/ZDF/zdfgls.F90
r14156 r14834 137 137 USE zdf_oce , ONLY : en, avtb, avmb ! ocean vertical physics 138 138 !! 139 INTEGER , INTENT(in ) :: kt ! ocean time step140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices141 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)139 INTEGER , INTENT(in ) :: kt ! ocean time step 140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 141 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 143 143 ! 144 144 INTEGER :: ji, jj, jk ! dummy loop arguments … … 151 151 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 152 152 REAL(wp) :: zmsku, zmskv ! - - 153 REAL(wp), DIMENSION( jpi,jpj) :: zdep154 REAL(wp), DIMENSION( jpi,jpj) :: zkar155 REAL(wp), DIMENSION( jpi,jpj) :: zflxs! Turbulence fluxed induced by internal waves156 REAL(wp), DIMENSION( jpi,jpj) :: zhsro! Surface roughness (surface waves)157 REAL(wp), DIMENSION( jpi,jpj) :: zice_fra! Tapering of wave breaking under sea ice158 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eb! tke at time before159 REAL(wp), DIMENSION( jpi,jpj,jpk) :: hmxl_b! mixing length at time before160 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eps! dissipation rate161 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwall_psi! Wall function use in the wb case (ln_sigpsi)162 REAL(wp), DIMENSION( jpi,jpj,jpk) :: psi! psi at time now163 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix164 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zstt, zstm! stability function on tracer and momentum153 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdep 154 REAL(wp), DIMENSION(A2D(nn_hls)) :: zkar 155 REAL(wp), DIMENSION(A2D(nn_hls)) :: zflxs ! Turbulence fluxed induced by internal waves 156 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhsro ! Surface roughness (surface waves) 157 REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra ! Tapering of wave breaking under sea ice 158 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eb ! tke at time before 159 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: hmxl_b ! mixing length at time before 160 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eps ! dissipation rate 161 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 162 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: psi ! psi at time now 163 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix 164 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zstt, zstm ! stability function on tracer and momentum 165 165 !!-------------------------------------------------------------------- 166 166 ! 167 167 ! Preliminary computing 168 169 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 170 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 171 ustar2_bot (:,:) = 0._wp 168 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 169 ustar2_surf(ji,jj) = 0._wp ; ustar2_top(ji,jj) = 0._wp ; ustar2_bot(ji,jj) = 0._wp 170 END_2D 171 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 172 psi(ji,jj,jk) = 0._wp ; zwall_psi(ji,jj,jk) = 0._wp 173 END_3D 172 174 173 175 SELECT CASE ( nn_z0_ice ) 174 176 CASE( 0 ) ; zice_fra(:,:) = 0._wp 175 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i( :,:) * 10._wp )176 CASE( 2 ) ; zice_fra(:,:) = fr_i( :,:)177 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i( :,:) , 1._wp )177 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 178 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 179 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 178 180 END SELECT 179 181 180 182 ! Compute surface, top and bottom friction at T-points 181 DO_2D ( 0, 0, 0, 0) !== surface ocean friction183 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== surface ocean friction 182 184 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) ! surface friction 183 185 END_2D … … 186 188 ! 187 189 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 DO_2D ( 0, 0, 0, 0 )! bottom friction (explicit before friction)190 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction (explicit before friction) 189 191 zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 190 192 zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) … … 193 195 END_2D 194 196 IF( ln_isfcav ) THEN 195 DO_2D ( 0, 0, 0, 0) ! top friction197 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 196 198 zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 197 199 zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) … … 206 208 zhsro(:,:) = rn_hsro 207 209 CASE ( 1 ) ! Standard Charnock formula 208 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf( :,:) , rn_hsro )210 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro ) 209 211 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 210 212 !!gm faster coding : the 2 comment lines should be used 211 213 !!gm zcof = 2._wp * 0.6_wp / 28._wp 212 214 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 213 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 214 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 215 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 216 zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) ) ! Wave age (eq. 10) 217 zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 218 END_2D 215 219 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 216 zhsro(:,:) = MAX(rn_frac_hs * hsw( :,:), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )220 zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 217 221 END SELECT 218 222 ! 219 223 ! adapt roughness where there is sea ice 220 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 221 ! 222 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! 224 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 225 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1) + & 226 & (1._wp - tmask(ji,jj,1))*rn_hsro 227 END_2D 228 ! 229 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !== Compute dissipation rate ==! 223 230 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 224 231 END_3D 225 232 226 233 ! Save tke at before time step 227 eb (:,:,:) = en (:,:,:) 228 hmxl_b(:,:,:) = hmxl_n(:,:,:) 234 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 235 eb (ji,jj,jk) = en (ji,jj,jk) 236 hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk) 237 END_3D 229 238 230 239 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 231 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )240 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 232 241 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 233 242 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) … … 250 259 ! Warning : after this step, en : right hand side of the matrix 251 260 252 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )261 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 253 262 ! 254 263 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 303 312 ! 304 313 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 305 ! First level 306 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 ) 307 zd_lw(:,:,1) = en(:,:,1) 308 zd_up(:,:,1) = 0._wp 309 zdiag(:,:,1) = 1._wp 310 ! 311 ! One level below 312 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 313 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 314 zd_lw(:,:,2) = 0._wp 315 zd_up(:,:,2) = 0._wp 316 zdiag(:,:,2) = 1._wp 317 ! 318 ! 314 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 315 ! First level 316 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 317 zd_lw(ji,jj,1) = en(ji,jj,1) 318 zd_up(ji,jj,1) = 0._wp 319 zdiag(ji,jj,1) = 1._wp 320 ! 321 ! One level below 322 en (ji,jj,2) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1 & 323 & * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf) )**r2_3 ) 324 zd_lw(ji,jj,2) = 0._wp 325 zd_up(ji,jj,2) = 0._wp 326 zdiag(ji,jj,2) = 1._wp 327 END_2D 328 ! 329 ! 319 330 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) 320 ! 321 ! Dirichlet conditions at k=1 322 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin ) 323 zd_lw(:,:,1) = en(:,:,1) 324 zd_up(:,:,1) = 0._wp 325 zdiag(:,:,1) = 1._wp 326 ! 327 ! at k=2, set de/dz=Fw 328 !cbr 329 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 330 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 331 zd_lw(ji,jj,2) = 0._wp 332 END_2D 333 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 334 zflxs(:,:) = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 335 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 331 ! 332 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 333 ! Dirichlet conditions at k=1 334 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 335 zd_lw(ji,jj,1) = en(ji,jj,1) 336 zd_up(ji,jj,1) = 0._wp 337 zdiag(ji,jj,1) = 1._wp 338 ! 339 ! at k=2, set de/dz=Fw 340 !cbr 341 ! zdiag zd_lw not defined/used on the halo 342 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 343 zd_lw(ji,jj,2) = 0._wp 344 ! 345 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) )) 346 zflxs(ji,jj) = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) & 347 & * ( ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj) )**(1.5_wp*ra_sf) 336 348 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 337 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 338 ! 339 ! 349 en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 350 END_2D 351 ! 352 ! 340 353 END SELECT 341 354 … … 348 361 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 349 362 ! ! Balance between the production and the dissipation terms 350 DO_2D ( 0, 0, 0, 0)363 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 351 364 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 352 365 !! With thick deep ocean level thickness, this may be quite large, no ??? … … 365 378 END_2D 366 379 ! 380 ! NOTE: ctl_stop with ln_isfcav when using GLS 367 381 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 368 DO_2D ( 0, 0, 0, 0)382 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 369 383 itop = mikt(ji,jj) ! k top w-point 370 384 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 384 398 CASE ( 1 ) ! Neumman boundary condition 385 399 ! 386 DO_2D ( 0, 0, 0, 0)400 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 387 401 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 388 402 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 398 412 en (ji,jj,ibot) = z_en 399 413 END_2D 414 ! NOTE: ctl_stop with ln_isfcav when using GLS 400 415 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 401 DO_2D ( 0, 0, 0, 0)416 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 402 417 itop = mikt(ji,jj) ! k top w-point 403 418 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 420 435 ! ---------------------------------------------------------- 421 436 ! 422 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1437 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 423 438 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 424 439 END_3D 425 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1440 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 426 441 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 427 442 END_3D 428 DO_3DS ( 0, 0, 0, 0, jpkm1, 2, -1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk443 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 429 444 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 430 445 END_3D 431 446 ! ! set the minimum value of tke 432 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 447 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 448 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) 449 END_3D 433 450 434 451 !!----------------------------------------!! … … 441 458 ! 442 459 CASE( 0 ) ! k-kl (Mellor-Yamada) 443 DO_3D( 0, 0, 0, 0, 2, jpkm1 )460 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 444 461 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 445 462 END_3D 446 463 ! 447 464 CASE( 1 ) ! k-eps 448 DO_3D( 0, 0, 0, 0, 2, jpkm1 )465 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 449 466 psi(ji,jj,jk) = eps(ji,jj,jk) 450 467 END_3D 451 468 ! 452 469 CASE( 2 ) ! k-w 453 DO_3D( 0, 0, 0, 0, 2, jpkm1 )470 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 454 471 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 455 472 END_3D 456 473 ! 457 474 CASE( 3 ) ! generic 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 )475 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 459 476 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 460 477 END_3D … … 469 486 ! Warning : after this step, en : right hand side of the matrix 470 487 471 DO_3D( 0, 0, 0, 0, 2, jpkm1 )488 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 472 489 ! 473 490 ! psi / k … … 516 533 CASE ( 0 ) ! Dirichlet boundary conditions 517 534 ! 518 ! Surface value 519 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 520 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 521 zd_lw(:,:,1) = psi(:,:,1) 522 zd_up(:,:,1) = 0._wp 523 zdiag(:,:,1) = 1._wp 524 ! 525 ! One level below 526 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 527 zdep (:,:) = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 528 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 529 zd_lw(:,:,2) = 0._wp 530 zd_up(:,:,2) = 0._wp 531 zdiag(:,:,2) = 1._wp 535 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 536 ! Surface value 537 zdep (ji,jj) = zhsro(ji,jj) * rl_sf ! Cosmetic 538 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 539 zd_lw(ji,jj,1) = psi(ji,jj,1) 540 zd_up(ji,jj,1) = 0._wp 541 zdiag(ji,jj,1) = 1._wp 542 ! 543 ! One level below 544 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) ))) 545 zdep (ji,jj) = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj) 546 psi (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 547 zd_lw(ji,jj,2) = 0._wp 548 zd_up(ji,jj,2) = 0._wp 549 zdiag(ji,jj,2) = 1._wp 550 END_2D 532 551 ! 533 552 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 534 553 ! 535 ! Surface value: Dirichlet536 zdep (:,:) = zhsro(:,:) * rl_sf537 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)538 zd_lw(:,:,1) = psi(:,:,1)539 zd_up(:,:,1) = 0._wp540 zdiag(:,:,1) = 1._wp541 !542 ! Neumann condition at k=2543 DO_2D( 0, 0, 0, 0 ) !zdiag zd_lw not defined/used on the halo554 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 555 ! Surface value: Dirichlet 556 zdep (ji,jj) = zhsro(ji,jj) * rl_sf 557 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 558 zd_lw(ji,jj,1) = psi(ji,jj,1) 559 zd_up(ji,jj,1) = 0._wp 560 zdiag(ji,jj,1) = 1._wp 561 ! 562 ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo 544 563 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 545 564 zd_lw(ji,jj,2) = 0._wp 565 ! 566 ! Set psi vertical flux at the surface: 567 zkar (ji,jj) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope 568 zdep (ji,jj) = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf) 569 zflxs(ji,jj) = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) & 570 & *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp) 571 zdep (ji,jj) = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * & 572 & ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.) 573 zflxs(ji,jj) = zdep(ji,jj) * zflxs(ji,jj) 574 psi (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 546 575 END_2D 547 !548 ! Set psi vertical flux at the surface:549 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope550 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)551 zflxs(:,:) = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &552 & *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)553 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &554 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)555 zflxs(:,:) = zdep(:,:) * zflxs(:,:)556 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)557 576 ! 558 577 END SELECT … … 569 588 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 570 589 ! ! Balance between the production and the dissipation terms 571 DO_2D( 0, 0, 0, 0)590 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 572 591 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 573 592 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 588 607 CASE ( 1 ) ! Neumman boundary condition 589 608 ! 590 DO_2D( 0, 0, 0, 0)609 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 591 610 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 592 611 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 616 635 ! ---------------- 617 636 ! 618 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1637 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 619 638 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 620 639 END_3D 621 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1640 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 622 641 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 623 642 END_3D 624 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk643 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 625 644 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 626 645 END_3D … … 632 651 ! 633 652 CASE( 0 ) ! k-kl (Mellor-Yamada) 634 DO_3D( 0, 0, 0, 0, 1, jpkm1 )653 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 635 654 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 636 655 END_3D 637 656 ! 638 657 CASE( 1 ) ! k-eps 639 DO_3D( 0, 0, 0, 0, 1, jpkm1 )658 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 640 659 eps(ji,jj,jk) = psi(ji,jj,jk) 641 660 END_3D 642 661 ! 643 662 CASE( 2 ) ! k-w 644 DO_3D( 0, 0, 0, 0, 1, jpkm1 )663 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 645 664 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 646 665 END_3D … … 650 669 zex1 = ( 1.5_wp + rmm/rnn ) 651 670 zex2 = -1._wp / rnn 652 DO_3D( 0, 0, 0, 0, 1, jpkm1 )671 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 653 672 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 654 673 END_3D … … 658 677 ! Limit dissipation rate under stable stratification 659 678 ! -------------------------------------------------- 660 DO_3D ( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time679 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 661 680 ! limitation 662 681 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 663 682 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 664 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 665 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 666 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 667 END_3D 683 END_3D 684 IF( ln_length_lim ) THEN ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 685 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 686 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 687 hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 688 END_3D 689 ENDIF 668 690 669 691 ! … … 674 696 ! 675 697 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 676 DO_3D( 0, 0, 0, 0, 2, jpkm1 )698 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 677 699 ! zcof = l²/q² 678 700 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) … … 691 713 ! 692 714 CASE ( 2, 3 ) ! Canuto stability functions 693 DO_3D( 0, 0, 0, 0, 2, jpkm1 )715 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 694 716 ! zcof = l²/q² 695 717 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) … … 723 745 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 724 746 zstm(:,:,jpk) = 0. 725 DO_2D( 0, 0, 0, 0) ! update bottom with good values747 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! update bottom with good values 726 748 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 727 749 END_2D 728 750 729 zstt(:,:, 1) = wmask( :,:, 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)730 zstt(:,:,jpk) = wmask( :,:,jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)751 zstt(:,:, 1) = wmask(A2D(nn_hls), 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 752 zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 731 753 732 754 !!gm should be done for ISF (top boundary cond.) … … 738 760 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 739 761 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 740 DO_3D ( 0, 0, 0, 0, 1, jpk )762 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 741 763 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 742 764 zavt = zsqen * zstt(ji,jj,jk) … … 745 767 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 746 768 END_3D 747 p_avt( :,:,1) = 0._wp769 p_avt(A2D(nn_hls),1) = 0._wp 748 770 ! 749 771 IF(sn_cfctl%l_prtctl) THEN -
NEMO/trunk/src/OCE/ZDF/zdfiwm.F90
r13497 r14834 125 125 ! 126 126 INTEGER :: ji, jj, jk ! dummy loop indices 127 REAL(wp) :: zztmp, ztmp1, ztmp2 ! scalar workspace 128 REAL(wp), DIMENSION(jpi,jpj) :: zfact ! Used for vertical structure 129 REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth 130 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwkb ! WKB-stretched height above bottom 131 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zweight ! Weight for high mode vertical distribution 132 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_t ! Molecular kinematic viscosity (T grid) 133 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_w ! Molecular kinematic viscosity (W grid) 134 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zReb ! Turbulence intensity parameter 135 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zemx_iwm ! local energy density available for mixing (W/kg) 136 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 137 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_wave ! Internal wave-induced diffusivity 127 REAL(wp), SAVE :: zztmp 128 REAL(wp) :: ztmp1, ztmp2 ! scalar workspace 129 REAL(wp), DIMENSION(A2D(nn_hls)) :: zfact ! Used for vertical structure 130 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhdep ! Ocean depth 131 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwkb ! WKB-stretched height above bottom 132 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zweight ! Weight for high mode vertical distribution 133 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_t ! Molecular kinematic viscosity (T grid) 134 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_w ! Molecular kinematic viscosity (W grid) 135 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zReb ! Turbulence intensity parameter 136 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zemx_iwm ! local energy density available for mixing (W/kg) 137 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 138 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zav_wave ! Internal wave-induced diffusivity 138 139 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d ! 3D workspace used for iom_put 139 140 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D - - - - … … 143 144 ! Set to zero the 1st and last vertical levels of appropriate variables 144 145 IF( iom_use("emix_iwm") ) THEN 145 DO_2D( 0, 0, 0, 0)146 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 146 147 zemx_iwm (ji,jj,1) = 0._wp ; zemx_iwm (ji,jj,jpk) = 0._wp 147 148 END_2D 148 149 ENDIF 149 150 IF( iom_use("av_ratio") ) THEN 150 DO_2D( 0, 0, 0, 0)151 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 151 152 zav_ratio(ji,jj,1) = 0._wp ; zav_ratio(ji,jj,jpk) = 0._wp 152 153 END_2D 153 154 ENDIF 154 155 IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN 155 DO_2D( 0, 0, 0, 0)156 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 156 157 zav_wave (ji,jj,1) = 0._wp ; zav_wave (ji,jj,jpk) = 0._wp 157 158 END_2D … … 164 165 ! !* Critical slope mixing: distribute energy over the time-varying ocean depth, 165 166 ! using an exponential decay from the seafloor. 166 DO_2D( 0, 0, 0, 0) ! part independent of the level167 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level 167 168 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 168 169 zfact(ji,jj) = rho0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) … … 170 171 END_2D 171 172 !!gm gde3w ==>>> check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 172 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! complete with the level-dependent part173 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 173 174 IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization 174 175 zemx_iwm(ji,jj,jk) = 0._wp … … 190 191 CASE ( 1 ) ! Dissipation scales as N (recommended) 191 192 ! 192 DO_2D( 0, 0, 0, 0)193 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 193 194 zfact(ji,jj) = 0._wp 194 195 END_2D 195 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! part independent of the level196 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level 196 197 zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk) 197 198 END_3D 198 199 ! 199 DO_2D( 0, 0, 0, 0)200 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 200 201 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 201 202 END_2D 202 203 ! 203 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! complete with the level-dependent part204 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 204 205 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk) 205 206 END_3D … … 207 208 CASE ( 2 ) ! Dissipation scales as N^2 208 209 ! 209 DO_2D( 0, 0, 0, 0)210 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 210 211 zfact(ji,jj) = 0._wp 211 212 END_2D 212 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! part independent of the level213 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level 213 214 zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 214 215 END_3D 215 216 ! 216 DO_2D( 0, 0, 0, 0)217 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 217 218 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 218 219 END_2D 219 220 ! 220 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! complete with the level-dependent part221 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 221 222 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 222 223 END_3D … … 227 228 ! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 228 229 ! 229 DO_2D( 0, 0, 0, 0)230 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 230 231 zwkb(ji,jj,1) = 0._wp 231 232 END_2D 232 DO_3D( 0, 0, 0, 0, 2, jpkm1 )233 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 233 234 zwkb(ji,jj,jk) = zwkb(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk) 234 235 END_3D 235 DO_2D( 0, 0, 0, 0)236 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 236 237 zfact(ji,jj) = zwkb(ji,jj,jpkm1) 237 238 END_2D 238 239 ! 239 DO_3D( 0, 0, 0, 0, 2, jpkm1 )240 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 240 241 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 241 242 & * wmask(ji,jj,jk) / zfact(ji,jj) 242 243 END_3D 243 DO_2D( 0, 0, 0, 0)244 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 244 245 zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1) 245 246 END_2D 246 247 ! 247 DO_3D( 0, 0, 0, 0, 2, jpkm1 )248 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 248 249 IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization: EXP coast a lot 249 250 zweight(ji,jj,jk) = 0._wp … … 254 255 END_3D 255 256 ! 256 DO_2D( 0, 0, 0, 0)257 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 257 258 zfact(ji,jj) = 0._wp 258 259 END_2D 259 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! part independent of the level260 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level 260 261 zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk) 261 262 END_3D 262 263 ! 263 DO_2D( 0, 0, 0, 0)264 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 264 265 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 265 266 END_2D 266 267 ! 267 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! complete with the level-dependent part268 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 268 269 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk) & 269 270 & / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) … … 273 274 !!gm this is to be replaced by just a constant value znu=1.e-6 m2/s 274 275 ! Calculate molecular kinematic viscosity 275 DO_3D( 0, 0, 0, 0, 1, jpkm1 )276 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 276 277 znu_t(ji,jj,jk) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm) & 277 278 & + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) & 278 279 & + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm) ) * tmask(ji,jj,jk) * r1_rho0 279 280 END_3D 280 DO_3D( 0, 0, 0, 0, 2, jpkm1 )281 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 281 282 znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk) 282 283 END_3D … … 284 285 ! 285 286 ! Calculate turbulence intensity parameter Reb 286 DO_3D( 0, 0, 0, 0, 2, jpkm1 )287 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 287 288 zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) ) 288 289 END_3D 289 290 ! 290 291 ! Define internal wave-induced diffusivity 291 DO_3D( 0, 0, 0, 0, 2, jpkm1 )292 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 292 293 zav_wave(ji,jj,jk) = znu_w(ji,jj,jk) * zReb(ji,jj,jk) * r1_6 ! This corresponds to a constant mixing efficiency of 1/6 293 294 END_3D 294 295 ! 295 296 IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the 296 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes297 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 297 298 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 298 299 zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) … … 303 304 ENDIF 304 305 ! 305 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Bound diffusivity by molecular value and 100 cm2/s306 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Bound diffusivity by molecular value and 100 cm2/s 306 307 zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk) 307 308 END_3D 308 309 ! 309 310 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave 310 zztmp = 0._wp311 IF( .NOT. l_istiled .OR. ntile == 1 ) zztmp = 0._wp ! Do only on the first tile 311 312 !!gm used of glosum 3D.... 312 313 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) … … 314 315 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 315 316 END_3D 316 CALL mpp_sum( 'zdfiwm', zztmp ) 317 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 318 ! 319 IF(lwp) THEN 320 WRITE(numout,*) 321 WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 322 WRITE(numout,*) '~~~~~~~ ' 323 WRITE(numout,*) 324 WRITE(numout,*) ' Total power consumption by av_wave = ', zztmp * 1.e-12_wp, 'TW' 317 318 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 319 CALL mpp_sum( 'zdfiwm', zztmp ) 320 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 321 ! 322 IF(lwp) THEN 323 WRITE(numout,*) 324 WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 325 WRITE(numout,*) '~~~~~~~ ' 326 WRITE(numout,*) 327 WRITE(numout,*) ' Total power consumption by av_wave = ', zztmp * 1.e-12_wp, 'TW' 328 ENDIF 325 329 ENDIF 326 330 ENDIF … … 332 336 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature 333 337 ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 334 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb338 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb 335 339 ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 336 340 IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN … … 341 345 END_3D 342 346 CALL iom_put( "av_ratio", zav_ratio ) 343 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing347 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing 344 348 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 345 349 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) … … 348 352 ! 349 353 ELSE !* update momentum & tracer diffusivity with wave-driven mixing 350 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )354 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 351 355 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 352 356 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) … … 361 365 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 362 366 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 363 ALLOCATE( z2d( jpi,jpj) , z3d(jpi,jpj,jpk) )367 ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) ) 364 368 ! Initialisation for iom_put 365 369 DO_2D( 0, 0, 0, 0 ) 366 370 z3d(ji,jj,1) = 0._wp ; z3d(ji,jj,jpk) = 0._wp 367 371 END_2D 368 z3d( 1:nn_hls,:,:) = 0._wp ; z3d(:, 1:nn_hls,:) = 0._wp369 z3d(jpi-nn_hls+1:jpi ,:,:) = 0._wp ; z3d(:,jpj-nn_hls+1: jpj,:) = 0._wp370 z2d( 1:nn_hls,: ) = 0._wp ; z2d(:, 1:nn_hls ) = 0._wp371 z2d(jpi-nn_hls+1:jpi ,: ) = 0._wp ; z2d(:,jpj-nn_hls+1: jpj ) = 0._wp372 372 373 373 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) -
NEMO/trunk/src/OCE/ZDF/zdfmfc.F90
r14433 r14834 96 96 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices 97 97 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 98 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: ztsp ! T/S of the plume99 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: ztse ! T/S at W point100 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zrwp !101 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zrwp2 !102 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zapp !103 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zedmf !104 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zepsT, zepsW !105 ! 106 REAL(wp), DIMENSION( jpi,jpj) :: zustar, zustar2 !107 REAL(wp), DIMENSION( jpi,jpj) :: zuws, zvws, zsws, zfnet !108 REAL(wp), DIMENSION( jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl109 REAL(wp), DIMENSION( jpi,jpj) :: zwpsurf !110 REAL(wp), DIMENSION( jpi,jpj) :: zop0 , zsp0 !111 REAL(wp), DIMENSION( jpi,jpj) :: zrwp_0, zrwp2_0 !112 REAL(wp), DIMENSION( jpi,jpj) :: zapp0 !113 REAL(wp), DIMENSION( jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro114 REAL(wp), DIMENSION( jpi,jpj) :: zhcmo !115 ! 116 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zn2 ! N^2117 REAL(wp), DIMENSION( jpi,jpj,2 ) :: zab, zabm1, zabp ! alpha and beta98 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: ztsp ! T/S of the plume 99 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: ztse ! T/S at W point 100 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp ! 101 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp2 ! 102 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zapp ! 103 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zedmf ! 104 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zepsT, zepsW ! 105 ! 106 REAL(wp), DIMENSION(A2D(nn_hls)) :: zustar, zustar2 ! 107 REAL(wp), DIMENSION(A2D(nn_hls)) :: zuws, zvws, zsws, zfnet ! 108 REAL(wp), DIMENSION(A2D(nn_hls)) :: zfbuo, zrautbm1, zrautb, zraupl 109 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwpsurf ! 110 REAL(wp), DIMENSION(A2D(nn_hls)) :: zop0 , zsp0 ! 111 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrwp_0, zrwp2_0 ! 112 REAL(wp), DIMENSION(A2D(nn_hls)) :: zapp0 ! 113 REAL(wp), DIMENSION(A2D(nn_hls)) :: zphp, zph, zphpm1, zphm1, zNHydro 114 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhcmo ! 115 ! 116 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zn2 ! N^2 117 REAL(wp), DIMENSION(A2D(nn_hls),2 ) :: zab, zabm1, zabp ! alpha and beta 118 118 119 119 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value … … 136 136 zcd = 1._wp 137 137 138 !------------------------------------------------------------------ 139 ! Surface boundary condition 140 !------------------------------------------------------------------ 141 ! surface Stress 142 !-------------------- 143 zuws(:,:) = utau(:,:) * r1_rho0 144 zvws(:,:) = vtau(:,:) * r1_rho0 145 zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:)) 146 zustar(:,:) = SQRT(zustar2(:,:)) 147 148 ! Heat Flux 149 !-------------------- 150 zfnet(:,:) = qns(:,:) + qsr(:,:) 151 zfnet(:,:) = zfnet(:,:) / (rho0 * rcp) 152 153 ! Water Flux 154 !--------------------- 155 zsws(:,:) = emp(:,:) 156 157 !------------------------------------------- 158 ! Initialisation of prognostic variables 159 !------------------------------------------- 160 zrwp (:,:,:) = 0._wp ; zrwp2(:,:,:) = 0._wp ; zedmf(:,:,:) = 0._wp 161 zph (:,:) = 0._wp ; zphm1(:,:) = 0._wp ; zphpm1(:,:) = 0._wp 162 ztsp(:,:,:,:)= 0._wp 163 164 ! Tracers inside plume (ztsp) and environment (ztse) 165 ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) 166 ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 167 ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) 168 ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 138 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 139 !------------------------------------------------------------------ 140 ! Surface boundary condition 141 !------------------------------------------------------------------ 142 ! surface Stress 143 !-------------------- 144 zuws(ji,jj) = utau(ji,jj) * r1_rho0 145 zvws(ji,jj) = vtau(ji,jj) * r1_rho0 146 zustar2(ji,jj) = SQRT(zuws(ji,jj)*zuws(ji,jj)+zvws(ji,jj)*zvws(ji,jj)) 147 zustar(ji,jj) = SQRT(zustar2(ji,jj)) 148 149 ! Heat Flux 150 !-------------------- 151 zfnet(ji,jj) = qns(ji,jj) + qsr(ji,jj) 152 zfnet(ji,jj) = zfnet(ji,jj) / (rho0 * rcp) 153 154 ! Water Flux 155 !--------------------- 156 zsws(ji,jj) = emp(ji,jj) 157 158 !------------------------------------------- 159 ! Initialisation of prognostic variables 160 !------------------------------------------- 161 zrwp (ji,jj,:) = 0._wp ; zrwp2(ji,jj,:) = 0._wp ; zedmf(ji,jj,:) = 0._wp 162 zph (ji,jj) = 0._wp ; zphm1(ji,jj) = 0._wp ; zphpm1(ji,jj) = 0._wp 163 ztsp(ji,jj,:,:)= 0._wp 164 165 ! Tracers inside plume (ztsp) and environment (ztse) 166 ztsp(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1) 167 ztsp(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1) 168 ztse(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1) 169 ztse(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1) 170 END_2D 169 171 170 172 CALL eos( ztse(:,:,1,:) , zrautb(:,:) ) … … 174 176 ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain) 175 177 !------------------------------------------- 176 zhcmo(:,:) = e3t( :,:,1,Kmm)178 zhcmo(:,:) = e3t(A1Di(nn_hls),A1Dj(nn_hls),1,Kmm) 177 179 zfbuo(:,:) = 0._wp 178 180 WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:) = & 179 & grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) 181 & grav * ( 2.e-4_wp *zfnet(:,:) & 182 & - 7.6E-4_wp*pts(A2D(nn_hls),1,jp_sal,Kmm) & 183 & * zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) 180 184 181 185 zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:)) … … 211 215 CALL eos( ztsp(:,:,jk-1,: ) , zraupl(:,:) ) 212 216 213 zphm1(:,:) = zphm1(:,:) + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm) 214 zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:) * e3t(:,:,jk-1, Kmm) 215 zph(:,:) = zphm1(:,:) + grav * zrautb(:,:) * e3t(:,:,jk , Kmm) 216 zph(:,:) = MAX( zph(:,:), zepsilon) 217 DO_2D( 0, 0, 0, 0 ) 218 zphm1(ji,jj) = zphm1(ji,jj) + grav * zrautbm1(ji,jj) * e3t(ji,jj,jk-1, Kmm) 219 zphpm1(ji,jj) = zphpm1(ji,jj) + grav * zraupl(ji,jj) * e3t(ji,jj,jk-1, Kmm) 220 zph(ji,jj) = zphm1(ji,jj) + grav * zrautb(ji,jj) * e3t(ji,jj,jk , Kmm) 221 zph(ji,jj) = MAX( zph(ji,jj), zepsilon) 222 END_2D 217 223 218 224 WHERE(zrautbm1 .NE. 0.) zfbuo(:,:) = grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:) … … 322 328 323 329 ! Compute Mass Flux on T-point 324 DO jk=1,jpk-1 325 edmfm(:,:,jk) = (zedmf(:,:,jk+1) + zedmf(:,:,jk) )*0.5_wp 326 END DO 327 edmfm(:,:,jpk) = zedmf(:,:,jpk) 330 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 331 edmfm(ji,jj,jk) = (zedmf(ji,jj,jk+1) + zedmf(ji,jj,jk) )*0.5_wp 332 END_3D 333 DO_2D( 0, 0, 0, 0 ) 334 edmfm(ji,jj,jpk) = zedmf(ji,jj,jpk) 335 END_2D 328 336 329 337 ! Save variable (on T point) … … 338 346 ! Computation of a tridiagonal matrix and right hand side terms of the linear system 339 347 !================================================================================= 340 edmfa(:,:,:) = 0._wp 341 edmfb(:,:,:) = 0._wp 342 edmfc(:,:,:) = 0._wp 343 edmftra(:,:,:,:) = 0._wp 348 DO_3D( 0, 0, 0, 0, 1, jpk ) 349 edmfa(ji,jj,jk) = 0._wp 350 edmfb(ji,jj,jk) = 0._wp 351 edmfc(ji,jj,jk) = 0._wp 352 edmftra(ji,jj,jk,:) = 0._wp 353 END_3D 344 354 345 355 !--------------------------------------------------------------- 346 356 ! Diagonal terms 347 357 !--------------------------------------------------------------- 348 DO jk=1,jpk-1 349 edmfa(:,:,jk) = 0._wp 350 edmfb(:,:,jk) = -edmfm(:,:,jk ) / e3w(:,:,jk+1,Kmm) 351 edmfc(:,:,jk) = edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm) 352 END DO 353 edmfa(:,:,jpk) = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm) 354 edmfb(:,:,jpk) = edmfm(:,:,jpk ) / e3w(:,:,jpk,Kmm) 355 edmfc(:,:,jpk) = 0._wp 358 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 359 edmfa(ji,jj,jk) = 0._wp 360 edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm) 361 edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 362 END_3D 363 DO_2D( 0, 0, 0, 0 ) 364 edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm) 365 edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm) 366 edmfc(ji,jj,jpk) = 0._wp 367 END_2D 356 368 357 369 !--------------------------------------------------------------- 358 370 ! right hand side term for Temperature 359 371 !--------------------------------------------------------------- 360 DO jk=1,jpk-1 361 edmftra(:,:,jk,1) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_tem) / e3w(:,:,jk+1,Kmm) & 362 & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm) 363 END DO 364 edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) & 365 & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_tem) / e3w(:,:,jpk,Kmm) 366 372 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 373 edmftra(ji,jj,jk,1) = - edmfm(ji,jj,jk ) * ztsp(ji,jj,jk ,jp_tem) / e3w(ji,jj,jk+1,Kmm) & 374 & + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_tem) / e3w(ji,jj,jk+1,Kmm) 375 END_3D 376 DO_2D( 0, 0, 0, 0 ) 377 edmftra(ji,jj,jpk,1) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_tem) / e3w(ji,jj,jpk,Kmm) & 378 & + edmfm(ji,jj,jpk ) * ztsp(ji,jj,jpk ,jp_tem) / e3w(ji,jj,jpk,Kmm) 379 END_2D 380 367 381 !--------------------------------------------------------------- 368 382 ! Right hand side term for Salinity 369 383 !--------------------------------------------------------------- 370 DO jk=1,jpk-1 371 edmftra(:,:,jk,2) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_sal) / e3w(:,:,jk+1,Kmm) & 372 & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm) 373 END DO 374 edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) & 375 & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_sal) / e3w(:,:,jpk,Kmm) 376 ! 377 ! 378 CALL lbc_lnk( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.) 384 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 385 edmftra(ji,jj,jk,2) = - edmfm(ji,jj,jk ) * ztsp(ji,jj,jk ,jp_sal) / e3w(ji,jj,jk+1,Kmm) & 386 & + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_sal) / e3w(ji,jj,jk+1,Kmm) 387 END_3D 388 DO_2D( 0, 0, 0, 0 ) 389 edmftra(ji,jj,jpk,2) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_sal) / e3w(ji,jj,jpk,Kmm) & 390 & + edmfm(ji,jj,jpk ) * ztsp(ji,jj,jpk ,jp_sal) / e3w(ji,jj,jpk,Kmm) 391 END_2D 379 392 ! 380 393 END SUBROUTINE tra_mfc … … 383 396 SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa ) 384 397 385 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms386 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step387 INTEGER , INTENT(in ) :: Kaa ! ocean time level indices398 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms 399 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 400 INTEGER , INTENT(in ) :: Kaa ! ocean time level indices 388 401 389 402 INTEGER :: ji, jj, jk ! dummy loop arguments -
NEMO/trunk/src/OCE/ZDF/zdfmxl.F90
r13497 r14834 26 26 PRIVATE 27 27 28 PUBLIC zdf_mxl ! called by zdfphy.F9028 PUBLIC zdf_mxl, zdf_mxl_turb ! called by zdfphy.F90 29 29 30 30 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP) … … 41 41 !!---------------------------------------------------------------------- 42 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 43 !! $Id$ 43 !! $Id$ 44 44 !! Software governed by the CeCILL license (see ./LICENSE) 45 45 !!---------------------------------------------------------------------- … … 65 65 !! *** ROUTINE zdfmxl *** 66 66 !! 67 !! ** Purpose : Compute the turbocline depth and the mixed layer depth 68 !! with density criteria. 67 !! ** Purpose : Compute the mixed layer depth with density criteria. 69 68 !! 70 69 !! ** Method : The mixed layer depth is the shallowest W depth with 71 70 !! the density of the corresponding T point (just bellow) bellow a 72 71 !! given value defined locally as rho(10m) + rho_c 73 !! The turbocline depth is the depth at which the vertical74 !! eddy diffusivity coefficient (resulting from the vertical physics75 !! alone, not the isopycnal part, see trazdf.F) fall below a given76 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default)77 72 !! 78 !! ** Action : nmln, hml d, hmlp, hmlpt73 !! ** Action : nmln, hmlp, hmlpt 79 74 !!---------------------------------------------------------------------- 80 75 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 82 77 ! 83 78 INTEGER :: ji, jj, jk ! dummy loop indices 84 INTEGER :: iik n, iiki, ikt! local integer79 INTEGER :: iik, ikt ! local integer 85 80 REAL(wp) :: zN2_c ! local scalar 86 INTEGER, DIMENSION(jpi,jpj) :: imld ! 2D workspace87 81 !!---------------------------------------------------------------------- 88 82 ! 89 IF( kt == nit000 ) THEN 90 IF(lwp) WRITE(numout,*) 91 IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 92 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 93 ! ! allocate zdfmxl arrays 94 IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' ) 83 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 84 IF( kt == nit000 ) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 87 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 88 ! ! allocate zdfmxl arrays 89 IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' ) 90 ENDIF 95 91 ENDIF 96 92 ! 97 93 ! w-level of the mixing and mixed layers 98 nmln(:,:) = nlb10 ! Initialization to the number of w ocean point 99 hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 94 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 95 nmln(ji,jj) = nlb10 ! Initialization to the number of w ocean point 96 hmlp(ji,jj) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 97 END_2D 100 98 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 101 DO_3D ( 1, 1, 1, 1, nlb10, jpkm1 ) ! Mixed layer level: w-level99 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) ! Mixed layer level: w-level 102 100 ikt = mbkt(ji,jj) 103 101 hmlp(ji,jj) = & … … 105 103 IF( hmlp(ji,jj) < zN2_c ) nmln(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 106 104 END_3D 107 ! 108 ! w-level of the turbocline and mixing layer (iom_use) 109 imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 110 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 111 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 112 END_3D 113 ! depth of the mixing and mixed layers 114 DO_2D( 1, 1, 1, 1 ) 115 iiki = imld(ji,jj) 116 iikn = nmln(ji,jj) 117 hmld (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth 118 hmlp (ji,jj) = gdepw(ji,jj,iikn ,Kmm) * ssmask(ji,jj) ! Mixed layer depth 119 hmlpt(ji,jj) = gdept(ji,jj,iikn-1,Kmm) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 105 ! depth of the mixed layer 106 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 107 iik = nmln(ji,jj) 108 hmlp (ji,jj) = gdepw(ji,jj,iik ,Kmm) * ssmask(ji,jj) ! Mixed layer depth 109 hmlpt(ji,jj) = gdept(ji,jj,iik-1,Kmm) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 120 110 END_2D 121 111 ! 122 IF( .NOT.l_offline ) THEN 123 IF( iom_use("mldr10_1") ) THEN 124 IF( ln_isfcav ) THEN ; CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 125 ELSE ; CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 126 END IF 112 IF( .NOT.l_offline .AND. iom_use("mldr10_1") ) THEN 113 IF( ln_isfcav ) THEN ; CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 114 ELSE ; CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 127 115 END IF 128 IF( iom_use("mldkz5") ) THEN129 IF( ln_isfcav ) THEN ; CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness130 ELSE ; CALL iom_put( "mldkz5" , hmld ) ! turbocline depth131 END IF132 ENDIF133 116 ENDIF 134 117 ! … … 137 120 END SUBROUTINE zdf_mxl 138 121 122 123 SUBROUTINE zdf_mxl_turb( kt, Kmm ) 124 !!---------------------------------------------------------------------- 125 !! *** ROUTINE zdf_mxl_turb *** 126 !! 127 !! ** Purpose : Compute the turbocline depth. 128 !! 129 !! ** Method : The turbocline depth is the depth at which the vertical 130 !! eddy diffusivity coefficient (resulting from the vertical physics 131 !! alone, not the isopycnal part, see trazdf.F) fall below a given 132 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default) 133 !! 134 !! ** Action : hmld 135 !!---------------------------------------------------------------------- 136 INTEGER, INTENT(in) :: kt ! ocean time-step index 137 INTEGER, INTENT(in) :: Kmm ! ocean time level index 138 ! 139 INTEGER :: ji, jj, jk ! dummy loop indices 140 INTEGER :: iik ! local integer 141 INTEGER, DIMENSION(A2D(nn_hls)) :: imld ! 2D workspace 142 !!---------------------------------------------------------------------- 143 ! 144 ! w-level of the turbocline and mixing layer (iom_use) 145 imld(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point 146 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 147 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 148 END_3D 149 ! depth of the mixing layer 150 DO_2D_OVR( 1, 1, 1, 1 ) 151 iik = imld(ji,jj) 152 hmld (ji,jj) = gdepw(ji,jj,iik ,Kmm) * ssmask(ji,jj) ! Turbocline depth 153 END_2D 154 ! 155 IF( .NOT.l_offline .AND. iom_use("mldkz5") ) THEN 156 IF( ln_isfcav ) THEN ; CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness 157 ELSE ; CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 158 END IF 159 ENDIF 160 ! 161 END SUBROUTINE zdf_mxl_turb 139 162 !!====================================================================== 140 163 END MODULE zdfmxl -
NEMO/trunk/src/OCE/ZDF/zdfphy.F90
r14433 r14834 12 12 !!---------------------------------------------------------------------- 13 13 USE oce ! ocean dynamics and tracers variables 14 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled) 15 USE domtile 14 16 USE zdf_oce ! vertical physics: shared variables 15 17 USE zdfdrg ! vertical physics: top/bottom drag coef. … … 54 56 INTEGER, PARAMETER :: np_OSM = 5 ! OSMOSIS-OBL closure scheme for Kz 55 57 56 LOGICAL :: l_zdfsh2 ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 57 58 LOGICAL, PUBLIC :: l_zdfsh2 ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 59 60 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm_k_n !: "Now" avm_k used for calculation of zsh2 with tiling 61 62 !! * Substitutions 63 # include "do_loop_substitute.h90" 58 64 !!---------------------------------------------------------------------- 59 65 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 180 186 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 181 187 IF( lk_top .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 188 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled) 189 IF( ln_tile .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis does not yet work with tiling' ) 182 190 IF( lk_top .AND. ln_zdfmfc ) CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' ) 183 191 IF(lwp) THEN … … 210 218 ENDIF 211 219 ! ! shear production term flag 212 IF( ln_zdfcst ) THEN ; l_zdfsh2 = .FALSE. 213 ELSE ; l_zdfsh2 = .TRUE. 214 ENDIF 220 IF( ln_zdfcst .OR. ln_zdfosm ) THEN ; l_zdfsh2 = .FALSE. 221 ELSE ; l_zdfsh2 = .TRUE. 222 ENDIF 223 IF( ln_tile .AND. l_zdfsh2 ) ALLOCATE( avm_k_n(jpi,jpj,jpk) ) 215 224 ! !== Mass Flux Convectiive algorithm ==! 216 225 IF( ln_zdfmfc ) CALL zdf_mfc_init ! Convection computed with eddy diffusivity mass flux … … 246 255 ! 247 256 INTEGER :: ji, jj, jk ! dummy loop indice 248 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! shear production 257 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zsh2 ! shear production 258 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled) 259 LOGICAL :: lskip 249 260 !! --------------------------------------------------------------------- 250 261 ! 251 262 IF( ln_timing ) CALL timing_start('zdf_phy') 263 264 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled) 265 lskip = .FALSE. 266 267 IF( ln_tile .AND. nzdf_phy == np_OSM ) THEN 268 IF( ntile == 1 ) THEN 269 CALL dom_tile_stop( ldhold=.TRUE. ) 270 ELSE 271 lskip = .TRUE. 272 ENDIF 273 ENDIF 252 274 ! 253 275 IF( l_zdfdrg ) THEN !== update top/bottom drag ==! (non-linear cases) … … 267 289 IF ( ln_drgice_imp) THEN 268 290 IF ( ln_isfcav ) THEN 269 rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:) 291 DO_2D_OVR( 1, 1, 1, 1 ) 292 rCdU_top(ji,jj) = rCdU_top(ji,jj) + ssmask(ji,jj) * tmask(ji,jj,1) * rCdU_ice(ji,jj) 293 END_2D 270 294 ELSE 271 rCdU_top(:,:) = rCdU_ice(:,:) 295 DO_2D_OVR( 1, 1, 1, 1 ) 296 rCdU_top(ji,jj) = rCdU_ice(ji,jj) 297 END_2D 272 298 ENDIF 273 299 ENDIF 274 300 #endif 275 301 ! 276 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 277 ! 278 IF( l_zdfsh2 ) & !* shear production at w-points (energy conserving form) 279 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 280 & zsh2 ) ! ==>> out : shear production 281 ! 282 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 283 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 284 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 285 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 286 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 287 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 288 ! ! avt_k and avm_k set one for all at initialisation phase 302 CALL zdf_mxl( kt, Kmm ) !* mixed layer depth, and level 303 304 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled) 305 IF( .NOT. lskip ) THEN 306 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 307 ! 308 ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2. 309 IF( l_zdfsh2 ) THEN !* shear production at w-points (energy conserving form) 310 IF( ln_tile ) THEN 311 IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:) ! Preserve "now" avm_k for calculation of zsh2 312 CALL zdf_sh2( Kbb, Kmm, avm_k_n, & ! <<== in 313 & zsh2 ) ! ==>> out : shear production 314 ELSE 315 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 316 & zsh2 ) ! ==>> out : shear production 317 ENDIF 318 ENDIF 319 ! 320 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 321 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 322 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 323 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 324 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 325 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 326 ! ! avt_k and avm_k set one for all at initialisation phase 289 327 !!gm avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 290 328 !!gm avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 291 END SELECT 329 END SELECT 330 331 IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. ) 332 ENDIF 292 333 ! 293 334 ! !== ocean Kz ==! (avt, avs, avm) 294 335 ! 295 336 ! !* start from turbulent closure values 296 avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1) 297 avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1) 337 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 338 avt(ji,jj,jk) = avt_k(ji,jj,jk) 339 avm(ji,jj,jk) = avm_k(ji,jj,jk) 340 END_3D 298 341 ! 299 342 IF( ln_rnf_mouth ) THEN !* increase diffusivity at rivers mouths 300 DO jk = 2, nkrnf301 avt( :,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)302 END DO343 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, nkrnf ) 344 avt(ji,jj,jk) = avt(ji,jj,jk) + 2._wp * rn_avt_rnf * rnfmsk(ji,jj) * wmask(ji,jj,jk) 345 END_3D 303 346 ENDIF 304 347 ! … … 309 352 CALL zdf_ddm( kt, Kmm, avm, avt, avs ) 310 353 ELSE ! same mixing on all tracers 311 avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 354 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 355 avs(ji,jj,jk) = avt(ji,jj,jk) 356 END_3D 312 357 ENDIF 313 358 ! … … 318 363 #if defined key_agrif 319 364 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 320 IF( l_zdfsh2 ) CALL Agrif_avm 365 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 366 IF( l_zdfsh2 ) CALL Agrif_avm 367 ENDIF 321 368 #endif 322 369 323 370 ! !* Lateral boundary conditions (sign unchanged) 324 IF( l_zdfsh2 ) THEN 325 CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, & 326 & avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 327 ELSE 328 CALL lbc_lnk( 'zdfphy', avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 329 ENDIF 330 ! 331 IF( l_zdfdrg ) THEN ! drag have been updated (non-linear cases) 332 IF( ln_isfcav ) THEN ; CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp ) ! top & bot drag 333 ELSE ; CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp ) ! bottom drag only 334 ENDIF 335 ENDIF 336 ! 337 CALL zdf_mxl( kt, Kmm ) !* mixed layer depth, and level 338 ! 339 IF( lrst_oce ) THEN !* write TKE, GLS or RIC fields in the restart file 340 IF( ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 341 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 342 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 343 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 371 IF(nn_hls==1) THEN 372 IF( l_zdfsh2 ) THEN 373 CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, & 374 & avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 375 ELSE 376 CALL lbc_lnk( 'zdfphy', avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 377 ENDIF 378 ! 379 IF( l_zdfdrg ) THEN ! drag have been updated (non-linear cases) 380 IF( ln_isfcav ) THEN ; CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp ) ! top & bot drag 381 ELSE ; CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp ) ! bottom drag only 382 ENDIF 383 ENDIF 384 ENDIF 385 ! 386 CALL zdf_mxl_turb( kt, Kmm ) !* turbocline depth 387 ! 388 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 389 IF( lrst_oce ) THEN !* write TKE, GLS or RIC fields in the restart file 390 IF( ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 391 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 392 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 393 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 394 ENDIF 344 395 ENDIF 345 396 ! -
NEMO/trunk/src/OCE/ZDF/zdfric.F90
r14072 r14834 145 145 !! PFJ Lermusiaux 2001. 146 146 !!---------------------------------------------------------------------- 147 INTEGER , INTENT(in ) :: kt ! ocean time-step148 INTEGER , INTENT(in ) :: Kmm ! ocean time level index149 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term150 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)147 INTEGER , INTENT(in ) :: kt ! ocean time-step 148 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 149 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 150 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 151 151 !! 152 152 INTEGER :: ji, jj, jk ! dummy loop indices 153 153 REAL(wp) :: zcfRi, zav, zustar, zhek ! local scalars 154 REAL(wp), DIMENSION( jpi,jpj) :: zh_ekm ! 2D workspace154 REAL(wp), DIMENSION(A2D(nn_hls)) :: zh_ekm ! 2D workspace 155 155 !!---------------------------------------------------------------------- 156 156 ! 157 157 ! !== avm and avt = F(Richardson number) ==! 158 DO_3D ( 1, 0, 1, 0, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri)158 DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri) 159 159 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) ) 160 160 zav = rn_avmri * zcfRi**nn_ric … … 169 169 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 170 170 ! 171 DO_2D( 0, 0, 0, 0 ) !* Ekman depth171 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 172 172 zustar = SQRT( taum(ji,jj) * r1_rho0 ) 173 173 zhek = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 174 174 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zhek , rn_mldmax ) ) ! set allowed range 175 175 END_2D 176 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* minimum mixing coeff. within the Ekman layer176 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* minimum mixing coeff. within the Ekman layer 177 177 IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN 178 178 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) -
NEMO/trunk/src/OCE/ZDF/zdfsh2.F90
r14072 r14834 55 55 !! References : Bruchard, OM 2002 56 56 !! --------------------------------------------------------------------- 57 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices58 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points)59 REAL(wp), DIMENSION( :,:,:) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points)57 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 58 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points) 59 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points) 60 60 ! 61 61 INTEGER :: ji, jj, jk ! dummy loop arguments 62 REAL(wp), DIMENSION( jpi,jpj) :: zsh2u, zsh2v ! 2D workspace62 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsh2u, zsh2v ! 2D workspace 63 63 !!-------------------------------------------------------------------- 64 64 ! 65 65 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 66 66 IF ( cpl_sdrftx .AND. ln_stshear ) THEN ! Surface Stokes Drift available ===>>> shear + stokes drift contibution 67 DO_2D( 1, 0, 1, 0)67 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 68 68 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 69 69 & * ( uu (ji,jj,jk-1,Kmm) - uu (ji,jj,jk,Kmm) & … … 78 78 END_2D 79 79 ELSE 80 DO_2D( 1, 0, 1, 0) !* 2 x shear production at uw- and vw-points (energy conserving form)80 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) !* 2 x shear production at uw- and vw-points (energy conserving form) 81 81 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 82 82 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & … … 91 91 END_2D 92 92 ENDIF 93 DO_2D( 0, 0, 0, 0) !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)93 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 94 94 p_sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) & 95 95 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) ) -
NEMO/trunk/src/OCE/ZDF/zdfswm.F90
r13295 r14834 63 63 ! 64 64 zcoef = 1._wp * 0.353553_wp 65 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )65 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 66 66 zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw(ji,jj,jk,Kmm) ) * wmask(ji,jj,jk) 67 67 ! -
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r14072 r14834 168 168 !! Bruchard OM 2002 169 169 !!---------------------------------------------------------------------- 170 INTEGER , INTENT(in ) :: kt ! ocean time step171 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices172 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term173 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)170 INTEGER , INTENT(in ) :: kt ! ocean time step 171 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 172 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 173 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 174 174 !!---------------------------------------------------------------------- 175 175 ! … … 201 201 USE zdf_oce , ONLY : en ! ocean vertical physics 202 202 !! 203 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices204 REAL(wp), DIMENSION( :,:,:) , INTENT(in ) :: p_sh2 ! shear production term205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points)203 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 204 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in ) :: p_sh2 ! shear production term 205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 206 206 ! 207 207 INTEGER :: ji, jj, jk ! dummy loop arguments … … 216 216 REAL(wp) :: zzd_up, zzd_lw ! - - 217 217 REAL(wp) :: ztaui, ztauj, z1_norm 218 INTEGER , DIMENSION( jpi,jpj) :: imlc219 REAL(wp), DIMENSION( jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2220 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw218 INTEGER , DIMENSION(A2D(nn_hls)) :: imlc 219 REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra, zhlc, zus3, zWlc2 220 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpelc, zdiag, zd_up, zd_lw 221 221 !!-------------------------------------------------------------------- 222 222 ! … … 232 232 SELECT CASE ( nn_eice ) 233 233 CASE( 0 ) ; zice_fra(:,:) = 0._wp 234 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i( :,:) * 10._wp )235 CASE( 2 ) ; zice_fra(:,:) = fr_i( :,:)236 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i( :,:) , 1._wp )234 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 235 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 236 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 237 237 END SELECT 238 238 ! … … 241 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 242 242 ! 243 DO_2D ( 0, 0, 0, 0)243 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 244 244 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) 245 245 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) … … 258 258 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 259 259 ! 260 DO_2D ( 0, 0, 0, 0) ! bottom friction260 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction 261 261 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 262 262 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) … … 267 267 END_2D 268 268 IF( ln_isfcav ) THEN 269 DO_2D ( 0, 0, 0, 0) ! top friction269 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 270 270 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 271 271 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) … … 294 294 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 295 295 !!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect ! 296 DO_2D( 0, 0, 0, 0)296 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 297 297 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) 298 298 zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) … … 301 301 ! Projection of Stokes drift in the wind stress direction 302 302 ! 303 DO_2D( 0, 0, 0, 0)303 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 304 304 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 305 305 ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) … … 307 307 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 308 308 END_2D 309 CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. )310 !311 309 ELSE ! Surface Stokes drift deduced from surface stress 312 310 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) … … 315 313 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 316 314 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 317 DO_2D( 1, 1, 1,1 )315 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 318 316 zWlc2(ji,jj) = zcof * taum(ji,jj) 319 317 END_2D … … 323 321 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 324 322 ! !- LHS of Eq.47 325 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 326 DO jk = 2, jpk 327 zpelc(:,:,jk) = zpelc(:,:,jk-1) + & 328 & MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 329 END DO 323 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 324 zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm) 325 END_2D 326 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk ) 327 zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + & 328 & MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 329 END_3D 330 330 ! 331 331 ! !- compare LHS to RHS of Eq.47 332 imlc(:,:) = mbkt( :,:) + 1 ! Initialization to the number of w ocean point (=2 over land)333 DO_3DS( 1, 1, 1,1, jpkm1, 2, -1 )332 imlc(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point (=2 over land) 333 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) 334 334 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 335 335 END_3D 336 336 ! ! finite LC depth 337 DO_2D( 1, 1, 1,1 )337 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 338 338 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 339 339 END_2D 340 340 ! 341 341 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 342 DO_2D( 0, 0, 0, 0)342 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 343 343 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 344 344 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 345 345 END_2D 346 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en346 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 347 347 IF ( zus3(ji,jj) /= 0._wp ) THEN 348 348 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN … … 365 365 ! 366 366 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 367 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )367 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 368 368 ! ! local Richardson number 369 369 IF (rn2b(ji,jj,jk) <= 0.0_wp) then … … 377 377 ENDIF 378 378 ! 379 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en379 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en 380 380 zcof = zfact1 * tmask(ji,jj,jk) 381 381 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical … … 406 406 407 407 CASE ( 0 ) ! Dirichlet BC 408 DO_2D ( 0, 0, 0, 0) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0)408 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 409 409 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 410 410 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) … … 413 413 414 414 CASE ( 1 ) ! Neumann BC 415 DO_2D ( 0, 0, 0, 0)415 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 416 416 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 417 417 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) … … 427 427 ! 428 428 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 429 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1429 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 430 430 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 431 431 END_3D … … 434 434 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 435 435 ! END_2D 436 DO_3D( 0, 0, 0, 0, 2, jpkm1 )436 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 437 437 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 438 438 END_3D 439 DO_2D ( 0, 0, 0, 0) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk439 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 440 440 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 441 441 END_2D 442 DO_3DS ( 0, 0, 0, 0, jpk-2, 2, -1 )442 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 443 443 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 444 444 END_3D 445 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) ! set the minimum value of tke445 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke 446 446 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 447 447 END_3D … … 456 456 ! 457 457 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 458 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )458 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 459 459 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 460 460 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 461 461 END_3D 462 462 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 463 DO_2D ( 0, 0, 0, 0)463 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 464 464 jk = nmln(ji,jj) 465 465 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & … … 467 467 END_2D 468 468 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 469 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )469 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 470 470 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 471 471 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 524 524 REAL(wp) :: zdku, zdkv, zsqen ! - - 525 525 REAL(wp) :: zemxl, zemlm, zemlp, zmaxice ! - - 526 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace526 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zmxlm, zmxld ! 3D workspace 527 527 !!-------------------------------------------------------------------- 528 528 ! … … 548 548 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 549 549 #if ! defined key_si3 && ! defined key_cice 550 DO_2D( 0, 0, 0, 0) ! No sea-ice550 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! No sea-ice 551 551 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 552 552 END_2D … … 555 555 ! 556 556 CASE( 0 ) ! No scaling under sea-ice 557 DO_2D( 0, 0, 0, 0)557 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 558 558 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 559 559 END_2D 560 560 ! 561 561 CASE( 1 ) ! scaling with constant sea-ice thickness 562 DO_2D( 0, 0, 0, 0)562 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 563 563 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 564 564 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) … … 566 566 ! 567 567 CASE( 2 ) ! scaling with mean sea-ice thickness 568 DO_2D( 0, 0, 0, 0)568 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 569 569 #if defined key_si3 570 570 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 578 578 ! 579 579 CASE( 3 ) ! scaling with max sea-ice thickness 580 DO_2D( 0, 0, 0, 0)580 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 581 581 zmaxice = MAXVAL( h_i(ji,jj,:) ) 582 582 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 587 587 #endif 588 588 ! 589 DO_2D( 0, 0, 0, 0)589 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 590 590 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 591 591 END_2D … … 596 596 ENDIF 597 597 ! 598 DO_3D( 0, 0, 0, 0, 2, jpkm1 )598 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 599 599 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 600 600 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 611 611 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 612 612 CASE ( 0 ) ! bounded by the distance to surface and bottom 613 DO_3D( 0, 0, 0, 0, 2, jpkm1 )613 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 614 614 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & 615 615 & gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) … … 622 622 ! 623 623 CASE ( 1 ) ! bounded by the vertical scale factor 624 DO_3D( 0, 0, 0, 0, 2, jpkm1 )624 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 625 625 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 626 626 zmxlm(ji,jj,jk) = zemxl … … 629 629 ! 630 630 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 631 DO_3D( 0, 0, 0, 0, 2, jpkm1 )! from the surface to the bottom :631 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : 632 632 zmxlm(ji,jj,jk) = & 633 633 & MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 634 634 END_3D 635 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface :635 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : 636 636 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 637 637 zmxlm(ji,jj,jk) = zemxl … … 640 640 ! 641 641 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 642 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! from the surface to the bottom : lup642 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : lup 643 643 zmxld(ji,jj,jk) = & 644 644 & MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 645 645 END_3D 646 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown646 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown 647 647 zmxlm(ji,jj,jk) = & 648 648 & MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 649 649 END_3D 650 DO_3D( 0, 0, 0, 0, 2, jpkm1 )650 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 651 651 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 652 652 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) … … 660 660 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 661 661 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 662 DO_3D ( 0, 0, 0, 0, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points662 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points 663 663 zsqen = SQRT( en(ji,jj,jk) ) 664 664 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 670 670 ! 671 671 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 672 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )672 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 673 673 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 674 674 END_3D … … 786 786 ! 787 787 ! !* Check of some namelist values 788 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2' )789 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 790 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 , 1 or 2' )788 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' ) 789 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' ) 790 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' ) 791 791 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 792 792 ! … … 796 796 rn_mxl0 = rmxl_min 797 797 ENDIF 798 799 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln800 801 798 ! !* depth of penetration of surface tke 802 799 IF( nn_etau /= 0 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
