- 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/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
Note: See TracChangeset
for help on using the changeset viewer.
