- Timestamp:
- 2021-04-27T17:33:44+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdftke.F90
r14601 r14757 241 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 242 242 ! 243 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )244 243 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 245 244 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) … … 259 258 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 260 259 ! 261 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! bottom friction262 260 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction 263 261 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) … … 269 267 END_2D 270 268 IF( ln_isfcav ) THEN 271 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! top friction272 269 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 273 270 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) … … 297 294 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 298 295 !!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect ! 299 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )300 296 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 301 297 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) … … 305 301 ! Projection of Stokes drift in the wind stress direction 306 302 ! 307 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )308 303 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 309 304 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) … … 312 307 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 313 308 END_2D 314 ! [comm_cleanup] 315 IF (nn_hls.eq.1) CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. ) 309 IF (nn_hls==1) CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. ) 316 310 ! 317 311 ELSE ! Surface Stokes drift deduced from surface stress … … 321 315 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 322 316 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 323 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )324 317 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 325 318 zWlc2(ji,jj) = zcof * taum(ji,jj) … … 338 331 ! !- compare LHS to RHS of Eq.47 339 332 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 340 ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )341 333 DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 ) 342 334 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 343 335 END_3D 344 336 ! ! finite LC depth 345 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )346 337 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 347 338 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) … … 349 340 ! 350 341 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 351 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )352 342 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 353 343 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 354 344 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 355 345 END_2D 356 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en357 346 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 358 347 IF ( zus3(ji,jj) /= 0._wp ) THEN … … 376 365 ! 377 366 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 378 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )379 367 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 380 368 ! ! local Richardson number … … 389 377 ENDIF 390 378 ! 391 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en392 379 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en 393 380 zcof = zfact1 * tmask(ji,jj,jk) … … 419 406 420 407 CASE ( 0 ) ! Dirichlet BC 421 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0)422 408 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 423 409 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp … … 427 413 428 414 CASE ( 1 ) ! Neumann BC 429 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )430 415 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 431 416 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp … … 442 427 ! 443 428 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 444 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1445 429 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 446 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) … … 450 434 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 451 435 ! END_2D 452 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )453 436 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 454 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) 455 438 END_3D 456 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk457 439 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 458 440 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 459 441 END_2D 460 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )461 442 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 462 443 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 463 444 END_3D 464 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! set the minimum value of tke465 445 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke 466 446 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) … … 476 456 ! 477 457 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 478 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )479 458 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 480 459 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & … … 482 461 END_3D 483 462 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 484 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )485 463 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 486 464 jk = nmln(ji,jj) … … 489 467 END_2D 490 468 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 491 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )492 469 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 493 470 ztx2 = utau(ji-1,jj ) + utau(ji,jj) … … 571 548 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 572 549 #if ! defined key_si3 && ! defined key_cice 573 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! No sea-ice574 550 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! No sea-ice 575 551 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) … … 579 555 ! 580 556 CASE( 0 ) ! No scaling under sea-ice 581 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )582 557 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 583 558 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) … … 585 560 ! 586 561 CASE( 1 ) ! scaling with constant sea-ice thickness 587 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )588 562 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 589 563 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 592 566 ! 593 567 CASE( 2 ) ! scaling with mean sea-ice thickness 594 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )595 568 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 596 569 #if defined key_si3 … … 605 578 ! 606 579 CASE( 3 ) ! scaling with max sea-ice thickness 607 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )608 580 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 609 581 zmaxice = MAXVAL( h_i(ji,jj,:) ) … … 615 587 #endif 616 588 ! 617 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 )618 589 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 619 590 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) … … 625 596 ENDIF 626 597 ! 627 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )628 598 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 629 599 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) … … 641 611 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 642 612 CASE ( 0 ) ! bounded by the distance to surface and bottom 643 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )644 613 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 645 614 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & … … 653 622 ! 654 623 CASE ( 1 ) ! bounded by the vertical scale factor 655 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )656 624 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 657 625 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) … … 661 629 ! 662 630 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 663 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! from the surface to the bottom :664 631 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : 665 632 zmxlm(ji,jj,jk) = & 666 633 & MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 667 634 END_3D 668 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface :669 635 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : 670 636 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) … … 674 640 ! 675 641 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 676 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! from the surface to the bottom : lup677 642 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : lup 678 643 zmxld(ji,jj,jk) = & 679 644 & MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 680 645 END_3D 681 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown682 646 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown 683 647 zmxlm(ji,jj,jk) = & 684 648 & MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 685 649 END_3D 686 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )687 650 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 688 651 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) … … 697 660 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 698 661 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 699 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points700 662 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points 701 663 zsqen = SQRT( en(ji,jj,jk) ) … … 708 670 ! 709 671 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 710 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 )711 672 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 712 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)
Note: See TracChangeset
for help on using the changeset viewer.
