- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdftke.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdftke.F90
r12178 r12928 89 89 90 90 !! * Substitutions 91 # include " vectopt_loop_substitute.h90"91 # include "do_loop_substitute.h90" 92 92 !!---------------------------------------------------------------------- 93 93 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 109 109 110 110 111 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt )111 SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 112 112 !!---------------------------------------------------------------------- 113 113 !! *** ROUTINE zdf_tke *** … … 155 155 !!---------------------------------------------------------------------- 156 156 INTEGER , INTENT(in ) :: kt ! ocean time step 157 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 157 158 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 158 159 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 159 160 !!---------------------------------------------------------------------- 160 161 ! 161 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en)162 ! 163 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl162 CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) ! now tke (en) 163 ! 164 CALL tke_avn( Kbb, Kmm, p_avm, p_avt ) ! now avt, avm, dissl 164 165 ! 165 166 END SUBROUTINE zdf_tke 166 167 167 168 168 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt )169 SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) 169 170 !!---------------------------------------------------------------------- 170 171 !! *** ROUTINE tke_tke *** … … 186 187 USE zdf_oce , ONLY : en ! ocean vertical physics 187 188 !! 188 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 189 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 189 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 190 190 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 191 191 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) … … 206 206 !!-------------------------------------------------------------------- 207 207 ! 208 zbbrau = rn_ebb / r au0 ! Local constant initialisation209 zfact1 = -.5_wp * r dt210 zfact2 = 1.5_wp * r dt * rn_ediss208 zbbrau = rn_ebb / rho0 ! Local constant initialisation 209 zfact1 = -.5_wp * rn_Dt 210 zfact2 = 1.5_wp * rn_Dt * rn_ediss 211 211 zfact3 = 0.5_wp * rn_ediss 212 212 ! … … 214 214 ! ! Surface/top/bottom boundary condition on tke 215 215 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 216 217 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 218 DO ji = fs_2, fs_jpim1 ! vector opt. 219 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 220 END DO 221 END DO 222 IF ( ln_isfcav ) THEN 223 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 224 DO ji = fs_2, fs_jpim1 ! vector opt. 225 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 226 END DO 227 END DO 228 ENDIF 216 ! 217 DO_2D_00_00 218 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 219 END_2D 229 220 ! 230 221 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 232 223 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 233 224 ! 234 ! en(bot) = (ebb0/r au0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)225 ! en(bot) = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 235 226 ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 236 227 ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 … … 238 229 IF( ln_drg ) THEN !== friction used as top/bottom boundary condition on TKE 239 230 ! 240 DO jj = 2, jpjm1 ! bottom friction 241 DO ji = fs_2, fs_jpim1 ! vector opt. 242 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 243 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 244 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 245 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 246 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 247 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 248 END DO 249 END DO 231 DO_2D_00_00 232 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 233 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 234 ! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 235 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 236 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 237 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 238 END_2D 250 239 IF( ln_isfcav ) THEN ! top friction 251 DO jj = 2, jpjm1252 DO ji = fs_2, fs_jpim1 ! vector opt.253 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )254 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)))255 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)256 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 &257 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 )258 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface259 END DO260 END DO240 DO_2D_00_00 241 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 242 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 243 ! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 244 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 245 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 246 ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present 247 en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) & 248 & + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 249 END_2D 261 250 ENDIF 262 251 ! … … 268 257 ! 269 258 ! !* total energy produce by LC : cumulative sum over jk 270 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1)259 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 271 260 DO jk = 2, jpk 272 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk)261 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 273 262 END DO 274 263 ! !* finite Langmuir Circulation depth 275 264 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 276 265 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 277 DO jk = jpkm1, 2, -1 278 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 279 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 280 zus = zcof * taum(ji,jj) 281 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 282 END DO 283 END DO 284 END DO 266 DO_3DS_11_11( jpkm1, 2, -1 ) 267 zus = zcof * taum(ji,jj) 268 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 269 END_3D 285 270 ! ! finite LC depth 286 DO jj = 1, jpj 287 DO ji = 1, jpi 288 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 289 END DO 290 END DO 271 DO_2D_11_11 272 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 273 END_2D 291 274 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 292 DO jj = 2, jpjm1 293 DO ji = fs_2, fs_jpim1 ! vector opt. 294 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 295 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 296 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 297 END DO 298 END DO 299 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 300 DO jj = 2, jpjm1 301 DO ji = fs_2, fs_jpim1 ! vector opt. 302 IF ( zfr_i(ji,jj) /= 0. ) THEN 303 ! vertical velocity due to LC 304 IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 305 ! ! vertical velocity due to LC 306 zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i 307 ! ! TKE Langmuir circulation source term 308 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 309 ENDIF 310 ENDIF 311 END DO 312 END DO 313 END DO 275 DO_2D_00_00 276 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 277 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 278 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 279 END_2D 280 DO_3D_00_00( 2, jpkm1 ) 281 IF ( zfr_i(ji,jj) /= 0. ) THEN 282 ! vertical velocity due to LC 283 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 284 ! ! vertical velocity due to LC 285 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i 286 ! ! TKE Langmuir circulation source term 287 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 288 ENDIF 289 ENDIF 290 END_3D 314 291 ! 315 292 ENDIF … … 323 300 ! 324 301 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 325 DO jk = 2, jpkm1 326 DO jj = 2, jpjm1 327 DO ji = 2, jpim1 328 ! ! local Richardson number 329 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 330 ! ! inverse of Prandtl number 331 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 332 END DO 333 END DO 334 END DO 302 DO_3D_00_00( 2, jpkm1 ) 303 ! ! local Richardson number 304 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 305 ! ! inverse of Prandtl number 306 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 307 END_3D 335 308 ENDIF 336 309 ! 337 DO jk = 2, jpkm1 !* Matrix and right hand side in en 338 DO jj = 2, jpjm1 339 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zcof = zfact1 * tmask(ji,jj,jk) 341 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 342 ! ! eddy coefficient (ensure numerical stability) 343 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 344 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 345 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 346 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 347 ! 348 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 349 zd_lw(ji,jj,jk) = zzd_lw 350 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 351 ! 352 ! ! right hand side in en 353 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 354 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 355 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 356 & ) * wmask(ji,jj,jk) 357 END DO 358 END DO 359 END DO 310 DO_3D_00_00( 2, jpkm1 ) 311 zcof = zfact1 * tmask(ji,jj,jk) 312 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 313 ! ! eddy coefficient (ensure numerical stability) 314 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 315 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) ) 316 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 317 & / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) ) 318 ! 319 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 320 zd_lw(ji,jj,jk) = zzd_lw 321 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 322 ! 323 ! ! right hand side in en 324 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * ( p_sh2(ji,jj,jk) & ! shear 325 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 326 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 327 & ) * wmask(ji,jj,jk) 328 END_3D 360 329 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 361 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 362 DO jj = 2, jpjm1 363 DO ji = fs_2, fs_jpim1 ! vector opt. 364 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 365 END DO 366 END DO 367 END DO 368 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 370 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 371 END DO 372 END DO 373 DO jk = 3, jpkm1 374 DO jj = 2, jpjm1 375 DO ji = fs_2, fs_jpim1 ! vector opt. 376 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) 377 END DO 378 END DO 379 END DO 380 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 381 DO ji = fs_2, fs_jpim1 ! vector opt. 382 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 383 END DO 384 END DO 385 DO jk = jpk-2, 2, -1 386 DO jj = 2, jpjm1 387 DO ji = fs_2, fs_jpim1 ! vector opt. 388 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 389 END DO 390 END DO 391 END DO 392 DO jk = 2, jpkm1 ! set the minimum value of tke 393 DO jj = 2, jpjm1 394 DO ji = fs_2, fs_jpim1 ! vector opt. 395 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 396 END DO 397 END DO 398 END DO 330 DO_3D_00_00( 3, jpkm1 ) 331 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 332 END_3D 333 DO_2D_00_00 334 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 335 END_2D 336 DO_3D_00_00( 3, jpkm1 ) 337 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) 338 END_3D 339 DO_2D_00_00 340 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 341 END_2D 342 DO_3DS_00_00( jpk-2, 2, -1 ) 343 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 344 END_3D 345 DO_3D_00_00( 2, jpkm1 ) 346 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 347 END_3D 399 348 ! 400 349 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 402 351 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 403 352 !!gm BUG : in the exp remove the depth of ssh !!! 404 !!gm i.e. use gde3w in argument ( pdepw)353 !!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 405 354 406 355 407 356 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 408 DO jk = 2, jpkm1 ! rn_eice =0 ON below sea-ice, =4 OFF when ice fraction > 0.25 409 DO jj = 2, jpjm1 410 DO ji = fs_2, fs_jpim1 ! vector opt. 411 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 412 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 413 END DO 414 END DO 415 END DO 357 DO_3D_00_00( 2, jpkm1 ) 358 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 359 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 360 END_3D 416 361 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 417 DO jj = 2, jpjm1 418 DO ji = fs_2, fs_jpim1 ! vector opt. 419 jk = nmln(ji,jj) 420 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 421 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 422 END DO 423 END DO 362 DO_2D_00_00 363 jk = nmln(ji,jj) 364 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 365 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 366 END_2D 424 367 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 425 DO jk = 2, jpkm1 426 DO jj = 2, jpjm1 427 DO ji = fs_2, fs_jpim1 ! vector opt. 428 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 429 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 430 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 431 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 432 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 433 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 434 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 435 END DO 436 END DO 437 END DO 368 DO_3D_00_00( 2, jpkm1 ) 369 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 370 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 371 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 372 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 373 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 374 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 375 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 376 END_3D 438 377 ENDIF 439 378 ! … … 441 380 442 381 443 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt )382 SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 444 383 !!---------------------------------------------------------------------- 445 384 !! *** ROUTINE tke_avn *** … … 477 416 USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics 478 417 !! 479 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 480 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 418 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 481 419 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 482 420 ! … … 498 436 zmxld(:,:,:) = rmxl_min 499 437 ! 500 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 501 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 502 DO jj = 2, jpjm1 503 DO ji = fs_2, fs_jpim1 504 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 505 END DO 506 END DO 438 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 439 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 440 DO_2D_00_00 441 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 442 END_2D 507 443 ELSE 508 444 zmxlm(:,:,1) = rn_mxl0 509 445 ENDIF 510 446 ! 511 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 512 DO jj = 2, jpjm1 513 DO ji = fs_2, fs_jpim1 ! vector opt. 514 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 515 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 516 END DO 517 END DO 518 END DO 447 DO_3D_00_00( 2, jpkm1 ) 448 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 449 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 450 END_3D 519 451 ! 520 452 ! !* Physical limits for the mixing length … … 526 458 ! 527 459 !!gm Not sure of that coding for ISF.... 528 ! where wmask = 0 set zmxlm == p_e3w460 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 529 461 CASE ( 0 ) ! bounded by the distance to surface and bottom 530 DO jk = 2, jpkm1 531 DO jj = 2, jpjm1 532 DO ji = fs_2, fs_jpim1 ! vector opt. 533 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 534 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 535 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 536 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 537 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 538 END DO 539 END DO 540 END DO 462 DO_3D_00_00( 2, jpkm1 ) 463 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & 464 & gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 465 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 466 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 467 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 468 END_3D 541 469 ! 542 470 CASE ( 1 ) ! bounded by the vertical scale factor 543 DO jk = 2, jpkm1 544 DO jj = 2, jpjm1 545 DO ji = fs_2, fs_jpim1 ! vector opt. 546 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 547 zmxlm(ji,jj,jk) = zemxl 548 zmxld(ji,jj,jk) = zemxl 549 END DO 550 END DO 551 END DO 471 DO_3D_00_00( 2, jpkm1 ) 472 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 473 zmxlm(ji,jj,jk) = zemxl 474 zmxld(ji,jj,jk) = zemxl 475 END_3D 552 476 ! 553 477 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 554 DO jk = 2, jpkm1 ! from the surface to the bottom : 555 DO jj = 2, jpjm1 556 DO ji = fs_2, fs_jpim1 ! vector opt. 557 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 558 END DO 559 END DO 560 END DO 561 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 562 DO jj = 2, jpjm1 563 DO ji = fs_2, fs_jpim1 ! vector opt. 564 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 565 zmxlm(ji,jj,jk) = zemxl 566 zmxld(ji,jj,jk) = zemxl 567 END DO 568 END DO 569 END DO 478 DO_3D_00_00( 2, jpkm1 ) 479 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 480 END_3D 481 DO_3DS_00_00( jpkm1, 2, -1 ) 482 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 483 zmxlm(ji,jj,jk) = zemxl 484 zmxld(ji,jj,jk) = zemxl 485 END_3D 570 486 ! 571 487 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 572 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 573 DO jj = 2, jpjm1 574 DO ji = fs_2, fs_jpim1 ! vector opt. 575 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 576 END DO 577 END DO 578 END DO 579 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 580 DO jj = 2, jpjm1 581 DO ji = fs_2, fs_jpim1 ! vector opt. 582 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 583 END DO 584 END DO 585 END DO 586 DO jk = 2, jpkm1 587 DO jj = 2, jpjm1 588 DO ji = fs_2, fs_jpim1 ! vector opt. 589 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 590 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 591 zmxlm(ji,jj,jk) = zemlm 592 zmxld(ji,jj,jk) = zemlp 593 END DO 594 END DO 595 END DO 488 DO_3D_00_00( 2, jpkm1 ) 489 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 490 END_3D 491 DO_3DS_00_00( jpkm1, 2, -1 ) 492 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 493 END_3D 494 DO_3D_00_00( 2, jpkm1 ) 495 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 496 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 497 zmxlm(ji,jj,jk) = zemlm 498 zmxld(ji,jj,jk) = zemlp 499 END_3D 596 500 ! 597 501 END SELECT … … 600 504 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 601 505 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 602 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 603 DO jj = 2, jpjm1 604 DO ji = fs_2, fs_jpim1 ! vector opt. 605 zsqen = SQRT( en(ji,jj,jk) ) 606 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 607 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 608 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 609 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 610 END DO 611 END DO 612 END DO 506 DO_3D_00_00( 1, jpkm1 ) 507 zsqen = SQRT( en(ji,jj,jk) ) 508 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 509 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 510 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 511 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 512 END_3D 613 513 ! 614 514 ! 615 515 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 616 DO jk = 2, jpkm1 617 DO jj = 2, jpjm1 618 DO ji = fs_2, fs_jpim1 ! vector opt. 619 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 620 END DO 621 END DO 622 END DO 623 ENDIF 624 ! 625 IF(ln_ctl) THEN 516 DO_3D_00_00( 2, jpkm1 ) 517 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) 518 END_3D 519 ENDIF 520 ! 521 IF(sn_cfctl%l_prtctl) THEN 626 522 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 627 523 CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke - m: ', kdim=jpk ) … … 631 527 632 528 633 SUBROUTINE zdf_tke_init 529 SUBROUTINE zdf_tke_init( Kmm ) 634 530 !!---------------------------------------------------------------------- 635 531 !! *** ROUTINE zdf_tke_init *** … … 647 543 USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag 648 544 !! 649 INTEGER :: ji, jj, jk ! dummy loop indices 650 INTEGER :: ios 545 INTEGER, INTENT(in) :: Kmm ! time level index 546 INTEGER :: ji, jj, jk ! dummy loop indices 547 INTEGER :: ios 651 548 !! 652 549 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & … … 656 553 !!---------------------------------------------------------------------- 657 554 ! 658 REWIND( numnam_ref ) ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy659 555 READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 660 556 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' ) 661 557 662 REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy663 558 READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 ) 664 559 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' ) … … 725 620 ENDIF 726 621 727 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln622 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln 728 623 729 624 ! !* depth of penetration of surface tke
Note: See TracChangeset
for help on using the changeset viewer.
