Changeset 8279 for branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO
- Timestamp:
- 2017-07-04T17:46:48+02:00 (7 years ago)
- Location:
- branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 10 edited
-
DOM/dom_omp.F90 (added)
-
ZDF/zdfddm.F90 (modified) (3 diffs)
-
ZDF/zdfevd.F90 (modified) (1 diff)
-
ZDF/zdfgls.F90 (modified) (28 diffs)
-
ZDF/zdfiwm.F90 (modified) (12 diffs)
-
ZDF/zdfphy.F90 (modified) (8 diffs)
-
ZDF/zdfric.F90 (modified) (3 diffs)
-
ZDF/zdfsh2.F90 (modified) (2 diffs)
-
ZDF/zdfswm.F90 (modified) (1 diff)
-
ZDF/zdftke.F90 (modified) (25 diffs)
-
domain_substitute.h90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r8055 r8279 95 95 !!gm and many acces in memory 96 96 97 DO jj = k_Jstr, k_Jend!== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==!98 DO ji = k_Jstr, k_Iend97 DO jj = tnldj, tnlej !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==! 98 DO ji = tnldi, tnlei 99 99 zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) & 100 100 !!gm please, use e3w_n below … … 113 113 END DO 114 114 115 DO jj = k_Jstr, k_Jend!== indicators ==!116 DO ji = k_Jstr, k_Iend115 DO jj = tnldj, tnlej !== indicators ==! 116 DO ji = tnldi, tnlei 117 117 ! stability indicator: msks=1 if rn2>0; 0 elsewhere 118 118 IF( rn2(ji,jj,jk) + 1.e-12 <= 0. ) THEN ; zmsks(ji,jj) = 0._wp … … 145 145 ! ------------------ 146 146 ! Constant eddy coefficient: reset to the background value 147 DO jj = k_Jstr, k_Jend148 DO ji = k_Jstr, k_Iend147 DO jj = tnldj, tnlej 148 DO ji = tnldi, tnlei 149 149 zinr = 1._wp / zrau(ji,jj) 150 150 ! salt fingering -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r8055 r8279 98 98 ! 99 99 DO jk = 1, jpkm1 100 DO jj = k_Jstr, k_Jend101 DO ji = k_Jstr, k_Iend100 DO jj = tnldj, tnlej 101 DO ji = tnldi, tnlei 102 102 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 103 103 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) -
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) -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
r8055 r8279 150 150 ! !* Critical slope mixing: distribute energy over the time-varying ocean depth, 151 151 ! using an exponential decay from the seafloor. 152 DO jj = k_Jstr, k_Jend153 DO ji = k_Istr, k_Iend! part independent of the level152 DO jj = tnldj, tnlej 153 DO ji = tnldi, tnlei ! part independent of the level 154 154 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 155 155 zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) … … 158 158 END DO 159 159 DO jk = 2, jpkm1 ! complete with the level-dependent part 160 DO jj = k_Jstr, k_Jend161 DO ji = k_Istr, k_Iend160 DO jj = tnldj, tnlej 161 DO ji = tnldi, tnlei 162 162 emix_iwm(ji,jj,jk) = zfact(ji,jj) * wmask(ji,jj,jk) & 163 163 & * ( EXP( ( gde3w_n(ji,jj,jk ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) & … … 185 185 ! 186 186 WHERE( zfact(WRK_2D) /= 0 ) zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 187 ! DO jj = k_Jstr, k_Jend188 ! DO ji = k_Istr, k_Iend187 ! DO jj = tnldj, tnlej 188 ! DO ji = tnldi, tnlei 189 189 ! IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 190 190 ! END DO … … 203 203 ! 204 204 WHERE( zfact(WRK_2D) /= 0 ) zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 205 ! DO jj = k_Jstr, k_Jend206 ! DO ji = k_Istr, k_Iend205 ! DO jj = tnldj, tnlej 206 ! DO ji = tnldi, tnlei 207 207 ! IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 208 208 ! END DO … … 224 224 ! zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT( MAX( 0._wp, rn2(WRK_2D,jk) ) ) 225 225 ! zwkb (WRK_2D,jk) = zfact(WRK_2D) 226 DO jj = k_Jstr, k_Jend227 DO ji = k_Istr, k_Iend226 DO jj = tnldj, tnlej 227 DO ji = tnldi, tnlei 228 228 zztemp = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) 229 229 zfact(ji,jj) = zztemp … … 234 234 ! 235 235 DO jk = 2, jpkm1 236 DO jj = k_Jstr, k_Jend237 DO ji = k_Istr, k_Iend236 DO jj = tnldj, tnlej 237 DO ji = tnldi, tnlei 238 238 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 239 239 & * wmask(ji,jj,jk) / zfact(ji,jj) … … 245 245 zweight(WRK_3D) = 0._wp 246 246 DO jk = 2, jpkm1 247 zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm( :,:)&247 zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(WRK_2D) & 248 248 & * ( EXP( -zwkb(WRK_2D,jk ) / hbot_iwm(WRK_2D) ) & 249 249 & - EXP( -zwkb(WRK_2D,jk-1) / hbot_iwm(WRK_2D) ) ) … … 256 256 ! 257 257 WHERE( zfact(WRK_2D) /= 0 ) zfact(WRK_2D) = ebot_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 258 ! DO jj = k_Jstr, k_Jend259 ! DO ji = k_Istr, k_Iend258 ! DO jj = tnldj, tnlej 259 ! DO ji = tnldi, tnlei 260 260 ! IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 261 261 ! END DO … … 290 290 IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the 291 291 DO jk = 2, jpkm1 ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 292 DO jj = k_Jstr, k_Jend293 DO ji = k_Istr, k_Iend292 DO jj = tnldj, tnlej 293 DO ji = tnldi, tnlei 294 294 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 295 295 zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) … … 310 310 !!gm used of glosum 3D.... 311 311 DO jk = 2, jpkm1 312 DO jj = k_Jstr, k_Jend313 DO ji = k_Istr, k_Iend312 DO jj = tnldj, tnlej 313 DO ji = tnldi, tnlei 314 314 zztemp = zztemp + e3w_n(ji,jj,jk) * e1e2t(ji,jj) & 315 315 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) … … 317 317 END DO 318 318 END DO 319 IF( lk_mpp ) CALL mpp_sum( zztemp )319 ! IF( lk_mpp ) CALL mpp_sum( zztemp ) 320 320 zztemp = rau0 * zztemp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 321 321 ! … … 335 335 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature 336 336 DO jk = 2, jpkm1 ! Calculate S/T diffusivity ratio as a function of Reb 337 DO jj = k_Jstr, k_Jend338 DO ji = k_Istr, k_Iend337 DO jj = tnldj, tnlej 338 DO ji = tnldi, tnlei 339 339 zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp * & 340 340 & TANH( 0.92_wp * ( LOG10( MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 ) ) - 0.60_wp ) ) & -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
r8055 r8279 32 32 USE lbclnk ! lateral boundary conditions 33 33 USE lib_mpp ! distribued memory computing 34 #if defined key_omp 35 USE dom_omp ! omp threads decomposition of MPI subdomain 36 #endif 34 37 35 38 IMPLICIT NONE … … 47 50 48 51 LOGICAL :: l_zdfsh2 ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 52 INTEGER :: ARG_2D 49 53 50 54 !! * Substitutions … … 216 220 INTEGER, INTENT(in) :: kt ! ocean time-step index 217 221 ! 218 INTEGER :: ji, jj,jk ! dummy loop indice222 INTEGER :: jk ! dummy loop indice 219 223 !!OMP REAL(wp), DIMENSION(WRK_3D) :: zsh2 ! shear production 220 224 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zsh2 ! shear production 221 225 !! --------------------------------------------------------------------- 222 226 ! 223 !!OMP ===>>>> Open-MP to be defined elsewhere 224 ! 225 INTEGER :: ARG_2D 226 ! 227 ! 228 IF (kt==nn_it000) THEN 229 #ifdef key_omp 230 CALL tile_dec_2D 231 #else 227 232 #if defined key_vectopt_loop 228 k_Istr = 1 ; k_Iend= jpi233 tnldi = 1 ; tnlei = jpi 229 234 #else 230 k_Istr = 2 ; k_Iend= jpim1235 tnldi = 2 ; tnlei = jpim1 231 236 #endif 232 k_Jstr = 2 ; k_Jend = jpjm1 233 ! 234 !!OMP <<<<=== end 235 237 tnldj = 2 ; tnlej = jpjm1 238 #endif 239 END IF 240 ! 241 CALL zdf_bfr( kt ) !* bottom friction (if quadratic) 242 ! 243 !$OMP PARALLEL PRIVATE(zsh2) 236 244 ALLOCATE( zsh2(WRK_3D) ) 237 238 239 !240 CALL zdf_bfr( kt ) !* bottom friction (if quadratic)241 !242 243 245 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 244 245 246 IF( l_zdfsh2 ) & !* shear production at w-points (energy conserving form) 246 247 CALL zdf_sh2( ARG_2D, ub, vb, un, vn, avm_k, & ! <<== in fields … … 248 249 249 250 ! 251 !$OMP BARRIER 250 252 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 251 253 CASE( np_RIC ) ; CALL zdf_ric( ARG_2D, kt, gdept_n, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz … … 257 259 !!gm avm_k(WRK_3D) = rn_avm0 * wmask(WRK_3D) 258 260 END SELECT 261 ! 262 DEALLOCATE( zsh2 ) 259 263 ! 260 264 ! !== ocean Kz ==! (avt, avs, avm) … … 271 275 ! 272 276 IF( ln_zdfevd ) CALL zdf_evd( ARG_2D, kt, avm, avt ) !* convection: enhanced vertical eddy diffusivity 273 ! 277 274 278 ! !* double diffusive mixing 275 279 IF( ln_zdfddm ) THEN ! update avt and compute avs … … 282 286 IF( ln_zdfswm ) CALL zdf_swm( ARG_2D, kt, avm, avt, avs ) ! surface wave (Qiao et al. 2004) 283 287 IF( ln_zdfiwm ) CALL zdf_iwm( ARG_2D, kt, avm, avt, avs ) ! internal wave (de Lavergne et al 2017) 284 288 !$OMP END PARALLEL 289 ! 285 290 286 291 ! !* Lateral boundary conditions (sign unchanged) … … 301 306 ENDIF 302 307 ! 303 DEALLOCATE( zsh2 )304 !305 308 END SUBROUTINE zdf_phy 306 309 -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r8056 r8279 164 164 ! !== avm and avt = F(Richardson number) ==! 165 165 DO jk = 2, jpkm1 166 DO jj = k_Jstr, k_Jend167 DO ji = k_Jstr, k_Iend! coefficient = F(richardson number) (avm-weighted Ri)166 DO jj = tnldj, tnlej 167 DO ji = tnldi, tnlei ! coefficient = F(richardson number) (avm-weighted Ri) 168 168 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , p_avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) ) 169 169 zav = rn_avmri * zcfRi**nn_ric … … 180 180 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 181 181 ! 182 DO jj = k_Jstr, k_Jend!* Ekman depth183 DO ji = k_Jstr, k_Iend182 DO jj = tnldj, tnlej !* Ekman depth 183 DO ji = tnldi, tnlei 184 184 zustar = SQRT( taum(ji,jj) * r1_rau0 ) 185 185 zed = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth … … 188 188 END DO 189 189 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer 190 DO jj = k_Jstr, k_Jend191 DO ji = k_Jstr, k_Iend190 DO jj = tnldj, tnlej 191 DO ji = tnldi, tnlei 192 192 IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 193 193 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90
r8055 r8279 62 62 ! 63 63 DO jk = 2, jpkm1 64 DO jj = k_Jstr-1, k_Jend!* 2 x shear production at uw- and vw-points (energy conserving form)65 DO ji = k_Istr-1, k_Iend64 DO jj = tnldj-1, tnlej !* 2 x shear production at uw- and vw-points (energy conserving form) 65 DO ji = tnldi-1, tnlei 66 66 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 67 67 & * ( pun(ji,jj,jk-1) - pun(ji,jj,jk) ) & … … 72 72 END DO 73 73 END DO 74 DO jj = k_Jstr, k_Jend!* shear production at w-point75 DO ji = k_Jstr, k_Iend! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)74 DO jj = tnldj, tnlej !* shear production at w-point 75 DO ji = tnldi, tnlei ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 76 76 p_sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) & 77 77 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) ) -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90
r8055 r8279 64 64 zcoef = 1._wp * 0.353553_wp 65 65 DO jk = 2, jpkm1 66 DO jj = k_Jstr, k_Jend67 DO ji = k_Istr, k_Iend66 DO jj = tnldj, tnlej 67 DO ji = tnldi, tnlei 68 68 zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk) 69 69 ! -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r8056 r8279 220 220 ! 221 221 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 222 222 223 ! 223 224 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 231 232 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 232 233 IF ( ln_isfcav ) THEN 233 DO jj = k_Jstr, k_Jend! en(mikt(ji,jj)) = rn_emin234 DO ji = k_Istr, k_Iend234 DO jj = tnldj, tnlej ! en(mikt(ji,jj)) = rn_emin 235 DO ji = tnldi, tnlei 235 236 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 236 237 END DO 237 238 END DO 238 239 END IF 239 DO jj = k_Jstr, k_Jend! en(1) = rn_ebb taum / rau0 (min value rn_emin0)240 DO ji = k_Istr, k_Iend240 DO jj = tnldj, tnlej ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 241 DO ji = tnldi, tnlei 241 242 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 242 243 END DO … … 254 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 255 256 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 256 !! DO jj = k_Jstr, k_Jend257 !! DO ji = k_Jstr, k_Iend257 !! DO jj = tnldj, tnlej 258 !! DO ji = tnldj, tnlei 258 259 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 259 260 !! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) ) … … 281 282 imlc(WRK_2D) = mbkt(WRK_2D) + 1 ! Initialization to the number of w ocean point (=2 over land) 282 283 DO jk = jpkm1, 2, -1 283 DO jj = k_Jstr, k_Jend! Last w-level at which zpelc>=0.5*us*us284 DO ji = k_Istr, k_Iend! with us=0.016*wind(starting from jpk-1)284 DO jj = tnldj, tnlej ! Last w-level at which zpelc>=0.5*us*us 285 DO ji = tnldi, tnlei ! with us=0.016*wind(starting from jpk-1) 285 286 zus = zcof * taum(ji,jj) 286 287 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk … … 289 290 END DO 290 291 ! ! finite LC depth 291 DO jj = k_Jstr, k_Jend292 DO ji = k_Istr, k_Iend292 DO jj = tnldj, tnlej 293 DO ji = tnldi, tnlei 293 294 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 294 295 END DO … … 296 297 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 297 298 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 298 DO jj = k_Jstr, k_Jend299 DO ji = k_Istr, k_Iend299 DO jj = tnldj, tnlej 300 DO ji = tnldi, tnlei 300 301 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 301 302 ! ! vertical velocity due to LC … … 320 321 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 321 322 DO jk = 2, jpkm1 322 DO jj = k_Jstr, k_Jend323 DO ji = k_Istr, k_Iend323 DO jj = tnldj, tnlej 324 DO ji = tnldi, tnlei 324 325 ! ! local Richardson number 325 326 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) … … 332 333 ! 333 334 DO jk = 2, jpkm1 !* Matrix and right hand side in en 334 DO jj = k_Jstr, k_Jend335 DO ji = k_Istr, k_Iend335 DO jj = tnldj, tnlej 336 DO ji = tnldi, tnlei 336 337 zcof = zfact1 * tmask(ji,jj,jk) 337 338 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical … … 356 357 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 357 358 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 358 DO jj = k_Jstr, k_Jend359 DO ji = k_Istr, k_Iend359 DO jj = tnldj, tnlej 360 DO ji = tnldi, tnlei 360 361 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 361 362 END DO 362 363 END DO 363 364 END DO 364 DO jj = k_Jstr, k_Jend! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1365 DO ji = k_Istr, k_Iend365 DO jj = tnldj, tnlej ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 366 DO ji = tnldi, tnlei 366 367 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 367 368 END DO 368 369 END DO 369 370 DO jk = 3, jpkm1 370 DO jj = k_Jstr, k_Jend371 DO ji = k_Istr, k_Iend371 DO jj = tnldj, tnlej 372 DO ji = tnldi, tnlei 372 373 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) 373 374 END DO 374 375 END DO 375 376 END DO 376 DO jj = k_Jstr, k_Jend! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk377 DO ji = k_Istr, k_Iend377 DO jj = tnldj, tnlej ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 378 DO ji = tnldi, tnlei 378 379 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 379 380 END DO 380 381 END DO 381 382 DO jk = jpk-2, 2, -1 382 DO jj = k_Jstr, k_Jend383 DO ji = k_Istr, k_Iend383 DO jj = tnldj, tnlej 384 DO ji = tnldi, tnlei 384 385 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 385 386 END DO … … 387 388 END DO 388 389 DO jk = 2, jpkm1 ! set the minimum value of tke 389 DO jj = k_Jstr, k_Jend390 DO ji = k_Istr, k_Iend390 DO jj = tnldj, tnlej 391 DO ji = tnldi, tnlei 391 392 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 392 393 END DO … … 402 403 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 403 404 DO jk = 2, jpkm1 404 DO jj = k_Jstr, k_Jend405 DO ji = k_Istr, k_Iend405 DO jj = tnldj, tnlej 406 DO ji = tnldi, tnlei 406 407 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 407 408 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) … … 410 411 END DO 411 412 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 412 DO jj = k_Jstr, k_Jend413 DO ji = k_Istr, k_Iend413 DO jj = tnldj, tnlej 414 DO ji = tnldi, tnlei 414 415 jk = nmln(ji,jj) 415 416 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & … … 419 420 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 420 421 DO jk = 2, jpkm1 421 DO jj = k_Jstr, k_Jend422 DO ji = k_Istr, k_Iend422 DO jj = tnldj, tnlej 423 DO ji = tnldi, tnlei 423 424 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 424 425 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 433 434 ENDIF 434 435 ! 436 435 437 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') 436 438 ! … … 498 500 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 499 501 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 500 DO jj = k_Jstr, k_Jend501 DO ji = k_Istr, k_Iend502 DO jj = tnldj, tnlej 503 DO ji = tnldi, tnlei 502 504 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 503 505 END DO … … 508 510 ! 509 511 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 510 DO jj = k_Jstr, k_Jend511 DO ji = k_Istr, k_Iend512 DO jj = tnldj, tnlej 513 DO ji = tnldi, tnlei 512 514 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 513 515 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 527 529 CASE ( 0 ) ! bounded by the distance to surface and bottom 528 530 DO jk = 2, jpkm1 529 DO jj = k_Jstr, k_Jend530 DO ji = k_Istr, k_Iend531 DO jj = tnldj, tnlej 532 DO ji = tnldi, tnlei 531 533 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 532 534 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) … … 540 542 CASE ( 1 ) ! bounded by the vertical scale factor 541 543 DO jk = 2, jpkm1 542 DO jj = k_Jstr, k_Jend543 DO ji = k_Istr, k_Iend544 DO jj = tnldj, tnlej 545 DO ji = tnldi, tnlei 544 546 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 545 547 zmxlm(ji,jj,jk) = zemxl … … 551 553 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 552 554 DO jk = 2, jpkm1 ! from the surface to the bottom : 553 DO jj = k_Jstr, k_Jend554 DO ji = k_Istr, k_Iend555 DO jj = tnldj, tnlej 556 DO ji = tnldi, tnlei 555 557 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 556 558 END DO … … 558 560 END DO 559 561 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 560 DO jj = k_Jstr, k_Jend561 DO ji = k_Istr, k_Iend562 DO jj = tnldj, tnlej 563 DO ji = tnldi, tnlei 562 564 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 563 565 zmxlm(ji,jj,jk) = zemxl … … 569 571 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 570 572 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 571 DO jj = k_Jstr, k_Jend572 DO ji = k_Istr, k_Iend573 DO jj = tnldj, tnlej 574 DO ji = tnldi, tnlei 573 575 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 574 576 END DO … … 576 578 END DO 577 579 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 578 DO jj = k_Jstr, k_Jend579 DO ji = k_Istr, k_Iend580 DO jj = tnldj, tnlej 581 DO ji = tnldi, tnlei 580 582 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 581 583 END DO … … 583 585 END DO 584 586 DO jk = 2, jpkm1 585 DO jj = k_Jstr, k_Jend586 DO ji = k_Istr, k_Iend587 DO jj = tnldj, tnlej 588 DO ji = tnldi, tnlei 587 589 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 588 590 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) … … 600 602 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 601 603 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 602 DO jj = k_Jstr, k_Jend603 DO ji = k_Istr, k_Iend604 DO jj = tnldj, tnlej 605 DO ji = tnldi, tnlei 604 606 zsqen = SQRT( en(ji,jj,jk) ) 605 607 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 614 616 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 615 617 DO jk = 2, jpkm1 616 DO jj = k_Jstr, k_Jend617 DO ji = k_Istr, k_Iend618 DO jj = tnldj, tnlej 619 DO ji = tnldi, tnlei 618 620 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 619 621 END DO -
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/domain_substitute.h90
r8055 r8279 4 4 !! ** purpose : local Sub-Domain: arguments, dimension & semi-column loops 5 5 !!---------------------------------------------------------------------- 6 #define ARG_2D k_Istr, k_Iend, k_Jstr, k_Jend6 #define ARG_2D tnldi, tnlei, tnldj, tnlej 7 7 8 #define WRK_2D k_Istr:k_Iend,k_Jstr:k_Jend9 #define WRK_2De(s,e) k_Istr s:k_Iend e,k_Jstr s:k_Jende8 #define WRK_2D tnldi:tnlei,tnldj:tnlej 9 #define WRK_2De(s,e) tnldi s:tnlei e,tnldj s:tnlej e 10 10 11 #define WRK_3D k_Istr:k_Iend,k_Jstr:k_Jend,1:jpk12 #define WRK_3De(s,e,t,b) k_Istr s:k_Iend e,k_Jstr s:k_Jende,1 t:jpk t11 #define WRK_3D tnldi:tnlei,tnldj:tnlej,1:jpk 12 #define WRK_3De(s,e,t,b) tnldi s:tnlei e,tnldj s:tnlej e,1 t:jpk t
Note: See TracChangeset
for help on using the changeset viewer.
