- Timestamp:
- 2017-07-04T17:46:48+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r8055 r8279 159 159 ! Preliminary computing 160 160 161 ustars2 = 0._wp ; ustarb2= 0._wp ; psi = 0._wp ; zwall_psi = 0._wp161 ustars2(WRK_2D) = 0._wp ; ustarb2(WRK_2D) = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 162 162 163 163 164 164 ! Compute surface and bottom friction at T-points 165 DO jj = k_Jstr, k_Jend166 DO ji = k_Istr, k_Iend165 DO jj = tnldj, tnlej 166 DO ji = tnldi, tnlei 167 167 ! 168 168 ! surface friction … … 195 195 ! 196 196 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 197 DO jj = k_Jstr, k_Jend198 DO ji = k_Istr, k_Iend197 DO jj = tnldj, tnlej 198 DO ji = tnldi, tnlei 199 199 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 200 200 END DO … … 208 208 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 209 209 DO jk = 2, jpkm1 210 DO jj = k_Jstr, k_Jend211 DO ji = k_Istr, k_Iend210 DO jj = tnldj, tnlej 211 DO ji = tnldi, tnlei 212 212 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 213 213 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) … … 233 233 234 234 DO jk = 2, jpkm1 235 DO jj = k_Jstr, k_Jend236 DO ji = k_Istr, k_Iend235 DO jj = tnldj, tnlej 236 DO ji = tnldi, tnlei 237 237 ! 238 238 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 329 329 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 330 330 ! ! Balance between the production and the dissipation terms 331 DO jj = k_Jstr, k_Jend332 DO ji = k_Istr, k_Iend331 DO jj = tnldj, tnlej 332 DO ji = tnldi, tnlei 333 333 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 334 334 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 350 350 CASE ( 1 ) ! Neumman boundary condition 351 351 ! 352 DO jj = k_Jstr, k_Jend353 DO ji = k_Istr, k_Iend352 DO jj = tnldj, tnlej 353 DO ji = tnldi, tnlei 354 354 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 355 355 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 373 373 ! 374 374 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 375 DO jj = k_Jstr, k_Jend376 DO ji = k_Istr, k_Iend375 DO jj = tnldj, tnlej 376 DO ji = tnldi, tnlei 377 377 z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1) 378 378 END DO … … 380 380 END DO 381 381 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 382 DO jj = k_Jstr, k_Jend383 DO ji = k_Istr, k_Iend382 DO jj = tnldj, tnlej 383 DO ji = tnldi, tnlei 384 384 z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1) 385 385 END DO … … 387 387 END DO 388 388 DO jk = jpk-1, 2, -1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 389 DO jj = k_Jstr, k_Jend390 DO ji = k_Istr, k_Iend389 DO jj = tnldj, tnlej 390 DO ji = tnldi, tnlei 391 391 en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk) 392 392 END DO … … 406 406 CASE( 0 ) ! k-kl (Mellor-Yamada) 407 407 DO jk = 2, jpkm1 408 DO jj = k_Jstr, k_Jend409 DO ji = k_Istr, k_Iend408 DO jj = tnldj, tnlej 409 DO ji = tnldi, tnlei 410 410 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 411 411 END DO … … 415 415 CASE( 1 ) ! k-eps 416 416 DO jk = 2, jpkm1 417 DO jj = k_Jstr, k_Jend418 DO ji = k_Istr, k_Iend417 DO jj = tnldj, tnlej 418 DO ji = tnldi, tnlei 419 419 psi(ji,jj,jk) = eps(ji,jj,jk) 420 420 END DO … … 424 424 CASE( 2 ) ! k-w 425 425 DO jk = 2, jpkm1 426 DO jj = k_Jstr, k_Jend427 DO ji = k_Istr, k_Iend426 DO jj = tnldj, tnlej 427 DO ji = tnldi, tnlei 428 428 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 429 429 END DO … … 433 433 CASE( 3 ) ! generic 434 434 DO jk = 2, jpkm1 435 DO jj = k_Jstr, k_Jend436 DO ji = k_Istr, k_Iend435 DO jj = tnldj, tnlej 436 DO ji = tnldi, tnlei 437 437 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 438 438 END DO … … 450 450 451 451 DO jk = 2, jpkm1 452 DO jj = k_Jstr, k_Jend453 DO ji = k_Istr, k_Iend452 DO jj = tnldj, tnlej 453 DO ji = tnldi, tnlei 454 454 ! 455 455 ! psi / k … … 548 548 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0 549 549 ! ! Balance between the production and the dissipation terms 550 DO jj = k_Jstr, k_Jend551 DO ji = k_Istr, k_Iend550 DO jj = tnldj, tnlej 551 DO ji = tnldi, tnlei 552 552 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 553 553 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 569 569 CASE ( 1 ) ! Neumman boundary condition 570 570 ! 571 DO jj = k_Jstr, k_Jend572 DO ji = k_Istr, k_Iend571 DO jj = tnldj, tnlej 572 DO ji = tnldi, tnlei 573 573 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 574 574 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 600 600 ! 601 601 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 602 DO jj = k_Jstr, k_Jend603 DO ji = k_Istr, k_Iend602 DO jj = tnldj, tnlej 603 DO ji = tnldi, tnlei 604 604 z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1) 605 605 END DO … … 607 607 END DO 608 608 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 609 DO jj = k_Jstr, k_Jend610 DO ji = k_Istr, k_Iend609 DO jj = tnldj, tnlej 610 DO ji = tnldi, tnlei 611 611 z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1) 612 612 END DO … … 614 614 END DO 615 615 DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 616 DO jj = k_Jstr, k_Jend617 DO ji = k_Istr, k_Iend616 DO jj = tnldj, tnlej 617 DO ji = tnldi, tnlei 618 618 psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk) 619 619 END DO … … 628 628 CASE( 0 ) ! k-kl (Mellor-Yamada) 629 629 DO jk = 1, jpkm1 630 DO jj = k_Jstr, k_Jend631 DO ji = k_Istr, k_Iend630 DO jj = tnldj, tnlej 631 DO ji = tnldi, tnlei 632 632 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) 633 633 END DO … … 637 637 CASE( 1 ) ! k-eps 638 638 DO jk = 1, jpkm1 639 DO jj = k_Jstr, k_Jend640 DO ji = k_Istr, k_Iend639 DO jj = tnldj, tnlej 640 DO ji = tnldi, tnlei 641 641 eps(ji,jj,jk) = psi(ji,jj,jk) 642 642 END DO … … 646 646 CASE( 2 ) ! k-w 647 647 DO jk = 1, jpkm1 648 DO jj = k_Jstr, k_Jend649 DO ji = k_Istr, k_Iend648 DO jj = tnldj, tnlej 649 DO ji = tnldi, tnlei 650 650 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 651 651 END DO … … 658 658 zex2 = -1._wp / rnn 659 659 DO jk = 1, jpkm1 660 DO jj = k_Jstr, k_Jend661 DO ji = k_Istr, k_Iend660 DO jj = tnldj, tnlej 661 DO ji = tnldi, tnlei 662 662 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 663 663 END DO … … 670 670 ! -------------------------------------------------- 671 671 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 672 DO jj = k_Jstr, k_Jend673 DO ji = k_Istr, k_Iend672 DO jj = tnldj, tnlej 673 DO ji = tnldi, tnlei 674 674 ! limitation 675 675 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) … … 690 690 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 691 691 DO jk = 2, jpkm1 692 DO jj = k_Jstr, k_Jend693 DO ji = k_Istr, k_Iend692 DO jj = tnldj, tnlej 693 DO ji = tnldi, tnlei 694 694 ! zcof = l²/q² 695 695 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) … … 711 711 CASE ( 2, 3 ) ! Canuto stability functions 712 712 DO jk = 2, jpkm1 713 DO jj = k_Jstr, k_Jend714 DO ji = k_Istr, k_Iend713 DO jj = tnldj, tnlej 714 DO ji = tnldi, tnlei 715 715 ! zcof = l²/q² 716 716 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) … … 745 745 746 746 !!gm should be done for ISF (top boundary cond.) 747 DO jj = k_Jstr, k_Jend748 DO ji = k_Istr, k_Iend747 DO jj = tnldj, tnlej 748 DO ji = tnldi, tnlei 749 749 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 750 750 END DO … … 754 754 ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used 755 755 DO jk = 1, jpk 756 DO jj = k_Jstr, k_Jend757 DO ji = k_Istr, k_Iend756 DO jj = tnldj, tnlej 757 DO ji = tnldi, tnlei 758 758 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 759 759 zav = zsqen * zstt(ji,jj,jk)
Note: See TracChangeset
for help on using the changeset viewer.
