- Timestamp:
- 2021-12-03T20:32:50+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdftke.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ZDF/zdftke.F90
r14072 r15574 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) :: imlc 219 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2 220 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 218 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 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztmp ! for diags 221 222 !!-------------------------------------------------------------------- 222 223 ! … … 232 233 SELECT CASE ( nn_eice ) 233 234 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 )235 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 236 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 237 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 237 238 END SELECT 238 239 ! … … 241 242 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 242 243 ! 243 DO_2D ( 0, 0, 0, 0)244 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 244 245 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) 245 246 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) … … 258 259 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 259 260 ! 260 DO_2D ( 0, 0, 0, 0) ! bottom friction261 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction 261 262 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 262 263 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) … … 267 268 END_2D 268 269 IF( ln_isfcav ) THEN 269 DO_2D ( 0, 0, 0, 0) ! top friction270 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 270 271 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 271 272 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) … … 294 295 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 295 296 !!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)297 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 297 298 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) 298 299 zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) … … 301 302 ! Projection of Stokes drift in the wind stress direction 302 303 ! 303 DO_2D( 0, 0, 0, 0)304 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 304 305 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 305 306 ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) … … 307 308 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 308 309 END_2D 309 CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. )310 !311 310 ELSE ! Surface Stokes drift deduced from surface stress 312 311 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) … … 315 314 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 316 315 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 )316 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 318 317 zWlc2(ji,jj) = zcof * taum(ji,jj) 319 318 END_2D … … 323 322 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 324 323 ! !- 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 324 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 325 zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm) 326 END_2D 327 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk ) 328 zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + & 329 & MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 330 END_3D 330 331 ! 331 332 ! !- 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 )333 imlc(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point (=2 over land) 334 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) 334 335 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 335 336 END_3D 336 337 ! ! finite LC depth 337 DO_2D( 1, 1, 1,1 )338 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 338 339 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 339 340 END_2D 340 341 ! 341 342 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 342 DO_2D( 0, 0, 0, 0)343 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 343 344 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 344 345 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 345 346 END_2D 346 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en347 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 348 IF ( zus3(ji,jj) /= 0._wp ) THEN 348 349 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN … … 365 366 ! 366 367 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 367 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )368 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 368 369 ! ! local Richardson number 369 370 IF (rn2b(ji,jj,jk) <= 0.0_wp) then … … 377 378 ENDIF 378 379 ! 379 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en380 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 381 zcof = zfact1 * tmask(ji,jj,jk) 381 382 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical … … 406 407 407 408 CASE ( 0 ) ! Dirichlet BC 408 DO_2D ( 0, 0, 0, 0) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0)409 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 410 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 410 411 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) … … 413 414 414 415 CASE ( 1 ) ! Neumann BC 415 DO_2D ( 0, 0, 0, 0)416 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 416 417 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 417 418 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) … … 427 428 ! 428 429 ! !* 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-1430 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 430 431 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 432 END_3D … … 434 435 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 435 436 ! END_2D 436 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )437 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 437 438 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 439 END_3D 439 DO_2D ( 0, 0, 0, 0) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk440 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 440 441 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 441 442 END_2D 442 DO_3DS ( 0, 0, 0, 0, jpk-2, 2, -1 )443 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 443 444 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 444 445 END_3D 445 DO_3D ( 0, 0, 0, 0, 2, jpkm1 ) ! set the minimum value of tke446 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke 446 447 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 447 448 END_3D 449 ! 450 ! Kolmogorov energy of dissipation (W/kg) 451 ! ediss = Ce*sqrt(en)/L*en 452 ! dissl = sqrt(en)/L 453 IF( iom_use('ediss_k') ) THEN 454 ALLOCATE( ztmp(A2D(nn_hls),jpk) ) 455 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 456 ztmp(ji,jj,jk) = zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) * wmask(ji,jj,jk) 457 END_3D 458 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 459 ztmp(ji,jj,jpk) = 0._wp 460 END_2D 461 CALL iom_put( 'ediss_k', ztmp ) 462 DEALLOCATE( ztmp ) 463 ENDIF 448 464 ! 449 465 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 456 472 ! 457 473 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 458 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )474 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 459 475 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 460 476 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 461 477 END_3D 462 478 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)479 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 464 480 jk = nmln(ji,jj) 465 481 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & … … 467 483 END_2D 468 484 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 469 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )485 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 470 486 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 471 487 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 524 540 REAL(wp) :: zdku, zdkv, zsqen ! - - 525 541 REAL(wp) :: zemxl, zemlm, zemlp, zmaxice ! - - 526 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace542 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zmxlm, zmxld ! 3D workspace 527 543 !!-------------------------------------------------------------------- 528 544 ! … … 548 564 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 549 565 #if ! defined key_si3 && ! defined key_cice 550 DO_2D( 0, 0, 0, 0) ! No sea-ice566 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! No sea-ice 551 567 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 552 568 END_2D … … 555 571 ! 556 572 CASE( 0 ) ! No scaling under sea-ice 557 DO_2D( 0, 0, 0, 0)573 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 558 574 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 559 575 END_2D 560 576 ! 561 577 CASE( 1 ) ! scaling with constant sea-ice thickness 562 DO_2D( 0, 0, 0, 0)578 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 563 579 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 564 580 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) … … 566 582 ! 567 583 CASE( 2 ) ! scaling with mean sea-ice thickness 568 DO_2D( 0, 0, 0, 0)584 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 569 585 #if defined key_si3 570 586 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 578 594 ! 579 595 CASE( 3 ) ! scaling with max sea-ice thickness 580 DO_2D( 0, 0, 0, 0)596 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 581 597 zmaxice = MAXVAL( h_i(ji,jj,:) ) 582 598 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 587 603 #endif 588 604 ! 589 DO_2D( 0, 0, 0, 0)605 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 590 606 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 591 607 END_2D … … 596 612 ENDIF 597 613 ! 598 DO_3D( 0, 0, 0, 0, 2, jpkm1 )614 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 599 615 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 600 616 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 611 627 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 612 628 CASE ( 0 ) ! bounded by the distance to surface and bottom 613 DO_3D( 0, 0, 0, 0, 2, jpkm1 )629 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 614 630 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & 615 631 & gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) … … 622 638 ! 623 639 CASE ( 1 ) ! bounded by the vertical scale factor 624 DO_3D( 0, 0, 0, 0, 2, jpkm1 )640 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 625 641 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 626 642 zmxlm(ji,jj,jk) = zemxl … … 629 645 ! 630 646 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 631 DO_3D( 0, 0, 0, 0, 2, jpkm1 )! from the surface to the bottom :647 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : 632 648 zmxlm(ji,jj,jk) = & 633 649 & MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 634 650 END_3D 635 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface :651 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : 636 652 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 637 653 zmxlm(ji,jj,jk) = zemxl … … 640 656 ! 641 657 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 : lup658 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : lup 643 659 zmxld(ji,jj,jk) = & 644 660 & MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 645 661 END_3D 646 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown662 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 663 zmxlm(ji,jj,jk) = & 648 664 & MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 649 665 END_3D 650 DO_3D( 0, 0, 0, 0, 2, jpkm1 )666 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 651 667 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 652 668 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) … … 660 676 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 661 677 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 662 DO_3D ( 0, 0, 0, 0, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points678 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 679 zsqen = SQRT( en(ji,jj,jk) ) 664 680 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 670 686 ! 671 687 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 672 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )688 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 673 689 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 690 END_3D … … 786 802 ! 787 803 ! !* 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' )804 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' ) 805 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' ) 806 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' ) 791 807 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 792 808 ! … … 796 812 rn_mxl0 = rmxl_min 797 813 ENDIF 798 799 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln800 801 814 ! !* depth of penetration of surface tke 802 815 IF( nn_etau /= 0 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
