- Timestamp:
- 2020-09-15T12:56:56+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdftke.F90
r13469 r13470 231 231 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 232 232 ! 233 DO_2D _00_00233 DO_2D( 0, 0, 0, 0 ) 234 234 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 235 235 !! one way around would be to increase zbbirau … … 249 249 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 250 250 ! 251 DO_2D _00_00251 DO_2D( 0, 0, 0, 0 ) 252 252 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 253 253 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) … … 258 258 END_2D 259 259 IF( ln_isfcav ) THEN ! top friction 260 DO_2D _00_00260 DO_2D( 0, 0, 0, 0 ) 261 261 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 262 262 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) … … 283 283 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 284 284 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 285 DO_3D _11_11(jpkm1, 2, -1 )285 DO_3D( 1, 1, 1, 1, jpkm1, 2, -1 ) 286 286 zus = zcof * taum(ji,jj) 287 287 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 288 288 END_3D 289 289 ! ! finite LC depth 290 DO_2D _11_11290 DO_2D( 1, 1, 1, 1 ) 291 291 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 292 292 END_2D 293 293 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 294 DO_2D _00_00294 DO_2D( 0, 0, 0, 0 ) 295 295 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 296 296 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 297 297 END_2D 298 DO_3D _00_00(2, jpkm1 )298 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 299 299 IF ( zus3(ji,jj) /= 0._wp ) THEN 300 300 ! vertical velocity due to LC … … 318 318 ! 319 319 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 320 DO_3D _00_00(2, jpkm1 )320 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 321 321 ! ! local Richardson number 322 322 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) … … 326 326 ENDIF 327 327 ! 328 DO_3D _00_00(2, jpkm1 )328 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 329 329 zcof = zfact1 * tmask(ji,jj,jk) 330 330 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical … … 346 346 END_3D 347 347 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 348 DO_3D _00_00(3, jpkm1 )348 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 349 349 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 350 350 END_3D 351 DO_2D _00_00351 DO_2D( 0, 0, 0, 0 ) 352 352 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 353 353 END_2D 354 DO_3D _00_00(3, jpkm1 )354 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 355 355 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) 356 356 END_3D 357 DO_2D _00_00357 DO_2D( 0, 0, 0, 0 ) 358 358 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 359 359 END_2D 360 DO_3D _00_00(jpk-2, 2, -1 )360 DO_3D( 0, 0, 0, 0, jpk-2, 2, -1 ) 361 361 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 362 362 END_3D 363 DO_3D _00_00(2, jpkm1 )363 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 364 364 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 365 365 END_3D … … 373 373 374 374 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 375 DO_3D _00_00(2, jpkm1 )375 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 376 376 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 377 377 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 378 378 END_3D 379 379 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 380 DO_2D _00_00380 DO_2D( 0, 0, 0, 0 ) 381 381 jk = nmln(ji,jj) 382 382 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & … … 384 384 END_2D 385 385 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 386 DO_3D _00_00(2, jpkm1 )386 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 387 387 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 388 388 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 459 459 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 460 460 #if ! defined key_si3 && ! defined key_cice 461 DO_2D _00_00461 DO_2D( 0, 0, 0, 0 ) 462 462 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 463 463 END_2D … … 467 467 ! 468 468 CASE( 0 ) ! No scaling under sea-ice 469 DO_2D _00_00469 DO_2D( 0, 0, 0, 0 ) 470 470 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 471 471 END_2D 472 472 ! 473 473 CASE( 1 ) ! scaling with constant sea-ice thickness 474 DO_2D _00_00474 DO_2D( 0, 0, 0, 0 ) 475 475 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 476 476 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) … … 478 478 ! 479 479 CASE( 2 ) ! scaling with mean sea-ice thickness 480 DO_2D _00_00480 DO_2D( 0, 0, 0, 0 ) 481 481 #if defined key_si3 482 482 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 490 490 ! 491 491 CASE( 3 ) ! scaling with max sea-ice thickness 492 DO_2D _00_00492 DO_2D( 0, 0, 0, 0 ) 493 493 zmaxice = MAXVAL( h_i(ji,jj,:) ) 494 494 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & … … 499 499 #endif 500 500 ! 501 DO_2D _00_00501 DO_2D( 0, 0, 0, 0 ) 502 502 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 503 503 END_2D … … 507 507 ENDIF 508 508 ! 509 DO_3D _00_00(2, jpkm1 )509 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 510 510 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 511 511 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 522 522 ! where wmask = 0 set zmxlm == p_e3w 523 523 CASE ( 0 ) ! bounded by the distance to surface and bottom 524 DO_3D _00_00(2, jpkm1 )524 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 525 525 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 526 526 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) … … 531 531 ! 532 532 CASE ( 1 ) ! bounded by the vertical scale factor 533 DO_3D _00_00(2, jpkm1 )533 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 534 534 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 535 535 zmxlm(ji,jj,jk) = zemxl … … 538 538 ! 539 539 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 540 DO_3D _00_00(2, jpkm1 )540 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 541 541 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 542 542 END_3D 543 DO_3D _00_00(jpkm1, 2, -1 )543 DO_3D( 0, 0, 0, 0, jpkm1, 2, -1 ) 544 544 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 545 545 zmxlm(ji,jj,jk) = zemxl … … 548 548 ! 549 549 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 550 DO_3D _00_00(2, jpkm1 )550 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 551 551 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 552 552 END_3D 553 DO_3D _00_00(jpkm1, 2, -1 )553 DO_3D( 0, 0, 0, 0, jpkm1, 2, -1 ) 554 554 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 555 555 END_3D 556 DO_3D _00_00(2, jpkm1 )556 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 557 557 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 558 558 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) … … 566 566 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 567 567 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 568 DO_3D _00_00(1, jpkm1 )568 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 569 569 zsqen = SQRT( en(ji,jj,jk) ) 570 570 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 576 576 ! 577 577 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 578 DO_3D _00_00(2, jpkm1 )578 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 579 579 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 580 580 END_3D
Note: See TracChangeset
for help on using the changeset viewer.
