- Timestamp:
- 2021-03-25T12:51:33+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
r14156 r14636 179 179 180 180 ! Compute surface, top and bottom friction at T-points 181 DO_2D( 0, 0, 0, 0 ) !== surface ocean friction 181 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) !== surface ocean friction 182 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== surface ocean friction 182 183 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) ! surface friction 183 184 END_2D … … 186 187 ! 187 188 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 DO_2D( 0, 0, 0, 0 ) ! bottom friction (explicit before friction) 189 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! bottom friction (explicit before friction) 190 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction (explicit before friction) 189 191 zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 190 192 zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) … … 193 195 END_2D 194 196 IF( ln_isfcav ) THEN 195 DO_2D( 0, 0, 0, 0 ) ! top friction 197 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! top friction 198 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 196 199 zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 197 200 zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) … … 220 223 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 221 224 ! 222 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! 225 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! 226 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !== Compute dissipation rate ==! 223 227 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 224 228 END_3D … … 229 233 230 234 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 231 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 235 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 236 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 232 237 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 233 238 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) … … 250 255 ! Warning : after this step, en : right hand side of the matrix 251 256 252 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 257 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 258 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 253 259 ! 254 260 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 327 333 ! at k=2, set de/dz=Fw 328 334 !cbr 329 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 335 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 336 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! zdiag zd_lw not defined/used on the halo 330 337 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 331 338 zd_lw(ji,jj,2) = 0._wp … … 348 355 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 349 356 ! ! Balance between the production and the dissipation terms 350 DO_2D( 0, 0, 0, 0 ) 357 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 358 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 351 359 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 352 360 !! With thick deep ocean level thickness, this may be quite large, no ??? … … 366 374 ! 367 375 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 368 DO_2D( 0, 0, 0, 0 ) 376 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 377 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 369 378 itop = mikt(ji,jj) ! k top w-point 370 379 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 384 393 CASE ( 1 ) ! Neumman boundary condition 385 394 ! 386 DO_2D( 0, 0, 0, 0 ) 395 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 396 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 387 397 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 388 398 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 399 409 END_2D 400 410 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 401 DO_2D( 0, 0, 0, 0 ) 411 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 412 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 402 413 itop = mikt(ji,jj) ! k top w-point 403 414 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 420 431 ! ---------------------------------------------------------- 421 432 ! 422 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 433 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 434 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 423 435 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 424 436 END_3D 425 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 437 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 438 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 426 439 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) 427 440 END_3D 428 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 441 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 442 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 429 443 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 430 444 END_3D … … 441 455 ! 442 456 CASE( 0 ) ! k-kl (Mellor-Yamada) 443 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 457 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 458 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 444 459 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 445 460 END_3D 446 461 ! 447 462 CASE( 1 ) ! k-eps 448 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 463 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 464 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 449 465 psi(ji,jj,jk) = eps(ji,jj,jk) 450 466 END_3D 451 467 ! 452 468 CASE( 2 ) ! k-w 453 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 469 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 470 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 454 471 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 455 472 END_3D 456 473 ! 457 474 CASE( 3 ) ! generic 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 475 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 476 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 459 477 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 460 478 END_3D … … 469 487 ! Warning : after this step, en : right hand side of the matrix 470 488 471 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 489 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 490 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 472 491 ! 473 492 ! psi / k … … 541 560 ! 542 561 ! Neumann condition at k=2 543 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 562 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 563 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! zdiag zd_lw not defined/used on the halo 544 564 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 545 565 zd_lw(ji,jj,2) = 0._wp … … 569 589 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 570 590 ! ! Balance between the production and the dissipation terms 571 DO_2D( 0, 0, 0, 0 ) 591 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 592 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 572 593 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 573 594 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 588 609 CASE ( 1 ) ! Neumman boundary condition 589 610 ! 590 DO_2D( 0, 0, 0, 0 ) 611 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 612 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 591 613 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 592 614 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 616 638 ! ---------------- 617 639 ! 618 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 640 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 641 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 619 642 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 620 643 END_3D 621 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 644 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 645 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 622 646 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 623 647 END_3D 624 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 648 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 649 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 625 650 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 626 651 END_3D … … 632 657 ! 633 658 CASE( 0 ) ! k-kl (Mellor-Yamada) 634 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 659 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 660 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 635 661 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 636 662 END_3D 637 663 ! 638 664 CASE( 1 ) ! k-eps 639 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 665 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 666 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 640 667 eps(ji,jj,jk) = psi(ji,jj,jk) 641 668 END_3D 642 669 ! 643 670 CASE( 2 ) ! k-w 644 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 671 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 672 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 645 673 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 646 674 END_3D … … 650 678 zex1 = ( 1.5_wp + rmm/rnn ) 651 679 zex2 = -1._wp / rnn 652 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 680 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 681 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 653 682 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 654 683 END_3D … … 658 687 ! Limit dissipation rate under stable stratification 659 688 ! -------------------------------------------------- 660 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 689 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 690 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 661 691 ! limitation 662 692 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) … … 674 704 ! 675 705 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 676 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 706 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 707 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 677 708 ! zcof = l²/q² 678 709 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) … … 691 722 ! 692 723 CASE ( 2, 3 ) ! Canuto stability functions 693 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 724 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 725 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 694 726 ! zcof = l²/q² 695 727 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) … … 723 755 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 724 756 zstm(:,:,jpk) = 0. 725 DO_2D( 0, 0, 0, 0 ) ! update bottom with good values 757 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! update bottom with good values 758 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! update bottom with good values 726 759 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 727 760 END_2D … … 738 771 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 739 772 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 740 DO_3D( 0, 0, 0, 0, 1, jpk ) 773 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk ) 774 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 741 775 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 742 776 zavt = zsqen * zstt(ji,jj,jk)
Note: See TracChangeset
for help on using the changeset viewer.
