Changeset 1050 for trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2008-06-03T12:25:26+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
r899 r1050 18 18 !! 9.0 ! 02-08 (G. Madec) autotasking optimization 19 19 !! 9.0 ! 06-07 (S. Masson) distributed restart using iom 20 !! - ! 2005-07 (C. Ethe, G.Madec) : update TKE physics: 21 !! - tke penetration (wind steering) 22 !! - suface condition for tke & mixing length 23 !! - Langmuir cells 24 !! 3.0 ! 2007-11 (G. Madec) avtb_2d, nn_avtb_2d 20 25 !!---------------------------------------------------------------------- 21 26 #if defined key_zdftke || defined key_esopa … … 33 38 USE sbc_oce ! surface boundary condition: ocean 34 39 USE phycst ! physical constants 40 USE zdfmxl ! mixed layer 35 41 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 36 42 USE prtctl ! Print control … … 67 73 & emin0 = 1.e-4_wp , & !: surface minimum value of tke (m2/s2) 68 74 & ri_c = 2._wp / 9._wp !: critic Richardson number 75 ! !!! ** Namelist (namtke) ** 76 INTEGER :: nn_avtb_2d = 1 !: = 0/1 horizontal shape for avtb 77 INTEGER :: nn_etau = 0 !: = 0/1/2 tke depth penetration 78 INTEGER :: nn_htau = 0 !: = 0/1/2 type of tke profile of penetration 79 REAL(wp) :: rn_lmin = 0.4_wp !: surface min value of mixing turbulent length scale 80 REAL(wp) :: rn_efr = 1.0_wp !: fraction of TKE surface value which penetrates in the ocean 81 LOGICAL :: ln_lsfc = .FALSE. !: mixing length scale surface value as function of wind stress or not 82 LOGICAL :: ln_lc = .FALSE. !: Langmuir cells (LC) as a source term of TKE or not 83 REAL(wp) :: rn_lc = 0.15_wp !: coef to compute vertical velocity of LC 84 85 REAL(wp), DIMENSION (jpi,jpj) :: avtb_2d ! set in tke_init, for other modif than ice 69 86 70 87 # if defined key_c1d … … 148 165 & zmxlm => ta, & ! use ta as workspace 149 166 & zmxld => sa ! use sa as workspace 167 USE oce , ztkelc => va ! use va as workspace 150 168 ! 151 169 INTEGER, INTENT(in) :: kt ! ocean time step … … 160 178 & zsqen, zesh2, & ! 161 179 & zemxl, zemlm, zemlp 180 INTEGER , DIMENSION(jpi,jpj) :: imlc ! 2D workspace 181 REAL(wp) :: zraug, zus, zwlc, zind ! temporary scalar 182 REAL(wp), DIMENSION(jpi,jpj) :: zhtau ! 2D workspace 183 REAL(wp), DIMENSION(jpi,jpj) :: zhlc ! 2D workspace 184 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc ! 3D workspace 162 185 !!-------------------------------------------------------------------- 163 186 … … 165 188 166 189 ! ! Local constant initialization 167 zmlmin = 1.e-8190 zmlmin = 0.4 168 191 zbbrau = .5 * ebb / rau0 169 192 zfact1 = -.5 * rdt * efave … … 177 200 ! Buoyancy length scale: l=sqrt(2*e/n**2) 178 201 ! --------------------- 179 zmxlm(:,:, 1 ) = zmlmin ! surface set to the minimum value 180 zmxlm(:,:,jpk) = zmlmin ! bottom set to the minimum value 202 IF( ln_lsfc ) THEN ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 203 zmxlm(:,:,1) = 0.e0 204 zraug = 0.5 * 2.e5 / ( rau0 * grav ) 205 DO jj = 2, jpjm1 206 !CDIR NOVERRCHK 207 DO ji = fs_2, fs_jpim1 ! vector opt. 208 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 209 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 210 zmxlm(ji,jj,1 ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 211 ! set the surface minimum value to rn_lmin 212 zmxlm(ji,jj,1 ) = MAX( zmxlm(ji,jj,1) , rn_lmin ) 213 END DO 214 END DO 215 ELSE ! surface set to the minimum value 216 zmxlm(:,:,1) = rn_lmin 217 ENDIF 218 zmxlm(:,:,jpk) = rn_lmin ! bottom set to the minimum value 181 219 !CDIR NOVERRCHK 182 220 DO jk = 2, jpkm1 … … 193 231 ! Physical limits for the mixing length 194 232 ! ------------------------------------- 195 zmxld(:,:, 1 ) = zm lmin! surface set to the minimum value196 zmxld(:,:,jpk) = zmlmin! bottom set to the minimum value233 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 234 zmxld(:,:,jpk) = rn_lmin ! bottom set to the minimum value 197 235 198 236 SELECT CASE ( nmxl ) … … 277 315 278 316 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 279 ! II Tubulent kinetic energy time stepping 317 ! II TKE Langmuir circulation source term 318 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 319 IF( ln_lc ) THEN 320 ! 321 ! Computation of total energy produce by LC : cumulative sum over jk 322 zpelc(:,:,1) = MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 323 DO jk = 2, jpk 324 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 325 END DO 326 ! 327 ! Computation of finite Langmuir Circulation depth 328 ! Initialization to the number of w ocean point mbathy 329 imlc(:,:) = mbathy(:,:) 330 DO jk = jpkm1, 2, -1 331 DO jj = 1, jpj 332 DO ji = 1, jpi 333 ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1) 334 zus = 0.000128 * wndm(ji,jj) * wndm(ji,jj) 335 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 336 END DO 337 END DO 338 END DO 339 ! 340 ! finite LC depth 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 344 END DO 345 END DO 346 ! 347 # if defined key_c1d 348 hlc(:,:) = zhlc(:,:) * tmask(:,:,1) ! save finite Langmuir Circulation depth 349 # endif 350 ! 351 ! TKE Langmuir circulation source term 352 !CDIR NOVERRCHK 353 DO jk = 2, jpkm1 354 !CDIR NOVERRCHK 355 DO jj = 2, jpjm1 356 !CDIR NOVERRCHK 357 DO ji = fs_2, fs_jpim1 ! vector opt. 358 ! Stokes drift 359 zus = 0.016 * wndm(ji,jj) 360 ! computation of vertical velocity due to LC 361 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 362 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 363 ! TKE Langmuir circulation source term 364 ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 365 END DO 366 END DO 367 END DO 368 ! 369 ENDIF 370 371 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 372 ! III Tubulent kinetic energy time stepping 280 373 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 281 374 … … 290 383 zsqen = SQRT( en(ji,jj,jk) ) 291 384 zav = ediff * zmxlm(ji,jj,jk) * zsqen 292 avt (ji,jj,jk) = MAX( zav, avtb (jk) ) * tmask(ji,jj,jk)385 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 293 386 zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 294 387 zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) … … 309 402 zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 310 403 en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1) 311 zmxlm(ji,jj,1 ) = avmb(1) * tmask(ji,jj,1) 404 zav = ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(ji,jj,1) 405 zmxlm(ji,jj,1) = MAX( zav, avmb(1) ) * tmask(ji,jj,1) 406 avt (ji,jj,1) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1) 312 407 zmxlm(ji,jj,jpk) = 0.e0 313 408 END DO … … 347 442 zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 348 443 ! right hand side in en 349 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en (ji,jj,jk) & 350 & + rdt * zmxlm(ji,jj,jk) * zesh2 444 IF( .NOT. ln_lc ) THEN 445 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk) & 446 & + rdt * zmxlm (ji,jj,jk) * zesh2 447 ELSE 448 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk) & 449 & + rdt * zmxlm (ji,jj,jk) * zesh2 & 450 & + rdt * ztkelc(ji,jj,jk) 451 ENDIF 351 452 END DO 352 453 END DO … … 390 491 zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 391 492 ! right hand side in en 392 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en (ji,jj,jk) & 393 & + rdt * zmxlm(ji,jj,jk) * zesh2 394 ! save masked Prandlt number in zmxlm array 493 IF( .NOT. ln_lc ) THEN 494 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en (ji,jj,jk) & 495 & + rdt * zmxlm (ji,jj,jk) * zesh2 496 ELSE 497 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en (ji,jj,jk) & 498 & + rdt * zmxlm (ji,jj,jk) * zesh2 & 499 & + rdt * ztkelc(ji,jj,jk) 500 ENDIF 501 ! save masked Prandlt number in zmxld array 395 502 zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk) 396 503 END DO … … 454 561 END DO 455 562 END DO 563 564 ! Modify the value of TKE by an exponential assumption 565 ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) ) 566 567 SELECT CASE( nn_htau ) ! Choice of H value 568 ! 569 CASE( 0 ) 570 DO jj = 2, jpjm1 571 DO ji = fs_2, fs_jpim1 ! vector opt. 572 zhtau(ji,jj) = 10. 573 END DO 574 END DO 575 ! 576 CASE( 1 ) 577 DO jj = 2, jpjm1 578 DO ji = fs_2, fs_jpim1 ! vector opt. 579 zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 580 END DO 581 END DO 582 ! 583 CASE( 2 ) 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 586 zhtau(ji,jj) = MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 587 END DO 588 END DO 589 ! 590 END SELECT 591 592 IF( nn_etau == 1 ) THEN 593 DO jk = 2, jpkm1 594 DO jj = 2, jpjm1 595 DO ji = fs_2, fs_jpim1 ! vector opt. 596 en(ji,jj,jk) = en(ji,jj,jk) & 597 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 598 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 599 END DO 600 END DO 601 END DO 602 ELSEIF( nn_etau == 2 ) THEN ! only at the base of the mixed layer 603 DO jj = 2, jpjm1 604 DO ji = fs_2, fs_jpim1 ! vector opt. 605 jk = nmln(ji,jj) 606 en(ji,jj,jk) = en(ji,jj,jk) & 607 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 608 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 609 END DO 610 END DO 611 ENDIF 612 456 613 457 614 ! Lateral boundary conditions on ( avt, en ) (sign unchanged) … … 459 616 460 617 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 461 ! I II. Before vertical eddy vicosity and diffusivity coefficients618 ! IV. Before vertical eddy vicosity and diffusivity coefficients 462 619 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 463 620 … … 563 720 DO ji = fs_2, fs_jpim1 ! vector opt. 564 721 zpdl = zmxld(ji,jj,jk) 565 avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb (jk) ) * tmask(ji,jj,jk)722 avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 566 723 END DO 567 724 END DO … … 621 778 INTEGER :: jk ! dummy loop indices 622 779 # endif 623 780 !! 624 781 NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, & 625 & ri_c, nitke, nmxl, npdl, nave, navb 782 ri_c, nitke, nmxl, npdl, nave, navb, & 783 ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr, & 784 ln_lc , rn_lc 626 785 !!---------------------------------------------------------------------- 627 786 … … 659 818 WRITE(numout,*) ' and its associated coeff. eboost = ', eboost 660 819 WRITE(numout,*) ' constant background or profile navb = ', navb 820 WRITE(numout,*) ' flag for compu.of bls as fun. of wind ln_lsfc = ', ln_lsfc 821 WRITE(numout,*) ' buoyancy lenght scale surface min value rn_lmin = ', rn_lmin 822 WRITE(numout,*) ' horizontal variation for avtb nn_avtb_2d = ', nn_avtb_2d 823 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 824 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 825 WRITE(numout,*) ' % of emin which pene. the thermocline rn_efr = ', rn_efr 826 WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc 827 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc 661 828 WRITE(numout,*) 662 829 ENDIF … … 665 832 IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( ' bad flag: nmxl is < 0 or > 3 ' ) 666 833 IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( ' bad flag: npdl is < 0 or > 1 ' ) 834 IF( nn_htau < 0 .OR. nn_htau > 2 ) CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' ) 835 IF( rn_lc < 0.15 .OR. rn_lc > 0.2 ) CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' ) 836 837 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln 667 838 668 839 ! Horizontal average : initialization of weighting arrays … … 743 914 ! Background eddy viscosity and diffusivity profil 744 915 ! ------------------------------------------------ 745 IF( navb == 0 ) THEN 746 ! Define avmb, avtb from namelist parameter 916 IF( navb == 0 ) THEN ! Define avmb, avtb from namelist parameter 747 917 avmb(:) = avm0 748 918 avtb(:) = avt0 749 ELSE 750 ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 919 ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 751 920 avmb(:) = avm0 752 921 !!bug this is not valide neither in scoord 753 922 IF(ln_sco .AND. lwp) WRITE(numout,cform_war) 754 IF(ln_sco .AND. lwp) WRITE(numout,*) ' avtb profile no rt valid in sco'923 IF(ln_sco .AND. lwp) WRITE(numout,*) ' avtb profile not valid in sco' 755 924 756 925 avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:) ! m2/s 757 926 ENDIF 758 927 928 ! ! 2D shape of the avtb 929 avtb_2d(:,:) = 1.e0 ! uniform 930 ! 931 IF( nn_avtb_2d == 1 ) THEN ! decrease avtb in the equatorial band 932 ! -15S -5S : linear decrease from avt0 to avt0/10. 933 ! -5S +5N : cst value avt0/10. 934 ! 5N 15N : linear increase from avt0/10, to avt0 935 WHERE(-15. <= gphit .AND. gphit < -5 ) avtb_2d = (1. - 0.09 * (gphit + 15.)) 936 WHERE( -5. <= gphit .AND. gphit < 5 ) avtb_2d = 0.1 937 WHERE( 5. <= gphit .AND. gphit < 15 ) avtb_2d = (0.1 + 0.09 * (gphit - 5.)) 938 ENDIF 939 759 940 ! Increase the background in the surface layers 941 !!gm the increase is only required when using cen2 advection scheme 942 !!gm for the equatorial upwelling. 760 943 avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1) 761 944 avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2)
Note: See TracChangeset
for help on using the changeset viewer.