Changeset 1050


Ignore:
Timestamp:
2008-06-03T12:25:26+02:00 (12 years ago)
Author:
ctlod
Message:

trunk: update TKE physics, see ticket: #183

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/ZDF/zdftke.F90

    r899 r1050  
    1818   !!             9.0  !  02-08  (G. Madec)  autotasking optimization 
    1919   !!             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  
    2025   !!---------------------------------------------------------------------- 
    2126#if defined key_zdftke   ||   defined key_esopa 
     
    3338   USE sbc_oce         ! surface boundary condition: ocean 
    3439   USE phycst          ! physical constants 
     40   USE zdfmxl          ! mixed layer 
    3541   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3642   USE prtctl          ! Print control 
     
    6773      &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2) 
    6874      &                  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 
    6986 
    7087# if defined key_c1d 
     
    148165         &          zmxlm => ta,  &  ! use ta as workspace 
    149166         &          zmxld => sa      ! use sa as workspace 
     167      USE oce     , ztkelc => va     ! use va as workspace 
    150168      ! 
    151169      INTEGER, INTENT(in) ::   kt ! ocean time step 
     
    160178         &          zsqen, zesh2,            &  ! 
    161179         &          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 
    162185      !!-------------------------------------------------------------------- 
    163186 
     
    165188 
    166189      !                                            ! Local constant initialization 
    167       zmlmin = 1.e-8 
     190      zmlmin = 0.4 
    168191      zbbrau =  .5 * ebb / rau0 
    169192      zfact1 = -.5 * rdt * efave 
     
    177200      ! Buoyancy length scale: l=sqrt(2*e/n**2) 
    178201      ! --------------------- 
    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 
    181219!CDIR NOVERRCHK 
    182220      DO jk = 2, jpkm1 
     
    193231      ! Physical limits for the mixing length 
    194232      ! ------------------------------------- 
    195       zmxld(:,:, 1 ) = zmlmin   ! surface set to the minimum value 
    196       zmxld(:,:,jpk) = zmlmin   ! bottom  set to the minimum value 
     233      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
     234      zmxld(:,:,jpk) = rn_lmin           ! bottom  set to the minimum value 
    197235 
    198236      SELECT CASE ( nmxl ) 
     
    277315 
    278316      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    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 
    280373      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    281374 
     
    290383               zsqen = SQRT( en(ji,jj,jk) ) 
    291384               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) 
    293386               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 
    294387               zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
     
    309402            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
    310403            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) 
    312407            zmxlm(ji,jj,jpk) = 0.e0 
    313408         END DO 
     
    347442                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 
    348443                  ! 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 
    351452               END DO 
    352453            END DO 
     
    390491                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 
    391492                  ! 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 
    395502                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk) 
    396503               END DO 
     
    454561         END DO 
    455562      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 
    456613 
    457614      ! Lateral boundary conditions on ( avt, en )   (sign unchanged) 
     
    459616 
    460617      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    461       ! III.  Before vertical eddy vicosity and diffusivity coefficients 
     618      ! IV.  Before vertical eddy vicosity and diffusivity coefficients 
    462619      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    463620 
     
    563720               DO ji = fs_2, fs_jpim1   ! vector opt. 
    564721                  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) 
    566723               END DO 
    567724            END DO 
     
    621778      INTEGER ::           jk   ! dummy loop indices 
    622779# endif 
    623  
     780      !! 
    624781      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  
    626785      !!---------------------------------------------------------------------- 
    627786 
     
    659818         WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost 
    660819         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 
    661828         WRITE(numout,*) 
    662829      ENDIF 
     
    665832      IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' ) 
    666833      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  
    667838 
    668839      ! Horizontal average : initialization of weighting arrays  
     
    743914      ! Background eddy viscosity and diffusivity profil 
    744915      ! ------------------------------------------------ 
    745       IF( navb == 0 ) THEN 
    746          !   Define avmb, avtb from namelist parameter 
     916      IF( navb == 0 ) THEN      ! Define avmb, avtb from namelist parameter 
    747917         avmb(:) = avm0 
    748918         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)  
    751920         avmb(:) = avm0 
    752921!!bug  this is not valide neither in scoord 
    753922         IF(ln_sco .AND. lwp) WRITE(numout,cform_war) 
    754          IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile nort valid in sco' 
     923         IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile not valid in sco' 
    755924 
    756925         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s 
    757926      ENDIF 
    758927 
     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 
    759940      !   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. 
    760943      avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1) 
    761944      avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2) 
Note: See TracChangeset for help on using the changeset viewer.