Changeset 13295 for NEMO/trunk/src/OCE/ZDF/zdfgls.F90
- Timestamp:
- 2020-07-10T20:24:21+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdfgls.F90
r13286 r13295 168 168 169 169 ! Compute surface, top and bottom friction at T-points 170 DO_2D _00_00170 DO_2D( 0, 0, 0, 0 ) 171 171 ! 172 172 ! surface friction … … 181 181 END_2D 182 182 IF( ln_isfcav ) THEN !top friction 183 DO_2D _00_00183 DO_2D( 0, 0, 0, 0 ) 184 184 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 185 185 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) … … 204 204 END SELECT 205 205 ! 206 DO_3D _10_10(2, jpkm1 )206 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 207 207 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 208 208 END_3D … … 213 213 214 214 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 215 DO_3D _00_00(2, jpkm1 )215 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 216 216 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 217 217 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) … … 234 234 ! Warning : after this step, en : right hand side of the matrix 235 235 236 DO_3D _00_00(2, jpkm1 )236 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 237 237 ! 238 238 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 330 330 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 331 331 ! ! Balance between the production and the dissipation terms 332 DO_2D _00_00332 DO_2D( 0, 0, 0, 0 ) 333 333 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 334 334 !! With thick deep ocean level thickness, this may be quite large, no ??? … … 348 348 ! 349 349 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 350 DO_2D _00_00350 DO_2D( 0, 0, 0, 0 ) 351 351 itop = mikt(ji,jj) ! k top w-point 352 352 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 366 366 CASE ( 1 ) ! Neumman boundary condition 367 367 ! 368 DO_2D _00_00368 DO_2D( 0, 0, 0, 0 ) 369 369 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 370 370 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 380 380 END_2D 381 381 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 382 DO_2D _00_00382 DO_2D( 0, 0, 0, 0 ) 383 383 itop = mikt(ji,jj) ! k top w-point 384 384 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 400 400 ! ---------------------------------------------------------- 401 401 ! 402 DO_3D _00_00(2, jpkm1 )402 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 403 403 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 404 404 END_3D 405 DO_3D _00_00(2, jpk )405 DO_3D( 0, 0, 0, 0, 2, jpk ) 406 406 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) 407 407 END_3D 408 DO_3DS _00_00(jpk-1, 2, -1 )408 DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 ) 409 409 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 410 410 END_3D … … 421 421 ! 422 422 CASE( 0 ) ! k-kl (Mellor-Yamada) 423 DO_3D _00_00(2, jpkm1 )423 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 424 424 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 425 425 END_3D 426 426 ! 427 427 CASE( 1 ) ! k-eps 428 DO_3D _00_00(2, jpkm1 )428 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 429 429 psi(ji,jj,jk) = eps(ji,jj,jk) 430 430 END_3D 431 431 ! 432 432 CASE( 2 ) ! k-w 433 DO_3D _00_00(2, jpkm1 )433 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 434 434 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 435 435 END_3D 436 436 ! 437 437 CASE( 3 ) ! generic 438 DO_3D _00_00(2, jpkm1 )438 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 439 439 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 440 440 END_3D … … 449 449 ! Warning : after this step, en : right hand side of the matrix 450 450 451 DO_3D _00_00(2, jpkm1 )451 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 452 452 ! 453 453 ! psi / k … … 546 546 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 547 547 ! ! Balance between the production and the dissipation terms 548 DO_2D _00_00548 DO_2D( 0, 0, 0, 0 ) 549 549 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 550 550 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 565 565 CASE ( 1 ) ! Neumman boundary condition 566 566 ! 567 DO_2D _00_00567 DO_2D( 0, 0, 0, 0 ) 568 568 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 569 569 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 593 593 ! ---------------- 594 594 ! 595 DO_3D _00_00(2, jpkm1 )595 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 596 596 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 597 597 END_3D 598 DO_3D _00_00(2, jpk )598 DO_3D( 0, 0, 0, 0, 2, jpk ) 599 599 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) 600 600 END_3D 601 DO_3DS _00_00(jpk-1, 2, -1 )601 DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 ) 602 602 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 603 603 END_3D … … 609 609 ! 610 610 CASE( 0 ) ! k-kl (Mellor-Yamada) 611 DO_3D _00_00(1, jpkm1 )611 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 612 612 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) 613 613 END_3D 614 614 ! 615 615 CASE( 1 ) ! k-eps 616 DO_3D _00_00(1, jpkm1 )616 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 617 617 eps(ji,jj,jk) = psi(ji,jj,jk) 618 618 END_3D 619 619 ! 620 620 CASE( 2 ) ! k-w 621 DO_3D _00_00(1, jpkm1 )621 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 622 622 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 623 623 END_3D … … 627 627 zex1 = ( 1.5_wp + rmm/rnn ) 628 628 zex2 = -1._wp / rnn 629 DO_3D _00_00(1, jpkm1 )629 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 630 630 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 631 631 END_3D … … 635 635 ! Limit dissipation rate under stable stratification 636 636 ! -------------------------------------------------- 637 DO_3D _00_00(1, jpkm1 )637 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 638 638 ! limitation 639 639 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) … … 651 651 ! 652 652 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 653 DO_3D _00_00(2, jpkm1 )653 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 654 654 ! zcof = l²/q² 655 655 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) … … 668 668 ! 669 669 CASE ( 2, 3 ) ! Canuto stability functions 670 DO_3D _00_00(2, jpkm1 )670 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 671 671 ! zcof = l²/q² 672 672 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) … … 700 700 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 701 701 zstm(:,:,jpk) = 0. 702 DO_2D _00_00702 DO_2D( 0, 0, 0, 0 ) 703 703 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 704 704 END_2D … … 715 715 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 716 716 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 717 DO_3D _00_00(1, jpk )717 DO_3D( 0, 0, 0, 0, 1, jpk ) 718 718 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 719 719 zavt = zsqen * zstt(ji,jj,jk)
Note: See TracChangeset
for help on using the changeset viewer.