Changeset 14834 for NEMO/trunk/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2021-05-11T11:24:44+02:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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.
