Changeset 12377 for NEMO/trunk/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r11536 r12377 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) … … 215 215 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 216 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 217 DO_2D_00_00 218 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 219 END_2D 222 220 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 221 DO_2D_00_00 222 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 223 END_2D 228 224 ENDIF 229 225 ! … … 238 234 IF( ln_drg ) THEN !== friction used as top/bottom boundary condition on TKE 239 235 ! 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 236 DO_2D_00_00 237 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 238 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 239 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 240 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 & 241 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 242 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 243 END_2D 250 244 IF( ln_isfcav ) THEN ! top friction 251 DO jj = 2, jpjm1 252 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 surface 259 END DO 260 END DO 245 DO_2D_00_00 246 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 247 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 248 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 249 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 & 250 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 251 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 252 END_2D 261 253 ENDIF 262 254 ! … … 268 260 ! 269 261 ! !* total energy produce by LC : cumulative sum over jk 270 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1)262 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 271 263 DO jk = 2, jpk 272 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk)264 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 273 265 END DO 274 266 ! !* finite Langmuir Circulation depth 275 267 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 276 268 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 269 DO_3DS_11_11( jpkm1, 2, -1 ) 270 zus = zcof * taum(ji,jj) 271 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 272 END_3D 285 273 ! ! 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 274 DO_2D_11_11 275 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 276 END_2D 291 277 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 278 DO_2D_00_00 279 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 280 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 281 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 282 END_2D 283 DO_3D_00_00( 2, jpkm1 ) 284 IF ( zfr_i(ji,jj) /= 0. ) THEN 285 ! vertical velocity due to LC 286 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 287 ! ! vertical velocity due to LC 288 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i 289 ! ! TKE Langmuir circulation source term 290 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 291 ENDIF 292 ENDIF 293 END_3D 314 294 ! 315 295 ENDIF … … 323 303 ! 324 304 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 305 DO_3D_00_00( 2, jpkm1 ) 306 ! ! local Richardson number 307 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 308 ! ! inverse of Prandtl number 309 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 310 END_3D 335 311 ENDIF 336 312 ! 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 313 DO_3D_00_00( 2, jpkm1 ) 314 zcof = zfact1 * tmask(ji,jj,jk) 315 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 316 ! ! eddy coefficient (ensure numerical stability) 317 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 318 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) ) 319 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 320 & / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) ) 321 ! 322 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 323 zd_lw(ji,jj,jk) = zzd_lw 324 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 325 ! 326 ! ! right hand side in en 327 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 328 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 329 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 330 & ) * wmask(ji,jj,jk) 331 END_3D 360 332 ! !* 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 333 DO_3D_00_00( 3, jpkm1 ) 334 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 335 END_3D 336 DO_2D_00_00 337 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 338 END_2D 339 DO_3D_00_00( 3, jpkm1 ) 340 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) 341 END_3D 342 DO_2D_00_00 343 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 344 END_2D 345 DO_3DS_00_00( jpk-2, 2, -1 ) 346 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 347 END_3D 348 DO_3D_00_00( 2, jpkm1 ) 349 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 350 END_3D 399 351 ! 400 352 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 402 354 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 403 355 !!gm BUG : in the exp remove the depth of ssh !!! 404 !!gm i.e. use gde3w in argument ( pdepw)356 !!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 405 357 406 358 407 359 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 360 DO_3D_00_00( 2, jpkm1 ) 361 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 362 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 363 END_3D 416 364 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 365 DO_2D_00_00 366 jk = nmln(ji,jj) 367 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 368 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 369 END_2D 424 370 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 371 DO_3D_00_00( 2, jpkm1 ) 372 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 373 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 374 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 375 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 376 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 377 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 378 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 379 END_3D 438 380 ENDIF 439 381 ! … … 441 383 442 384 443 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt )385 SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 444 386 !!---------------------------------------------------------------------- 445 387 !! *** ROUTINE tke_avn *** … … 477 419 USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics 478 420 !! 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) 421 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 481 422 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 482 423 ! … … 500 441 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 501 442 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 443 DO_2D_00_00 444 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 445 END_2D 507 446 ELSE 508 447 zmxlm(:,:,1) = rn_mxl0 509 448 ENDIF 510 449 ! 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 450 DO_3D_00_00( 2, jpkm1 ) 451 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 452 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 453 END_3D 519 454 ! 520 455 ! !* Physical limits for the mixing length … … 526 461 ! 527 462 !!gm Not sure of that coding for ISF.... 528 ! where wmask = 0 set zmxlm == p_e3w463 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 529 464 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 465 DO_3D_00_00( 2, jpkm1 ) 466 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & 467 & gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 468 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 469 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 470 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 471 END_3D 541 472 ! 542 473 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 474 DO_3D_00_00( 2, jpkm1 ) 475 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 476 zmxlm(ji,jj,jk) = zemxl 477 zmxld(ji,jj,jk) = zemxl 478 END_3D 552 479 ! 553 480 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 481 DO_3D_00_00( 2, jpkm1 ) 482 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 483 END_3D 484 DO_3DS_00_00( jpkm1, 2, -1 ) 485 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 486 zmxlm(ji,jj,jk) = zemxl 487 zmxld(ji,jj,jk) = zemxl 488 END_3D 570 489 ! 571 490 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 491 DO_3D_00_00( 2, jpkm1 ) 492 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 493 END_3D 494 DO_3DS_00_00( jpkm1, 2, -1 ) 495 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 496 END_3D 497 DO_3D_00_00( 2, jpkm1 ) 498 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 499 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 500 zmxlm(ji,jj,jk) = zemlm 501 zmxld(ji,jj,jk) = zemlp 502 END_3D 596 503 ! 597 504 END SELECT … … 600 507 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 601 508 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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 509 DO_3D_00_00( 1, jpkm1 ) 510 zsqen = SQRT( en(ji,jj,jk) ) 511 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 512 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 513 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 514 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 515 END_3D 613 516 ! 614 517 ! 615 518 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 519 DO_3D_00_00( 2, jpkm1 ) 520 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) 521 END_3D 522 ENDIF 523 ! 524 IF(sn_cfctl%l_prtctl) THEN 626 525 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 627 526 CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke - m: ', kdim=jpk ) … … 631 530 632 531 633 SUBROUTINE zdf_tke_init 532 SUBROUTINE zdf_tke_init( Kmm ) 634 533 !!---------------------------------------------------------------------- 635 534 !! *** ROUTINE zdf_tke_init *** … … 647 546 USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag 648 547 !! 649 INTEGER :: ji, jj, jk ! dummy loop indices 650 INTEGER :: ios 548 INTEGER, INTENT(in) :: Kmm ! time level index 549 INTEGER :: ji, jj, jk ! dummy loop indices 550 INTEGER :: ios 651 551 !! 652 552 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & … … 656 556 !!---------------------------------------------------------------------- 657 557 ! 658 REWIND( numnam_ref ) ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy659 558 READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 660 559 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' ) 661 560 662 REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy663 561 READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 ) 664 562 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' ) … … 725 623 ENDIF 726 624 727 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln625 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln 728 626 729 627 ! !* depth of penetration of surface tke
Note: See TracChangeset
for help on using the changeset viewer.