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 13403 – NEMO

Changeset 13403


Ignore:
Timestamp:
2020-08-14T17:08:06+02:00 (4 years ago)
Author:
dancopsey
Message:

Merge in George Nurser's improvements to Stokes Drift in OSMOSIS.
Custom merge from revision 13392 to revision 13400 of:
http://forge.ipsl.jussieu.fr/nemo/log/NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0

Location:
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/field_def_nemo-oce.xml

    r13402 r13403  
    789789 
    790790     <!-- OMIP  layer-integrated trends --> 
    791      <field id="sfd"      long_name="Total downward salt flux"                                                     unit="kg/(m^2 s)" > saltflux * 1026.0 * 0.001 </field> 
     791     <field id="sfd"      long_name="Total downward salt flux"                                                     unit="kg/(m^2 s)" > saltflx * 1026.0 * 0.001 </field> 
    792792     <field id="wfo"      long_name="Total downward FW mass flux"                                                   unit="kg/(m^2 s)" > -empmr </field> 
    793793 
  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/namelist_ref

    r13402 r13403  
    55!! namelists    2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 
    66!!                                    namsbc_sas, namtra_qsr, namsbc_rnf, 
    7 !!                                    namsbc_isf, namsbc_iscpl, namsbc_apr,  
     7!!                                    namsbc_isf, namsbc_iscpl, namsbc_apr, 
    88!!                                    namsbc_ssr, namsbc_wave, namberg) 
    99!!              3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) 
     
    8888      cn_domcfg = "domain_cfg"  ! domain configuration filename 
    8989      ! 
    90       ln_closea    = .false.    !  T => keep closed seas (defined by closea_mask field) in the   
     90      ln_closea    = .false.    !  T => keep closed seas (defined by closea_mask field) in the 
    9191      !                         !       domain and apply special treatment of freshwater fluxes. 
    92       !                         !  F => suppress closed seas (defined by closea_mask field)  
     92      !                         !  F => suppress closed seas (defined by closea_mask field) 
    9393      !                         !       from the bathymetry at runtime. 
    9494      !                         !  If closea_mask field doesn't exist in the domain_cfg file 
     
    106106   ln_tsd_init = .false.         !  ocean initialisation 
    107107   ln_tsd_dmp  = .false.         !  T-S restoring   (see namtra_dmp) 
    108     
     108 
    109109   cn_dir      = './'      !  root directory for the T-S data location 
    110110   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
     
    195195   nn_fsbc     = 2         !  frequency of SBC module call 
    196196      !                    !  (control sea-ice & iceberg model call) 
    197                      ! Type of air-sea fluxes  
     197                     ! Type of air-sea fluxes 
    198198   ln_usr      = .false.   !  user defined formulation                  (T => check usrdef_sbc) 
    199199   ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx ) 
     
    205205      !                    !  =0 no opa-sas OASIS coupling: default single executable config. 
    206206      !                    !  =1 opa-sas OASIS coupling: multi executable config., OPA component 
    207       !                    !  =2 opa-sas OASIS coupling: multi executable config., SAS component  
     207      !                    !  =2 opa-sas OASIS coupling: multi executable config., SAS component 
    208208                     ! Sea-ice : 
    209    nn_ice      = 0         !  =0 no ice boundary condition     
     209   nn_ice      = 0         !  =0 no ice boundary condition 
    210210      !                    !  =1 use observed ice-cover                 (  => fill namsbc_iif ) 
    211211      !                    !  =2 or 3 automatically for SI3 or CICE    ("key_si3" or "key_cice") 
     
    213213   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges) 
    214214      !                    !  =F levitating ice (no pressure, mass and salt exchanges) 
    215                      ! Misc. options of sbc :  
     215                     ! Misc. options of sbc : 
    216216   ln_traqsr   = .false.   !  Light penetration in the ocean            (T => fill namtra_qsr) 
    217217   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
     
    225225   ln_wave     = .false.   !  Activate coupling with wave  (T => fill namsbc_wave) 
    226226   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave) 
    227    ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)  
     227   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 
    228228   nn_sdrift   =  0        !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 
    229229      !                    !   = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
     
    379379   nn_chldta   =      0       !  RGB : Chl data (=1) or cst value (=0) 
    380380   rn_si1      =   23.0       !  2BD : longest depth of extinction 
    381     
     381 
    382382   cn_dir      = './'      !  root directory for the chlorophyl data location 
    383383   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
     
    443443/ 
    444444!----------------------------------------------------------------------- 
    445 &namsbc_isf    !  Top boundary layer (ISF)                              (ln_isfcav =T : read (ln_read_cfg=T)  
     445&namsbc_isf    !  Top boundary layer (ISF)                              (ln_isfcav =T : read (ln_read_cfg=T) 
    446446!-----------------------------------------------------------------------             or set or usr_def_zgr ) 
    447    !                 ! type of top boundary layer  
     447   !                 ! type of top boundary layer 
    448448   nn_isf      = 1         !  ice shelf melting/freezing 
    449                            !  1 = presence of ISF   ;  2 = bg03 parametrisation  
     449                           !  1 = presence of ISF   ;  2 = bg03 parametrisation 
    450450                           !  3 = rnf file for ISF  ;  4 = ISF specified freshwater flux 
    451451                           !  options 1 and 4 need ln_isfcav = .true. (domzgr) 
     
    470470!* nn_isf = 3 case 
    471471   sn_rnfisf   = 'rnfisf'    ,         -12.      ,'sofwfisf' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    472 !* nn_isf = 2 and 3 cases  
     472!* nn_isf = 2 and 3 cases 
    473473   sn_depmax_isf ='rnfisf'   ,         -12.      ,'sozisfmax',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    474474   sn_depmin_isf ='rnfisf'   ,         -12.      ,'sozisfmin',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     
    477477/ 
    478478!----------------------------------------------------------------------- 
    479 &namsbc_iscpl  !   land ice / ocean coupling option                     (ln_isfcav =T : read (ln_read_cfg=T)  
     479&namsbc_iscpl  !   land ice / ocean coupling option                     (ln_isfcav =T : read (ln_read_cfg=T) 
    480480!-----------------------------------------------------------------------             or set or usr_def_zgr ) 
    481481   nn_drown    = 10        ! number of iteration of the extrapolation loop (fill the new wet cells) 
     
    577577         ln_read_load  = .false.               ! Or read load potential from file 
    578578            cn_tide_load = 'tide_LOAD_grid_T.nc'  ! filename for load potential 
    579             !       
     579            ! 
    580580      ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup 
    581581         rdttideramp   =    0.                 !  ramp duration in days 
     
    656656   filtide          = 'bdydta/amm12_bdytide_'   !  file name root of tidal forcing files 
    657657   ln_bdytide_2ddta = .false.                   ! 
    658    ln_bdytide_conj  = .false.                   !  
     658   ln_bdytide_conj  = .false.                   ! 
    659659/ 
    660660 
     
    683683!----------------------------------------------------------------------- 
    684684   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
    685    rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     685   rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0) 
    686686   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
    687687   rn_ke0      =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     
    694694!----------------------------------------------------------------------- 
    695695   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
    696    rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     696   rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0) 
    697697   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
    698698   rn_ke0      =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     
    761761      nn_cen_v   =  4            !  =2/4, vertical   2nd order CEN / 4th order COMPACT 
    762762   ln_traadv_fct = .false. !  FCT scheme 
    763       nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order  
    764       nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order  
     763      nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order 
     764      nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order 
    765765   ln_traadv_mus = .false. !  MUSCL scheme 
    766766      ln_mus_ups = .false.       !  use upstream scheme near river mouths 
     
    783783   ln_traldf_triad = .false.   !  iso-neutral (triad    operator) 
    784784   ! 
    785    !                       !  iso-neutral options:         
     785   !                       !  iso-neutral options: 
    786786   ln_traldf_msc   = .false.   !  Method of Stabilizing Correction      (both operators) 
    787787   rn_slpmax       =  0.01     !  slope limit                           (both operators) 
     
    793793   nn_aht_ijk_t    = 0         !  space/time variation of eddy coefficient: 
    794794      !                             !   =-20 (=-30)    read in eddy_diffusivity_2D.nc (..._3D.nc) file 
    795       !                             !   =  0           constant  
    796       !                             !   = 10 F(k)      =ldf_c1d  
    797       !                             !   = 20 F(i,j)    =ldf_c2d  
     795      !                             !   =  0           constant 
     796      !                             !   = 10 F(k)      =ldf_c1d 
     797      !                             !   = 20 F(i,j)    =ldf_c2d 
    798798      !                             !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation 
    799799      !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d 
    800800      !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
    801       !                        !  time invariant coefficients:  aht0 = 1/2  Ud*Ld   (lap case)  
     801      !                        !  time invariant coefficients:  aht0 = 1/2  Ud*Ld   (lap case) 
    802802      !                             !                           or   = 1/12 Ud*Ld^3 (blp case) 
    803803      rn_Ud        = 0.01           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) 
     
    825825      nn_aei_ijk_t    = 0           !  space/time variation of eddy coefficient: 
    826826      !                             !   =-20 (=-30)    read in eddy_induced_velocity_2D.nc (..._3D.nc) file 
    827       !                             !   =  0           constant  
    828       !                             !   = 10 F(k)      =ldf_c1d  
    829       !                             !   = 20 F(i,j)    =ldf_c2d  
     827      !                             !   =  0           constant 
     828      !                             !   = 10 F(k)      =ldf_c1d 
     829      !                             !   = 20 F(i,j)    =ldf_c2d 
    830830      !                             !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation 
    831831      !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d 
    832       !                        !  time invariant coefficients:  aei0 = 1/2  Ue*Le  
     832      !                        !  time invariant coefficients:  aei0 = 1/2  Ue*Le 
    833833      rn_Ue        = 0.02           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) 
    834834      rn_Le        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10) 
     
    890890   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) 
    891891   ln_dynvor_een = .false. !  energy & enstrophy scheme 
    892       nn_een_e3f = 0          ! =0  e3f = mi(mj(e3t))/4  
     892      nn_een_e3f = 0          ! =0  e3f = mi(mj(e3t))/4 
    893893      !                       ! =1  e3f = mi(mj(e3t))/mi(mj( tmask)) 
    894894   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T)        ==>>> PLEASE DO NOT ACTIVATE 
     
    935935      !                             !  =-30  read in eddy_viscosity_3D.nc file 
    936936      !                             !  =-20  read in eddy_viscosity_2D.nc file 
    937       !                             !  =  0  constant  
     937      !                             !  =  0  constant 
    938938      !                             !  = 10  F(k)=c1d 
    939939      !                             !  = 20  F(i,j)=F(grid spacing)=c2d 
     
    941941      !                             !  = 31  F(i,j,k)=F(grid spacing and local velocity) 
    942942      !                             !  = 32  F(i,j,k)=F(local gridscale and deformation rate) 
    943       !                        !  time invariant coefficients :  ahm = 1/2  Uv*Lv   (lap case)  
     943      !                        !  time invariant coefficients :  ahm = 1/2  Uv*Lv   (lap case) 
    944944      !                             !                            or  = 1/12 Uv*Lv^3 (blp case) 
    945945      rn_Uv      = 0.1              !  lateral viscous velocity [m/s] (nn_ahm_ijk_t= 0, 10, 20, 30) 
     
    10651065                              !        = 0  constant 10 m length scale 
    10661066                              !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
    1067       rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4    
     1067      rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 
    10681068/ 
    10691069!----------------------------------------------------------------------- 
     
    10911091   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    10921092   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    1093    rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     1093   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor 
     1094   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4 
    10941095   nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv 
    10951096   ln_dia_osm = .true.         ! output OSMOSIS-OBL variables 
     
    10991100   rn_difri  =  0.005          ! max Ri# diffusivity at Ri_g = 0 (m^2/s) 
    11001101   ln_convmix  = .true.        ! Use convective instability mixing below BL 
    1101    rn_difconv = 1.             ! diffusivity when unstable below BL  (m2/s) 
     1102   rn_difconv = 1. !0.01 !1.             ! diffusivity when unstable below BL  (m2/s) 
     1103   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
    11021104   nn_osm_wave = 0             ! Method used to calculate Stokes drift 
    11031105      !                        !  = 2: Use ECMWF wave fields 
    11041106      !                        !  = 1: Pierson Moskowitz wave spectrum 
    11051107      !                        !  = 0: Constant La# = 0.3 
     1108   nn_osm_SD_reduce = 0        ! Method used to get active Stokes drift from surface value 
     1109      !                        !  = 0: No reduction 
     1110                               !  = 1: use SD avged over top 10% hbl 
     1111                               !  = 2:use surface value of SD fit to slope at rn_osm_hblfrac*hbl below surface 
     1112   ln_zdfosm_ice_shelter = .true.  ! reduce surface SD and depth scale under ice 
    11061113   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
    11071114/ 
     
    11161123   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
    11171124   rn_osm_mle_thresh = 0.0005      ! delta b criterion used for FK MLD criterion 
    1118    rn_osm_mle_tau     = 172800.   ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
    1119 / 
     1125   rn_osm_mle_tau     = 172800.    ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
     1126   ln_osm_hmle_limit   = .false.   ! limit hmle to rn_osm_hmle_limit*hbl 
     1127   rn_osm_hmle_limit   = 1.2 
     1128   / 
    11201129!----------------------------------------------------------------------- 
    11211130&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/ref_cfgs.txt

    r13402 r13403  
    99ORCA2_ICE_PISCES OCE TOP ICE NST 
    1010SPITZ12 OCE ICE 
    11 eORCA1 OCE ICE 
  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90

    r13402 r13403  
    103103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    104104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
     105   REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by 
     106   REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl 
     107   LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering 
    105108   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs 
    106109   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt 
    107110   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave 
     111   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value 
    108112   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la 
    109  
    110113 
    111114   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing 
     
    123126   !                                        ! parameters used in nn_osm_mle = 1 case 
    124127   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
     128   LOGICAL  ::   ln_osm_hmle_limit           ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
     129   REAL(wp) ::   rn_osm_hmle_limit           ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
    125130   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
    126131   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
     
    321326      REAL(wp) :: zzeta_v = 0.46 
    322327      REAL(wp) :: zabsstke 
     328      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 
     329      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 
    323330 
    324331      ! For debugging 
     
    407414           zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    408415           ! Surface upward velocity fluxes 
    409            zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
    410            zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
     416           zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rau0 * tmask(ji,jj,1) 
     417           zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rau0 * tmask(ji,jj,1) 
    411418           ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    412419           zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
     
    423430              zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
    424431              zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
     432              ! Linearly 
    425433              zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 
    426               ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 
     434              dstokes(ji,jj) = rn_osm_dstokes 
    427435           END DO 
    428436        END DO 
     
    432440           DO ji = 2, jpim1 
    433441              ! Use wind speed wndm included in sbc_oce module 
    434               zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    435               dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
     442              zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     443              dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    436444           END DO 
    437445        END DO 
     
    439447     CASE(2) 
    440448        zfac =  2.0_wp * rpi / 16.0_wp 
     449 
    441450        DO jj = 2, jpjm1 
    442451           DO ji = 2, jpim1 
    443               ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 
    444               !    The coefficient 0.8 gives La=0.3  in this situation. 
    445               ! It could represent the effects of the spread of wave directions 
    446               ! around the mean wind. The effect of this adjustment needs to be tested. 
    447               zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
    448               zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
    449               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
     452              IF (hsw(ji,jj) > 1.e-4) THEN 
     453                 ! Use  wave fields 
     454                 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     455                 zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     456                 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) 
     457              ELSE 
     458                 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 
     459                 ! .. so default to Pierson-Moskowitz 
     460                 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     461                 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
     462              END IF 
     463            END DO 
     464        END DO 
     465     END SELECT 
     466 
     467     IF (ln_zdfosm_ice_shelter) THEN 
     468        ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 
     469        DO jj = 2, jpjm1 
     470           DO ji = 2, jpim1 
     471              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     472              dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     473           END DO 
     474        END DO 
     475     END IF 
     476 
     477     SELECT CASE (nn_osm_SD_reduce) 
     478     ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van  Roekel (2012) or Grant (2020). 
     479     CASE(0) 
     480        ! The Langmur number from the ECMWF model (or from PM)  appears to give La<0.3 for wind-driven seas. 
     481        !    The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3  in this situation. 
     482        ! It could represent the effects of the spread of wave directions 
     483        ! around the mean wind. The effect of this adjustment needs to be tested. 
     484        IF(nn_osm_wave > 0) THEN 
     485           zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) 
     486        END IF 
     487     CASE(1) 
     488        ! van  Roekel (2012): consider average SD over top 10% of boundary layer 
     489        ! assumes approximate depth profile of SD from Breivik (2016) 
     490        zsqrtpi = SQRT(rpi) 
     491        z_two_thirds = 2.0_wp / 3.0_wp 
     492 
     493        DO jj = 2, jpjm1 
     494           DO ji = 2, jpim1 
     495              zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     496              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     497              zsqrt_depth = SQRT(z2k_times_thickness) 
     498              zexp_depth  = EXP(-z2k_times_thickness) 
     499              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth  & 
     500                   &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & 
     501                   &              + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness 
     502 
     503           END DO 
     504        END DO 
     505     CASE(2) 
     506        ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer 
     507        ! assumes approximate depth profile of SD from Breivik (2016) 
     508        zsqrtpi = SQRT(rpi) 
     509 
     510        DO jj = 2, jpjm1 
     511           DO ji = 2, jpim1 
     512              zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     513              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     514 
     515              IF(z2k_times_thickness < 50._wp) THEN 
     516                 zsqrt_depth = SQRT(z2k_times_thickness) 
     517                 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) 
     518              ELSE 
     519                 ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness 
     520                 ! See Abramowitz and Stegun, Eq. 7.1.23 
     521                 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness)  + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 
     522                 zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp 
     523              END IF 
     524              zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 
     525              dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 
     526              zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 
    450527           END DO 
    451528        END DO 
     
    486563     ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 
    487564     ! previously exist for hbl also. 
     565 
     566     ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 
     567     ! ########################################################################## 
    488568      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
    489569      ibld(:,:) = 4 
     
    497577         END DO 
    498578      END DO 
     579     ! ########################################################################## 
    499580 
    500581      DO jj = 2, jpjm1 
     
    10701151 
    10711152 
    1072   
     1153 
    10731154       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
    10741155          DO jj = 2 , jpjm1 
     
    11631244            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    11641245         ! Stokes drift read in from sbcwave  (=2). 
    1165          CASE(2) 
     1246         CASE(2:3) 
    11661247            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
    11671248            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     
    15181599             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
    15191600             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
    1520              zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
    1521                   &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1601             IF (nn_osm_SD_reduce > 0 ) THEN 
     1602                ! Effective Stokes drift already reduced from surface value 
     1603                zr_stokes = 1.0_wp 
     1604             ELSE 
     1605                ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 
     1606                ! requires further reduction where BL is deep 
     1607                zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1608                     &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1609             END IF 
    15221610 
    15231611             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     
    15561644 
    15571645                  zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    1558                   &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1646                  &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 
    15591647                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    15601648                ! added ajgn 23 July as temporay fix 
     
    15891677                    zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    15901678                ENDIF 
    1591                 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 
     1679                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
    15921680            ENDIF 
    15931681         END DO 
     
    17221810                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
    17231811                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
    1724                         IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1812                        IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
    17251813                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    1726                                 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1814                                (1.0 + zdb_bl(ji,jj)**2 / MAX( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
    17271815                        ELSE                                                     ! unstable 
    17281816                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    17291817                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
    17301818                        ENDIF 
    1731                         ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1819                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    17321820                        zddhdt = zari * hbl(ji,jj) 
    17331821                     ELSE 
    1734                         ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1822                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    17351823                        zddhdt = 0.2 * hbl(ji,jj) 
    17361824                     ENDIF 
    17371825                  ELSE 
    1738                      ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1826                     ztau = 0.2 * hbl(ji,jj) /  MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    17391827                     zddhdt = 0.2 * hbl(ji,jj) 
    17401828                  ENDIF 
    17411829               ELSE ! ln_osm_mle 
    17421830                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
    1743                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    1744                         zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    1745                              (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1831                     IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1832                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 
     1833                             (1.0 + zdb_bl(ji,jj)**2 /  MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
    17461834                     ELSE                                                     ! unstable 
    1747                         zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    1748                              (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1835                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 
     1836                             (1.0 + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
    17491837                     ENDIF 
    1750                      ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1838                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    17511839                     zddhdt = zari * hbl(ji,jj) 
    17521840                  ELSE 
    1753                      ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1841                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    17541842                     zddhdt = 0.2 * hbl(ji,jj) 
    17551843                  ENDIF 
     
    17601848               ! Alan: this hml is never defined or used 
    17611849            ELSE   ! IF (lconv) 
    1762                ztau = hbl(ji,jj) / zvstr(ji,jj) 
     1850               ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 
    17631851               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    17641852                  ! boundary layer deepening 
     
    17661854                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    17671855                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    1768                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     1856                          & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 ) 
    17691857                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
    17701858                  ELSE 
     
    18211909      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
    18221910      DO jk = nlb10, jpkm1 
    1823          DO jj = 1, jpj                ! Mixed layer level: w-level  
     1911         DO jj = 1, jpj                ! Mixed layer level: w-level 
    18241912            DO ji = 1, jpi 
    18251913               ikt = mbkt(ji,jj) 
     
    19402028          ENDIF 
    19412029          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
    1942           hmle(ji,jj) = MIN(hmle(ji,jj), 1.2*zmld(ji,jj))  
     2030          IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*zmld(ji,jj)) 
    19432031        END DO 
    19442032      END DO 
     
    19702058                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
    19712059             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 
    1972       ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
    1973              zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1)  
     2060      ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
     2061             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
    19742062             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf 
    19752063          ENDIF 
     
    19982086     !! 
    19992087     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    2000           & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
     2088          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 
    20012089          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
    2002           & ,ln_osm_mle 
     2090          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, ln_zdfosm_ice_shelter 
    20032091! Namelist for Fox-Kemper parametrization. 
    20042092      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
    2005            & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 
     2093           & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit 
    20062094 
    20072095     !!---------------------------------------------------------------------- 
     
    20242112        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    20252113        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
     2114        WRITE(numout,*) '     Stokes drift reduction factor                 rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd 
    20262115        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    20272116        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
     
    20352124        CASE(2) 
    20362125           WRITE(numout,*) '     calculated from ECMWF wave fields' 
     2126         END SELECT 
     2127        WRITE(numout,*) '     Stokes drift reduction                        nn_osm_SD_reduce', nn_osm_SD_reduce 
     2128        WRITE(numout,*) '     fraction of hbl to average SD over/fit' 
     2129        WRITE(numout,*) '     exponential with nn_osm_SD_reduce = 1 or 2    rn_osm_hblfrac =  ', rn_osm_hblfrac 
     2130        SELECT CASE (nn_osm_SD_reduce) 
     2131        CASE(0) 
     2132           WRITE(numout,*) '     No reduction' 
     2133        CASE(1) 
     2134           WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL' 
     2135        CASE(2) 
     2136           WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL' 
    20372137        END SELECT 
     2138        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 
    20382139        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm 
    20392140        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix 
     
    20722173            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
    20732174            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
     2175            WRITE(numout,*) '         switch to limit hmle                                      ln_osm_hmle_limit  = ', ln_osm_hmle_limit 
     2176            WRITE(numout,*) '         fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T.  rn_osm_hmle_limit = ', rn_osm_hmle_limit 
    20742177         ENDIF         ! 
    20752178     ENDIF 
Note: See TracChangeset for help on using the changeset viewer.