New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 1705 for trunk/NEMO/OPA_SRC/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2009-11-03T17:37:00+01:00 (14 years ago)
Author:
smasson
Message:

impact of HF winds in TKE, see ticket:585

File:
1 edited

Legend:

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

    r1695 r1705  
    6565 
    6666   !                                       !!! ** Namelist  namzdf_tke  ** 
    67    LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not 
    68    INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3) 
    69    REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length   [m] 
    70    REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length   [m] 
    71    INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1) 
    72    REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
    73    REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation  
    74    REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke 
    75    REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke           [m2/s2] 
    76    REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke   [m2/s2] 
    77    REAL(wp) ::   rn_bshear= 1.e-20          ! background shear (>0) 
    78    INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2) 
    79    INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1) 
    80    REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean 
    81    LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not 
    82    REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells 
     67   LOGICAL  ::   ln_mxl0   = .FALSE.         ! mixing length scale surface value as function of wind stress or not 
     68   INTEGER  ::   nn_mxl    =  2              ! type of mixing length (=0/1/2/3) 
     69   REAL(wp) ::   rn_lmin0  = 0.4_wp          ! surface  min value of mixing length   [m] 
     70   REAL(wp) ::   rn_lmin   = 0.1_wp          ! interior min value of mixing length   [m] 
     71   INTEGER  ::   nn_pdl    =  1              ! Prandtl number or not (ratio avt/avm) (=0/1) 
     72   REAL(wp) ::   rn_ediff  = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     73   REAL(wp) ::   rn_ediss  = 0.7_wp          ! coefficient of the Kolmogoroff dissipation  
     74   REAL(wp) ::   rn_ebb    = 3.75_wp         ! coefficient of the surface input of tke 
     75   REAL(wp) ::   rn_emin   = 0.7071e-6_wp    ! minimum value of tke           [m2/s2] 
     76   REAL(wp) ::   rn_emin0  = 1.e-4_wp        ! surface minimum value of tke   [m2/s2] 
     77   REAL(wp) ::   rn_bshear = 1.e-20          ! background shear (>0) 
     78   INTEGER  ::   nn_etau   = 0               ! type of depth penetration of surface tke (=0/1/2) 
     79   INTEGER  ::   nn_htau   = 0               ! type of tke profile of penetration (=0/1) 
     80   REAL(wp) ::   rn_efr    = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean 
     81   REAL(wp) ::   rn_addhft = 0.0_wp          ! add offset   applied to HF tau 
     82   REAL(wp) ::   rn_sclhft = 1.0_wp          ! scale factor applied to HF tau 
     83   LOGICAL  ::   ln_lc     = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not 
     84   REAL(wp) ::   rn_lc     = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells 
    8385 
    8486   REAL(wp) ::   ri_cri   ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    177179      USE oce,   zd_lw  =>   ta   ! use ta as workspace 
    178180      !! 
    179       INTEGER  ::   ji, jj, jk                ! dummy loop arguments 
    180       REAL(wp) ::   zbbrau, zesh2             ! temporary scalars 
    181       REAL(wp) ::   zfact1, zfact2, zfact3    !    -         - 
    182       REAL(wp) ::   ztx2  , zty2  , zcof      !    -         - 
    183       REAL(wp) ::   zus   , zwlc  , zind      !    -         - 
    184       REAL(wp) ::   zzd_up, zzd_lw            !    -         - 
     181      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
    185182      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar 
    186183      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar 
    187       REAL(wp) ::   zebot                           ! temporary scalars 
     184      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
     185      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
     186      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
     187      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
     188      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
     189      REAL(wp) ::   ztau  , zdif                    !    -         - 
     190      REAL(wp) ::   zus   , zwlc  , zind            !    -         - 
     191      REAL(wp) ::   zzd_up, zzd_lw                  !    -         - 
     192      REAL(wp) ::   zebot                           !    -         - 
    188193      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace 
    189194      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      - 
     
    242247         END DO 
    243248         !                        !* finite Langmuir Circulation depth 
     249         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    244250         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy 
    245251         DO jk = jpkm1, 2, -1 
    246252            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us  
    247253               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1) 
    248                   zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj) 
     254                  zus  = zcof * taum(ji,jj) 
    249255                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
    250256               END DO 
     
    265271         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth 
    266272# endif 
     273         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    267274!CDIR NOVERRCHK 
    268275         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
     
    271278!CDIR NOVERRCHK 
    272279               DO ji = fs_2, fs_jpim1   ! vector opt. 
    273 !!gm replace here wndn by a formulation with the stress module 
    274                   zus  = 0.016 * wndm(ji,jj)                  ! Stokes drift 
     280                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    275281                  !                                           ! vertical velocity due to LC 
    276282                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 
     
    371377      !                            !  TKE due to surface and internal wave breaking 
    372378      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    373       IF( nn_etau == 1 ) THEN           !*  penetration throughout the water column 
     379      IF( nn_etau == 1 ) THEN           !* penetration throughout the water column 
    374380         DO jk = 2, jpkm1 
    375381            DO jj = 2, jpjm1 
     
    386392               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    387393                  &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     394            END DO 
     395         END DO 
     396      ELSEIF( nn_etau == 3 ) THEN       !* penetration throughout the water column 
     397!CDIR NOVERRCHK 
     398         DO jk = 2, jpkm1 
     399!CDIR NOVERRCHK 
     400            DO jj = 2, jpjm1 
     401!CDIR NOVERRCHK 
     402               DO ji = fs_2, fs_jpim1   ! vector opt. 
     403                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
     404                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
     405                  ztau = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 )      ! module of the mean stress  
     406                  zdif = taum(ji,jj) - ztau                           ! mean of the module - module of the mean  
     407                  zdif = rn_sclhft * MAX( 0.e0, zdif + rn_addhft )    ! apply some modifications... 
     408                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     409                     &                                        * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 
     410               END DO 
    388411            END DO 
    389412         END DO 
     
    642665      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    643666      !! 
    644       NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb, rn_emin,   & 
    645          &                 rn_emin0, rn_bshear, nn_mxl, ln_mxl0,   & 
    646          &                 rn_lmin , rn_lmin0 , nn_pdl, nn_etau,   & 
    647          &                 nn_htau , rn_efr   , ln_lc , rn_lc  
     667      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,   & 
     668         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,   & 
     669         &                 rn_lmin , rn_lmin0 , nn_pdl   , nn_etau  ,   & 
     670         &                 nn_htau , rn_efr   , rn_addhft, rn_sclhft,   & 
     671         &                 ln_lc   , rn_lc  
    648672      !!---------------------------------------------------------------------- 
    649673 
     
    658682         WRITE(numout,*) '~~~~~~~~' 
    659683         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters' 
    660          WRITE(numout,*) '      coef. to compute avt                        rn_ediff = ', rn_ediff 
    661          WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss = ', rn_ediss 
    662          WRITE(numout,*) '      tke surface input coef.                     rn_ebb   = ', rn_ebb 
    663          WRITE(numout,*) '      minimum value of tke                        rn_emin  = ', rn_emin 
    664          WRITE(numout,*) '      surface minimum value of tke                rn_emin0 = ', rn_emin0 
    665          WRITE(numout,*) '      background shear (>0)                       rn_bshear= ', rn_bshear 
    666          WRITE(numout,*) '      mixing length type                          nn_mxl   = ', nn_mxl 
    667          WRITE(numout,*) '      prandl number flag                          nn_pdl   = ', nn_pdl 
    668          WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0  = ', ln_mxl0 
    669          WRITE(numout,*) '      surface  mixing length minimum value        rn_lmin0 = ', rn_lmin0 
    670          WRITE(numout,*) '      interior mixing length minimum value        rn_lmin0 = ', rn_lmin0 
    671          WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau  = ', nn_etau 
    672          WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau  = ', nn_htau 
    673          WRITE(numout,*) '      % of rn_emin0 which pene. the thermocline   rn_efr   = ', rn_efr 
    674          WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc    = ', ln_lc 
    675          WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc    = ', rn_lc 
     684         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff 
     685         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss 
     686         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb 
     687         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin 
     688         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0 
     689         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear 
     690         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl 
     691         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl 
     692         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0 
     693         WRITE(numout,*) '      surface  mixing length minimum value        rn_lmin0  = ', rn_lmin0 
     694         WRITE(numout,*) '      interior mixing length minimum value        rn_lmin0  = ', rn_lmin0 
     695         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau 
     696         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau 
     697         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr 
     698         WRITE(numout,*) '      add offset   applied to HF tau              rn_addhft = ', rn_addhft 
     699         WRITE(numout,*) '      scale factor applied to HF tau              rn_sclhft = ', rn_sclhft 
     700         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc 
     701         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc 
    676702         WRITE(numout,*) 
    677703         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri 
     
    679705 
    680706      !                               !* Check of some namelist values 
    681       IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    682       IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    683       IF( nn_htau < 0    .OR. nn_htau > 1   )   CALL ctl_stop( 'bad flag: nn_htau is 0 or 1    ' ) 
    684       IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' ) 
     707      IF( nn_mxl  < 0     .OR.  nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
     708      IF( nn_pdl  < 0     .OR.  nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
     709      IF( nn_htau < 0     .OR.  nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     710      IF( rn_lc   < 0.15  .OR.  rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' ) 
     711      IF( rn_sclhft == 0. .AND. nn_etau == 3  )   CALL ctl_stop( 'force null HF tau to penetrate the thermocline...' ) 
     712      IF( .NOT. lhftau    .AND. nn_etau == 3  )   CALL ctl_stop( 'bad flag: nn_etau == 3 must be used with HF tau' ) 
    685713 
    686714      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     
    697725               END DO 
    698726            END DO 
     727         CASE( 2 )                                    ! constant depth penetration (here 30 meters) 
     728            htau(:,:) = 30.e0 
    699729         END SELECT 
    700730      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.