Changeset 12912


Ignore:
Timestamp:
2020-05-12T14:51:37+02:00 (7 months ago)
Author:
cguiavarch
Message:

OSMOSIS and Fox-Kemper changes (merged from NEMO_4.0.1_FKOSM)

Location:
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM
Files:
14 edited
2 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/C1D_PAPA/cpp_C1D_PAPA.fcm

    r9799 r12912  
    1  bld::tool::fppkeys   key_c1d key_mpp_mpi key_iomput 
     1 bld::tool::fppkeys   key_c1d key_iomput 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/GYRE_PISCES/EXPREF/namelist_cfg

    r12206 r12912  
    142142      rn_Ld        = 100.e+3        !  lateral diffusive length   [m] 
    143143/ 
     144!----------------------------------------------------------------------- 
     145&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     146!----------------------------------------------------------------------- 
     147/ 
    144148!!====================================================================== 
    145149!!                      ***  Dynamics namelists  ***                  !! 
     
    212216!----------------------------------------------------------------------- 
    213217   nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves 
     218/ 
     219!----------------------------------------------------------------------- 
     220&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
     221!----------------------------------------------------------------------- 
     222/ 
     223!----------------------------------------------------------------------- 
     224&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     225!----------------------------------------------------------------------- 
    214226/ 
    215227!!====================================================================== 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r12206 r12912  
    135135/ 
    136136!----------------------------------------------------------------------- 
     137&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     138!----------------------------------------------------------------------- 
     139/ 
     140!----------------------------------------------------------------------- 
    137141&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T) 
    138142!----------------------------------------------------------------------- 
     
    377381/ 
    378382!----------------------------------------------------------------------- 
     383&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
     384!----------------------------------------------------------------------- 
     385/ 
     386!----------------------------------------------------------------------- 
     387&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     388!----------------------------------------------------------------------- 
     389/ 
     390!----------------------------------------------------------------------- 
    379391&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
    380392!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/SHARED/domain_def_nemo.xml

    r12276 r12912  
    1010     </domain> 
    1111 
    12       <domain id="1point" domain_ref="grid_T" > 
    13          <zoom_domain ibegin="139" jbegin="119" ni="1" nj="1"/> 
    14       </domain> 
    15  
    16  
    17      <!--   Eq section --> 
     12     <domain id="1point" domain_ref="grid_T" > 
     13       <zoom_domain ibegin="1" jbegin="1" ni="1" nj="1"/> 
     14     </domain> 
     15     <domain id="danger_T" domain_ref="grid_T" > 
     16       <zoom_domain id="danger_T"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     17     </domain> 
     18     <domain id="danger_U" domain_ref="grid_U" > 
     19       <zoom_domain id="danger_U"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     20     </domain> 
     21      <domain id="danger_V" domain_ref="grid_V" > 
     22       <zoom_domain id="danger_V"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     23     </domain> 
     24     <domain id="danger_W" domain_ref="grid_W" > 
     25       <zoom_domain id="danger_W"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     26    </domain> 
     27    <!--   Eq section --> 
    1828     <domain id="EqT" domain_ref="grid_T" > <zoom_domain id="EqT"/> </domain> 
    1929     <!--   TAO : see example above   --> 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml

    r12288 r12912  
    215215 
    216216      <field_group id="OSMOSIS_T" grid_ref="grid_T_2D"> 
     217        <field id="hml"                 long_name="mixed layr depth"                         unit="m"       /> 
     218        <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
     219        <field id="dh"                  long_name="Pycnocline thickness"                     unit=" m"      /> 
     220        <field id="ibld"                long_name="index of boundary layer depth"            unit="#"       /> 
     221        <field id="imld"                long_name="index of mixed layer depth"            unit="#"       /> 
     222        <field id="zhbl"                long_name="boundary layer depth -grid"                     unit="m"       /> 
     223        <field id="zhml"                long_name="mixed layer depth - grid"                        unit="m"       /> 
     224        <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
     225        <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
     226        <field id="us_x"        long_name="i component of active Stokes drift"                      unit="m/s"     /> 
     227        <field id="us_y"        long_name="j component of active Stokes drift"                      unit="m/s"     /> 
     228        <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    217229        <field id="zwth0"               long_name="surface non-local temperature flux"       unit="deg m/s" /> 
    218230        <field id="zws0"                long_name="surface non-local salinity flux"          unit="psu m/s" /> 
    219         <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
    220         <field id="hbli"                long_name="initial boundary layer depth"             unit="m"       /> 
    221         <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    222         <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
    223231        <field id="zwstrc"              long_name="convective velocity scale"                unit="m/s"     /> 
     232        <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    224233        <field id="zwstrl"              long_name="langmuir velocity scale"                  unit="m/s"     /> 
    225         <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    226         <field id="zhbl"                long_name="boundary layer depth"                     unit="m"       /> 
    227         <field id="zhml"                long_name="mixed layer depth"                        unit="m"       /> 
     234        <field id="zvstr"               long_name="mixed velocity scale"                     unit="m/s"     /> 
     235        <field id="zla"                 long_name="langmuir number"                          unit="m/s"     /> 
    228236        <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2"                        unit="mW"      /> 
    229237        <field id="wind_wave_power"     long_name="U_s \dot  tau"                            unit="mW"      /> 
    230238        <field id="wind_power"          long_name="\rho  u*^3"                               unit="mW"      /> 
    231239 
    232         <!-- extra OSMOSIS diagnostics --> 
     240       <!-- interior BL OSMOSIS diagnostics --> 
    233241        <field id="zwthav"              long_name="av turb flux of T in ml"                  unit="deg m/s" /> 
    234242        <field id="zt_ml"               long_name="av T in ml"                               unit="deg"     /> 
     243        <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
     244        <field id="zws_ent"            long_name="entrainment turb flux of S"                unit="10^-3 m/s" /> 
    235245        <field id="zwth_ent"            long_name="entrainment turb flux of T"               unit="deg m/s" /> 
    236         <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
    237         <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
    238       </field_group> 
    239  
    240       <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" operation="instant" > 
     246        <field id="zwb_ent"            long_name="entrainment turb flux of buoyancy"         unit="m^2/s^-3" /> 
     247  
     248        <field id="zdt_bl"             long_name="temperature jump at base of BL"                 unit="deg"      /> 
     249        <field id="zds_bl"             long_name="salinity jump at base of BL"                 unit="10^-3"      /> 
     250        <field id="zdb_bl"             long_name="buoyancy jump at base of BL"                 unit="m/s^2"      /> 
     251        <field id="zdu_bl"             long_name="u jump at base of BL"                       unit="m/s"      /> 
     252        <field id="zdv_bl"             long_name="v jump at base of BL"                       unit="m/s"      /> 
     253 
     254        <!-- extra OSMOSIS diagnostics for debugging --> 
     255       <field id="zsc_uw_1_0"       long_name="zsc u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 
     256        <field id="zsc_uw_1_f"       long_name="zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     257        <field id="zsc_vw_1_f"       long_name="zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     258        <field id="zsc_uw_2_f"       long_name="2nd zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     259        <field id="zsc_vw_2_f"       long_name="2nd zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     260        <field id="zuw_bse"       long_name="base u-flux T-points"                          unit="m^2/s^2" /> 
     261        <field id="zvw_bse"       long_name="base v-flux T-points"                          unit="m^2/s^2" /> 
     262 
     263       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     264         <field id="hmle"          long_name="OBL FK-layer thickness"                                     unit="m"        /> 
     265        <field id="mld_prof"              long_name="FK-layer depth index"                  unit="#" /> 
     266        <field id="zmld"          long_name="target FK-layer thickness"                                     unit="m"        /> 
     267        <field id="zwb_fk"          long_name="FK b-flux"                                     unit="m^2 s^-3"        /> 
     268        <field id="zwb_fk_b"          long_name="layer averaged FK b-flux"                 unit="m^2 s^-3"       /> 
     269        <field id="zdiff_mle"          long_name="max FK diffusivity in MLE"       unit=" 10^-4 m^2 s^-1"       /> 
     270        <field id="zvel_mle"          long_name="FK velocity scale in MLE"       unit=" m s^-1"       /> 
     271    </field_group> 
     272 
     273      <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" > 
     274        <field id="zviscos"       long_name="BL viscosity"   unit="m^2/s" /> 
    241275        <field id="ghamt"       long_name="non-local temperature flux"                       unit="deg m/s" /> 
    242276        <field id="ghams"       long_name="non-local salinity flux"                          unit="psu m/s" /> 
    243277        <field id="zdtdz_pyc"   long_name="Pycnocline temperature gradient"                  unit=" deg/m"  /> 
    244       </field_group> 
     278        <field id="zdsdz_pyc"   long_name="Pycnocline salinity gradient"                  unit=" 10^-3/m"  /> 
     279        <field id="zdbdz_pyc"   long_name="Pycnocline buoyancy gradient"                  unit=" s^-2"  /> 
     280        <field id="zdudz_pyc"   long_name="Pycnocline u gradient"                  unit=" s^-2"  /> 
     281        <field id="zdvdz_pyc"   long_name="Pycnocline v gradient"                  unit=" s^-2"  /> 
     282 
     283        <!-- extra OSMOSIS diagnostics for debugging --> 
     284         <field id="ghamu_00"       long_name="initial non-local u-momentum flux"   unit="m^2/s^2" /> 
     285        <field id="ghamv_00"       long_name="initial non-local v-momentum flux"   unit="m^2/s^2" /> 
     286        <field id="ghamu_0"       long_name="after dstokes non-local u-momentum flux"   unit="m^2/s^2" /> 
     287        <field id="ghamu_f"       long_name="after Coriolis non-local u-momentum flux"   unit="m^2/s^2" /> 
     288        <field id="ghamv_f"       long_name="after Coriolis  non-local v-momentum flux"   unit="m^2/s^2" /> 
     289        <field id="ghamu_b"       long_name="after buoyancy added non-local u-momentum flux"   unit="m^2/s^2" /> 
     290        <field id="ghamv_b"       long_name="after buoyancy added  non-local v-momentum flux"  unit="m^2/s^2" /> 
     291        <field id="ghamu_1"       long_name="after entrainment non-local u-momentum flux"   unit="m^2/s^2" /> 
     292        <field id="ghamv_1"       long_name="after entrainment  non-local v-momentum flux"  unit="m^2/s^2" /> 
     293     </field_group> 
    245294 
    246295      <field_group id="OSMOSIS_U" grid_ref="grid_U_2D" > 
    247296        <field id="ghamu"       long_name="non-local u-momentum flux"   grid_ref="grid_U_3D" unit="m^2/s^2" /> 
    248         <field id="us_x"        long_name="i component of Stokes drift"                      unit="m/s"     /> 
    249       </field_group> 
     297       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     298       <field id="zdtdx"          long_name="FK  T x-gradient"                                     unit=" deg C m^-1"        /> 
     299        <field id="zdsdx"          long_name="FK  S x-gradient"                                     unit=" 10^-3 m^-1"        /> 
     300        <field id="dbdx_mle"          long_name="FK  B x-gradient"                                     unit=" s^-2"        /> 
     301     </field_group> 
    250302 
    251303      <field_group id="OSMOSIS_V" grid_ref="grid_V_2D" > 
    252304        <field id="ghamv"       long_name="non-local v-momentum flux"   grid_ref="grid_V_3D" unit="m^2/s^2" /> 
    253         <field id="us_y"        long_name="j component of Stokes drift"                      unit="m/s"     /> 
     305        <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     306        <field id="zdtdy"          long_name="FK T y-gradient"                                     unit=" deg C m^-1"        /> 
     307        <field id="zdsdy"          long_name="FK S y-gradient"                                     unit=" 10^-3 m^-1"        /> 
     308        <field id="dbdy_mle"          long_name="FK B y-gradient"                                     unit=" s^-2"        /> 
    254309      </field_group> 
    255310       
     
    263318          <field id="emp_oce"      long_name="Evap minus Precip over ocean"         standard_name="evap_minus_precip_over_sea_water"                                     unit="kg/m2/s"   /> 
    264319          <field id="emp_ice"      long_name="Evap minus Precip over ice"           standard_name="evap_minus_precip_over_sea_ice"                                       unit="kg/m2/s"   /> 
    265           <field id="saltflx"      long_name="Downward salt flux"                                                                                                        unit="1e-3/m2/s" /> 
    266           <field id="fmmflx"       long_name="Water flux due to freezing/melting"                                                                                        unit="kg/m2/s"   /> 
     320          <field id="saltflx"      long_name="Total downward salinity flux"             standard_name="total_downward_salinity_flux"                                   unit="1e-3/m2/s" /> 
     321          <field id="fmmflx"       long_name="Water flux due to freezing/melting"     standard_name="water_flux_due_to_freezing/melting"                  unit="kg/m2/s"   /> 
    267322          <field id="snowpre"      long_name="Snow precipitation"                   standard_name="snowfall_flux"                                                        unit="kg/m2/s"   /> 
    268323          <field id="runoffs"      long_name="River Runoffs"                        standard_name="water_flux_into_sea_water_from_rivers"                                unit="kg/m2/s"   /> 
     
    699754     <field id="strd_zdfp"     long_name="salinity   -trend: pure vert. diffusion"   unit="1e-3/s" /> 
    700755 
    701      <!-- --> 
     756     <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     757     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s" /> 
     758     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s" /> 
     759 
     760 
     761    <!-- --> 
    702762     <field id="ttrd_dmp"      long_name="temperature-trend: interior restoring"        unit="degC/s" /> 
    703763     <field id="strd_dmp"      long_name="salinity   -trend: interior restoring"        unit="1e-3/s" /> 
     
    735795     <field id="strd_zdfp_e3t"     unit="1e-3/s * m"  >  strd_zdfp * e3t </field> 
    736796 
     797          <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     798     <field id="ttrd_osm_e3t"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field> 
     799     <field id="strd_osm_e3t"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field> 
     800      
    737801     <!-- --> 
    738802     <field id="ttrd_dmp_e3t"      unit="degC/s * m"  >  ttrd_dmp * e3t </field> 
     
    742806     <field id="ttrd_npc_e3t"      unit="degC/s * m"  >  ttrd_npc * e3t </field> 
    743807     <field id="strd_npc_e3t"      unit="1e-3/s * m"  >  strd_npc * e3t </field> 
    744      <field id="ttrd_qns_e3t"      unit="degC/s * m"  >  ttrd_qns * e3ts </field> 
    745      <field id="strd_cdt_e3t"      unit="degC/s * m"  >  strd_cdt * e3ts </field> 
     808     <field id="ttrd_qns_e3t"      unit="degC/s * m"  >  ttrd_qns * e3t_surf </field> 
     809     <!-- <field id="strd_cdt_e3t"      unit="degC/s * m"  >  strd_cdt * e3t_surf </field> --> 
    746810     <field id="ttrd_qsr_e3t"      unit="degC/s * m"  >  ttrd_qsr * e3t </field> 
    747811     <field id="ttrd_bbc_e3t"      unit="degC/s * m"  >  ttrd_bbc * e3t </field> 
    748812 
    749813     <!-- OMIP  layer-integrated trends --> 
     814     <field id="sfd"      long_name="Total downward salt flux"                                                     unit="kg/(m^2 s)" > saltflux * 1026.0 * 0.001 </field> 
     815     <field id="wfo"      long_name="Total downward FW mass flux"                                                   unit="kg/(m^2 s)" > -empmr </field> 
     816 
    750817     <field id="ttrd_totad_li"    long_name="layer integrated heat-trend: total advection"         unit="W/m^2"     > ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field> 
    751818     <field id="strd_totad_li"    long_name="layer integrated salt-trend: total advection"         unit="kg/(m^2 s)"    > strd_totad_e3t * 1026.0 * 0.001  </field> 
     819     <field id="ttrd_osm_li"    long_name="layer integrated heat-trend: non-local OSM"         unit="W/m^2"     > ttrd_osm_e3t * 1026.0 * 3991.86795711963 </field> 
     820     <field id="strd_osm_li"    long_name="layer integrated salt-trend: non-local OSM"         unit="kg/(m^2 s)"    > strd_osm_e3t * 1026.0 * 0.001  </field> 
    752821     <field id="ttrd_evd_li"      long_name="layer integrated heat-trend: EVD convection"          unit="W/m^2"    > ttrd_evd_e3t * 1026.0 * 3991.86795711963 </field> 
    753822     <field id="strd_evd_li"      long_name="layer integrated salt-trend: EVD convection"          unit="kg/(m^2 s)"  > strd_evd_e3t * 1026.0 * 0.001  </field> 
     823     <field id="ttrd_npc_li"      long_name="layer integrated heat-trend: NPC convection"          unit="W/m^2"    > ttrd_npc_e3t * 1026.0 * 3991.86795711963 </field> 
     824     <field id="strd_npc_li"      long_name="layer integrated salt-trend: NPC convection"          unit="kg/(m^2 s)"  > strd_npc_e3t * 1026.0 * 0.001  </field> 
    754825     <field id="ttrd_iso_li"      long_name="layer integrated heat-trend: isopycnal diffusion"     unit="W/m^2" > ttrd_iso_e3t * 1026.0 * 3991.86795711963 </field> 
    755826     <field id="strd_iso_li"      long_name="layer integrated salt-trend: isopycnal diffusion"     unit="kg/(m^2 s)" > strd_iso_e3t * 1026.0 * 0.001  </field> 
     
    758829     <field id="ttrd_qns_li"      long_name="layer integrated heat-trend: non-solar flux + runoff" unit="W/m^2" grid_ref="grid_T_2D"> ttrd_qns_e3t * 1026.0 * 3991.86795711963 </field> 
    759830     <field id="ttrd_qsr_li"      long_name="layer integrated heat-trend: solar flux"              unit="W/m^2"  grid_ref="grid_T_3D"> ttrd_qsr_e3t * 1026.0 * 3991.86795711963 </field> 
     831     <field id="ttrd_bbc_li"      long_name="layer integrated heat-trend: geothermal heating "  unit="W/m^2" > ttrd_bbc_e3t * 1026.0 * 3991.86795711963 </field> 
    760832     <field id="ttrd_bbl_li"      long_name="layer integrated heat-trend: bottom boundary layer "  unit="W/m^2" > ttrd_bbl_e3t * 1026.0 * 3991.86795711963 </field> 
    761833     <field id="strd_bbl_li"      long_name="layer integrated salt-trend: bottom boundary layer "  unit="kg/(m^2 s)" > strd_bbl_e3t * 1026.0 * 0.001  </field> 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/SHARED/namelist_ref

    r12658 r12912  
    257257   ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    258258   ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
     259 
     260   ln_humi_dpt    = .false.   !  Supply dewpoint tempearture instead of specific humidity (is true for ERA5) 
     261 
    259262      ! 
    260263      rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
     
    10881091&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
    10891092!----------------------------------------------------------------------- 
    1090    ln_use_osm_la = .false.      !  Use namelist  rn_osm_la 
     1093   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    10911094   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    10921095   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     
    11031106      !                        !  = 1: Pierson Moskowitz wave spectrum 
    11041107      !                        !  = 0: Constant La# = 0.3 
     1108   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
     1109/ 
     1110!----------------------------------------------------------------------- 
     1111&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     1112!----------------------------------------------------------------------- 
     1113   rn_osm_mle_ce       = 0.06      ! magnitude of the MLE (typical value: 0.06 to 0.08) 
     1114   nn_osm_mle          = 0         ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 
     1115   rn_osm_mle_lf       = 5.e+3     ! typical scale of mixed layer front (meters)                      (case rn_osm_mle=0) 
     1116   rn_osm_mle_time     = 172800.   ! time scale for mixing momentum across the mixed layer (seconds)  (case rn_osm_mle=0) 
     1117   rn_osm_mle_lat      = 20.       ! reference latitude (degrees) of MLE coef.                        (case rn_mle=1) 
     1118   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
     1119   rn_osm_mle_thresh = 0.0005      ! delta b criterion used for FK MLD criterion 
     1120   rn_osm_mle_tau     = 172800.   ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
    11051121/ 
    11061122!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/cfgs/ref_cfgs.txt

    r12658 r12912  
    99ORCA2_ICE_PISCES OCE TOP ICE NST 
    1010SPITZ12 OCE ICE 
     11eORCA1 OCE ICE 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/DIA/diawri.F90

    r12658 r12912  
    4343   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
    4444   USE zdfmxl         ! mixed layer 
     45   USE zdfosm         ! mixed layer 
    4546   ! 
    4647   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    923924         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    924925      ENDIF 
     926 
     927      IF( ln_zdfosm ) THEN 
     928         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)   )    ! now boundary-layer depth 
     929         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)    )    ! now mixed-layer depth 
     930         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask       )    ! w-level diffusion 
     931         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask       )    ! now w-level viscosity 
     932         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask       )    ! non-local t forcing 
     933         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask       )    ! non-local s forcing 
     934         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask       )    ! non-local u forcing 
     935         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask       )    ! non-local v forcing 
     936         IF( ln_osm_mle ) THEN 
     937            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)   )    ! now transition-layer depth 
     938         END IF 
     939      ENDIF 
    925940  
    926941#if defined key_si3 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/SBC/sbcblk.F90

    r12658 r12912  
    124124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    125125 
     126   LOGICAL  ::   ln_humi_dpt = .FALSE.                                        ! calculate specific hunidity from dewpoint 
     127   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qair                      ! specific humidity of air at input height 
     128 
    126129   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    127130   !                            ! associated indices: 
     
    145148      !!------------------------------------------------------------------- 
    146149      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    147          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     150         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), qair(jpi,jpj), STAT=sbc_blk_alloc ) 
    148151      ! 
    149152      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 
     
    171174      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    172175         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    173          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
     176         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, ln_humi_dpt,&   ! bulk algorithm 
    174177         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    175178         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
     
    323326      ! 
    324327      !                                            ! compute the surface ocean fluxes using bulk formulea 
     328      ! ..... if dewpoint supplied instead of specific humidaity, calculate specific humidity 
     329      IF(ln_humi_dpt) THEN 
     330         qair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     331      ELSE 
     332         qair(:,:) = sf(jp_humi)%fnow(:,:,1) 
     333      END IF 
     334       
    325335      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 
    326336 
     
    332342         ENDIF  
    333343         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    334          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
     344         qatm_ice(:,:)    = qair(:,:) 
    335345         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    336346         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     
    434444      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    435445      !!    (since reanalysis products provide T at z, not theta !) 
    436       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 
     446      ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), qair(:,:) ) * rn_zqt 
    437447 
    438448      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    439449      ! 
    440       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
     450      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! NCAR-COREv2 
    441451         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
     452      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.0 
    443453         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
     454      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.5 
    445455         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
     456      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! ECMWF 
    447457         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    448458      CASE DEFAULT 
     
    454464         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    455465      ELSE                                      ! At zt: 
    456          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     466         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    457467      END IF 
    458468 
     
    495505      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    496506         !! q_air and t_air are given at 10m (wind reference height) 
    497          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    498          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
     507         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - qair(:,:)) ) ! Evaporation, using bulk wind speed 
     508         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    499509      ELSE 
    500510         !! q_air and t_air are not given at 10m (wind reference height) 
    501511         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    502512         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
    503          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
     513         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    504514      ENDIF 
    505515 
     
    742752      ! local scalars ( place there for vector optimisation purposes) 
    743753      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    744       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
     754      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1)) 
    745755 
    746756      !!gm brutal.... 
     
    807817      zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    808818      ! 
    809       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     819      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    810820      ! 
    811821      zztmp = 1. / ( 1. - albo ) 
     
    838848               ! Latent Heat 
    839849               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    840                   &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     850                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - qair(ji,jj) ) ) 
    841851               ! Latent heat sensitivity for ice (Dqla/Dt) 
    842852               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRA/tramle.F90

    r12658 r12912  
    2121   USE lbclnk         ! lateral boundary condition / mpp link 
    2222 
     23   ! where OSMOSIS_OBL is used with integrated FK 
     24   USE zdf_oce, ONLY : ln_zdfosm 
     25   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 
     26 
    2327   IMPLICIT NONE 
    2428   PRIVATE 
     
    5660CONTAINS 
    5761 
    58    SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 
    59       !!---------------------------------------------------------------------- 
    60       !!                  ***  ROUTINE tra_mle_trp  *** 
    61       !! 
    62       !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport 
    63       !! 
    64       !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced 
    65       !!              transport are computed as follows : 
    66       !!                zu_mle = dk[ zpsi_uw ] 
    67       !!                zv_mle = dk[ zpsi_vw ] 
    68       !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 
    69       !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc) 
    70       !!              and added to the input velocity : 
    71       !!                p.n = p.n + z._mle 
    72       !! 
    73       !! ** Action  : - (pun,pvn,pwn) increased by the mle transport 
    74       !!                CAUTION, the transport is not updated at the last line/raw 
    75       !!                         this may be a problem for some advection schemes 
    76       !! 
    77       !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
    78       !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
    79       !!---------------------------------------------------------------------- 
    80       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    81       INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
    82       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    83       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
    84       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
    85       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
    86       ! 
    87       INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    88       INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
    89       REAL(wp) ::   zcuw, zmuw, zc      ! local scalar 
    90       REAL(wp) ::   zcvw, zmvw          !   -      - 
    91       INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
    92       REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 
    94       !!---------------------------------------------------------------------- 
    95       ! 
    96       !                                      !==  MLD used for MLE  ==! 
    97       !                                                ! compute from the 10m density to deal with the diurnal cycle 
    98       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1) 
    99       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
    100          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m) 
    101             DO jj = 1, jpj 
    102                DO ji = 1, jpi                          ! index of the w-level at the ML based 
    103                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
    104                END DO 
    105             END DO 
    106          END DO 
    107       ENDIF 
    108       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
    109       ! 
    110       ! 
    111       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
    112       zbm (:,:) = 0._wp 
    113       zn2 (:,:) = 0._wp 
    114       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
    115          DO jj = 1, jpj 
    116             DO ji = 1, jpi 
    117                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
    118                zmld(ji,jj) = zmld(ji,jj) + zc 
    119                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 
    120                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
    121             END DO 
    122          END DO 
    123       END DO 
    124  
    125       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    126       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
    127          DO jj = 1, jpjm1 
    128             DO ji = 1, fs_jpim1   ! vector opt. 
    129                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
    130                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
    131             END DO 
    132          END DO 
    133       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
    134          DO jj = 1, jpjm1 
    135             DO ji = 1, fs_jpim1   ! vector opt. 
    136                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
    137                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
    138             END DO 
    139          END DO 
    140       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
    141          DO jj = 1, jpjm1 
    142             DO ji = 1, fs_jpim1   ! vector opt. 
    143                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
    144                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
    145             END DO 
    146          END DO 
    147       END SELECT 
    148       !                                                ! convert density into buoyancy 
    149       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
    150       ! 
    151       ! 
    152       !                                      !==  Magnitude of the MLE stream function  ==! 
    153       ! 
    154       !                 di[bm]  Ds 
    155       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
    156       !                  e1u   Lf fu                      and the e2u for the "transport" 
    157       !                                                      (not *e3u as divided by e3u at the end) 
    158       ! 
    159       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
    160          DO jj = 1, jpjm1 
    161             DO ji = 1, fs_jpim1   ! vector opt. 
    162                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
    163                   &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
    164                   &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
    165                   ! 
    166                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
    167                   &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
    168                   &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
    169             END DO 
    170          END DO 
    171          ! 
    172       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
    173          DO jj = 1, jpjm1 
    174             DO ji = 1, fs_jpim1   ! vector opt. 
    175                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
    176                   &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    177                   ! 
    178                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
    179                   &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
    180             END DO 
    181          END DO 
    182       ENDIF 
    183       ! 
    184       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
    185          DO jj = 1, jpjm1 
    186             DO ji = 1, fs_jpim1   ! vector opt. 
    187                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
    188                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
    189             END DO 
    190          END DO 
    191       ENDIF 
    192       ! 
    193       !                                      !==  structure function value at uw- and vw-points  ==! 
    194       DO jj = 1, jpjm1 
    195          DO ji = 1, fs_jpim1   ! vector opt. 
    196             zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu 
    197             zhv(ji,jj) = 1._wp / zhv(ji,jj) 
    198          END DO 
    199       END DO 
    200       ! 
    201       zpsi_uw(:,:,:) = 0._wp 
    202       zpsi_vw(:,:,:) = 0._wp 
    203       ! 
    204       DO jk = 2, ikmax                                ! start from 2 : surface value = 0 
    205          DO jj = 1, jpjm1 
    206             DO ji = 1, fs_jpim1   ! vector opt. 
    207                zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 
    208                zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 
    209                zcuw = zcuw * zcuw 
    210                zcvw = zcvw * zcvw 
    211                zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  ) 
    212                zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  ) 
    213                ! 
    214                zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 
    215                zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 
    216             END DO 
    217          END DO 
    218       END DO 
    219       ! 
    220       !                                      !==  transport increased by the MLE induced transport ==! 
    221       DO jk = 1, ikmax 
    222          DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1 
    223             DO ji = 1, fs_jpim1   ! vector opt. 
    224                pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
    225                pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
    226             END DO 
    227          END DO 
    228          DO jj = 2, jpjm1 
    229             DO ji = fs_2, fs_jpim1   ! vector opt. 
    230                pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
     62  SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 
     63    !!---------------------------------------------------------------------- 
     64    !!                  ***  ROUTINE tra_mle_trp  *** 
     65    !! 
     66    !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport 
     67    !! 
     68    !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced 
     69    !!              transport are computed as follows : 
     70    !!                zu_mle = dk[ zpsi_uw ] 
     71    !!                zv_mle = dk[ zpsi_vw ] 
     72    !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 
     73    !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc) 
     74    !!              and added to the input velocity : 
     75    !!                p.n = p.n + z._mle 
     76    !! 
     77    !! ** Action  : - (pun,pvn,pwn) increased by the mle transport 
     78    !!                CAUTION, the transport is not updated at the last line/raw 
     79    !!                         this may be a problem for some advection schemes 
     80    !! 
     81    !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     82    !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     83    !!---------------------------------------------------------------------- 
     84    INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     85    INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
     86    CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     87    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
     88    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
     89    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
     90    ! 
     91    INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     92    INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     93    REAL(wp) ::   zcuw, zmuw, zc      ! local scalar 
     94    REAL(wp) ::   zcvw, zmvw          !   -      - 
     95    INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     96    REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 
     97    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 
     98    !!---------------------------------------------------------------------- 
     99    ! 
     100    ! 
     101    IF(ln_osm_mle.and.ln_zdfosm) THEN 
     102       ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     103       ! 
     104       ! 
     105       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     106       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     107          DO jj = 1, jpjm1 
     108             DO ji = 1, fs_jpim1   ! vector opt. 
     109                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 
     110                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 
     111             END DO 
     112          END DO 
     113       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     114          DO jj = 1, jpjm1 
     115             DO ji = 1, fs_jpim1   ! vector opt. 
     116                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
     117                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
     118             END DO 
     119          END DO 
     120       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     121          DO jj = 1, jpjm1 
     122             DO ji = 1, fs_jpim1   ! vector opt. 
     123                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
     124                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
     125             END DO 
     126          END DO 
     127       END SELECT 
     128       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     129          DO jj = 1, jpjm1 
     130             DO ji = 1, fs_jpim1   ! vector opt. 
     131                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            & 
     132                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     133                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     134                ! 
     135                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            & 
     136                     &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     137                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     138             END DO 
     139          END DO 
     140          ! 
     141       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     142          DO jj = 1, jpjm1 
     143             DO ji = 1, fs_jpim1   ! vector opt. 
     144                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               & 
     145                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     146                ! 
     147                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               & 
     148                     &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     149             END DO 
     150          END DO 
     151       ENDIF 
     152 
     153    ELSE !do not use osn_mle 
     154       !                                      !==  MLD used for MLE  ==! 
     155       !                                                ! compute from the 10m density to deal with the diurnal cycle 
     156       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1) 
     157       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
     158          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m) 
     159             DO jj = 1, jpj 
     160                DO ji = 1, jpi                          ! index of the w-level at the ML based 
     161                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
     162                END DO 
     163             END DO 
     164          END DO 
     165       ENDIF 
     166       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
     167 
     168       ! 
     169       ! 
     170       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
     171       zbm (:,:) = 0._wp 
     172       zn2 (:,:) = 0._wp 
     173       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
     174          DO jj = 1, jpj 
     175             DO ji = 1, jpi 
     176                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     177                zmld(ji,jj) = zmld(ji,jj) + zc 
     178                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 
     179                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
     180             END DO 
     181          END DO 
     182       END DO 
     183 
     184       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     185       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     186          DO jj = 1, jpjm1 
     187             DO ji = 1, fs_jpim1   ! vector opt. 
     188                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
     189                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
     190             END DO 
     191          END DO 
     192       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     193          DO jj = 1, jpjm1 
     194             DO ji = 1, fs_jpim1   ! vector opt. 
     195                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
     196                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
     197             END DO 
     198          END DO 
     199       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     200          DO jj = 1, jpjm1 
     201             DO ji = 1, fs_jpim1   ! vector opt. 
     202                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
     203                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
     204             END DO 
     205          END DO 
     206       END SELECT 
     207       !                                                ! convert density into buoyancy 
     208       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     209       ! 
     210       ! 
     211       !                                      !==  Magnitude of the MLE stream function  ==! 
     212       ! 
     213       !                 di[bm]  Ds 
     214       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
     215       !                  e1u   Lf fu                      and the e2u for the "transport" 
     216       !                                                      (not *e3u as divided by e3u at the end) 
     217       ! 
     218       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     219          DO jj = 1, jpjm1 
     220             DO ji = 1, fs_jpim1   ! vector opt. 
     221                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
     222                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     223                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     224                ! 
     225                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
     226                     &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     227                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     228             END DO 
     229          END DO 
     230          ! 
     231       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     232          DO jj = 1, jpjm1 
     233             DO ji = 1, fs_jpim1   ! vector opt. 
     234                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
     235                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     236                ! 
     237                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
     238                     &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     239             END DO 
     240          END DO 
     241       ENDIF 
     242       ! 
     243       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
     244          DO jj = 1, jpjm1 
     245             DO ji = 1, fs_jpim1   ! vector opt. 
     246                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
     247                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
     248             END DO 
     249          END DO 
     250       ENDIF 
     251       ! 
     252    ENDIF  ! end of lm_osm_mle loop 
     253    !                                      !==  structure function value at uw- and vw-points  ==! 
     254    DO jj = 1, jpjm1 
     255       DO ji = 1, fs_jpim1   ! vector opt. 
     256          zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu 
     257          zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall)  
     258       END DO 
     259    END DO 
     260    ! 
     261    zpsi_uw(:,:,:) = 0._wp 
     262    zpsi_vw(:,:,:) = 0._wp 
     263    ! 
     264    DO jk = 2, ikmax                                ! start from 2 : surface value = 0 
     265       DO jj = 1, jpjm1 
     266          DO ji = 1, fs_jpim1   ! vector opt. 
     267             zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 
     268             zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 
     269             zcuw = zcuw * zcuw 
     270             zcvw = zcvw * zcvw 
     271             zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  ) 
     272             zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  ) 
     273             ! 
     274             zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 
     275             zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 
     276          END DO 
     277       END DO 
     278    END DO 
     279    ! 
     280    !                                      !==  transport increased by the MLE induced transport ==! 
     281    DO jk = 1, ikmax 
     282       DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1 
     283          DO ji = 1, fs_jpim1   ! vector opt. 
     284             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
     285             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
     286          END DO 
     287       END DO 
     288       DO jj = 2, jpjm1 
     289          DO ji = fs_2, fs_jpim1   ! vector opt. 
     290             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
    231291                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 
    232             END DO 
    233          END DO 
    234       END DO 
    235  
    236       IF( cdtype == 'TRA') THEN              !==  outputs  ==! 
    237          ! 
    238          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
    239          CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
    240          ! 
    241          ! divide by cross distance to give streamfunction with dimensions m^2/s 
    242          DO jk = 1, ikmax+1 
    243             zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 
    244             zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 
    245          END DO 
    246          CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
    247          CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction 
    248       ENDIF 
    249       ! 
    250    END SUBROUTINE tra_mle_trp 
     292          END DO 
     293       END DO 
     294    END DO 
     295 
     296    IF( cdtype == 'TRA') THEN              !==  outputs  ==! 
     297       ! 
     298       IF (ln_osm_mle.and.ln_zdfosm) THEN 
     299          zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
     300       ELSE 
     301          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
     302       END IF 
     303       CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
     304       ! 
     305       ! divide by cross distance to give streamfunction with dimensions m^2/s 
     306       DO jk = 1, ikmax+1 
     307          zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 
     308          zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 
     309       END DO 
     310       CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
     311       CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction 
     312    ENDIF 
     313    ! 
     314  END SUBROUTINE tra_mle_trp 
    251315 
    252316 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRD/trd_oce.F90

    r12658 r12912  
    3333# endif 
    3434   !                                                  !!!* Active tracers trends indexes 
    35    INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 20     !: Total trend nb: change it when adding/removing one indice below 
     35   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 21     !: Total trend nb: change it when adding/removing one indice below 
    3636   !                               ===============     !   
    3737   INTEGER, PUBLIC, PARAMETER ::   jptra_xad  =  1     !: x- horizontal advection 
     
    4646   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  = 10     !: Bottom Boundary Condition (geoth. heating)  
    4747   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  = 11     !: Bottom Boundary Layer (diffusive and/or advective) 
     48   INTEGER, PUBLIC, PARAMETER ::   jptra_osm  = 21     !: Non-local terms from OSMOSIS OBL model 
    4849   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 12     !: non-penetrative convection treatment 
    4950   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 13     !: internal restoring (damping) 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRD/trdtra.F90

    r12658 r12912  
    346346         CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl"  , ptrdx )        ! bottom boundary layer 
    347347                                  CALL iom_put( "strd_bbl"  , ptrdy ) 
     348         CASE( jptra_osm  )   ;   CALL iom_put( "ttrd_osm"  , ptrdx )        ! OSMOSIS non-local forcing 
     349                                  CALL iom_put( "strd_osm"  , ptrdy ) 
    348350         CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc"  , ptrdx )        ! static instability mixing 
    349351                                  CALL iom_put( "strd_npc"  , ptrdy ) 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF/zdfosm.F90

    r12658 r12912  
    2525   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 
    2626   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 
    27    !!            (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out). 
     27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 
    2828   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 
    2929   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 
    3030   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 
     31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, 
     32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers 
     33   !!                  (b) The stable OSBL depth parametrization changed. 
     34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. 
     35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 
    3136   !!---------------------------------------------------------------------- 
    3237 
     
    4045   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD) 
    4146   !!   dyn_osm       : compute and add to u & v trensd the non-local flux 
     47   !! 
     48   !! Subroutines in revised code. 
    4249   !!---------------------------------------------------------------------- 
    4350   USE oce            ! ocean dynamics and active tracers 
     
    6976   PUBLIC   tra_osm       ! routine called by step.F90 
    7077   PUBLIC   trc_osm       ! routine called by trcstp.F90 
    71    PUBLIC   dyn_osm       ! routine called by 'step.F90' 
     78   PUBLIC   dyn_osm       ! routine called by step.F90 
     79 
     80   PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90 
    7281 
    7382   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux 
     
    7786   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt 
    7887   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth 
    79    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbli     !: intial boundary layer depth for stable blayer 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline 
     89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth 
    8090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift. 
     91 
     92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts 
     93   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML 
     96   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer. 
    8197 
    8298   !                      !!** Namelist  namzdf_osm  ** 
    8399   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la 
     100 
     101   LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 
     102 
    84103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    85104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
     
    96115   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s) 
    97116 
     117! OSMOSIS mixed layer eddy parametrization constants 
     118   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt 
     119   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient 
     120   !                                        ! parameters used in nn_osm_mle = 0 case 
     121   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front 
     122   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer 
     123   !                                        ! parameters used in nn_osm_mle = 1 case 
     124   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
     125   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
     126   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
     127   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 
     128   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
     129   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
     130   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
     131 
     132 
    98133   !                                    !!! ** General constants  ** 
    99    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number 
     134   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
     135   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    100136   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    101137   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
     
    105141   !!---------------------------------------------------------------------- 
    106142   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    107    !! $Id$ 
     143   !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $ 
    108144   !! Software governed by the CeCILL license (see ./LICENSE) 
    109145   !!---------------------------------------------------------------------- 
     
    114150      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    115151      !!---------------------------------------------------------------------- 
    116      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 
    117           &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     & 
    118           &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc ) 
     152     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
     153          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     154          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     155 
     156     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 
     157          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 
     158 
     159!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ? 
     160!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     161!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     162!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
     163!     IF ( ln_osm_mle ) THEN 
     164!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 
     165!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 
     166!     ENDIF 
     167 
    119168     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
    120169     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     
    161210      !! 
    162211      INTEGER ::   ji, jj, jk                   ! dummy loop indices 
     212 
     213      INTEGER ::   jl                   ! dummy loop indices 
     214 
    163215      INTEGER ::   ikbot, jkmax, jkm1, jkp2     ! 
    164216 
     
    166218      REAL(wp) ::   zbeta, zthermal                                  ! 
    167219      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 
    168       REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       ! 
     220      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      ! 
     221 
    169222      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density 
    170223      INTEGER  ::   jm                          ! dummy loop indices 
     
    191244      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    192245      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     246 
     247 
     248      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL 
     249      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux 
     250      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 
     251      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK 
     252 
    193253      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift 
    194254      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number 
     
    196256      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 
    197257      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    198       LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 
     258      LOGICAL, DIMENSION(jpi,jpj)  :: lconv    ! unstable/stable bl 
    199259 
    200260      ! mixed-layer variables 
     
    208268      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    209269      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid 
     270 
     271      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 
     272      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid 
     273 
    210274      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    211275      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
     276      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure. 
     277      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients 
     278 
     279      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
     280 
    212281      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    213282      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
     
    238307      ! Temporary variables 
    239308      INTEGER :: inhml 
    240       INTEGER :: i_lconv_alloc 
    241309      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 
    242310      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables 
     
    248316      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    249317 
     318      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
     319      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 
     320      REAL(wp) :: zzeta_s = 0._wp 
     321      REAL(wp) :: zzeta_v = 0.46 
     322      REAL(wp) :: zabsstke 
     323 
    250324      ! For debugging 
    251325      INTEGER :: ikt 
    252326      !!-------------------------------------------------------------------- 
    253327      ! 
    254       ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc ) 
    255       IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv') 
    256  
    257328      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    258329      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
     
    268339      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    269340      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     341 
    270342      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    271343      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     
    277349      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    278350      ! 
     351      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 
     352 
     353      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
     354         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 
     355         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 
     356         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 
     357         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp 
     358      ENDIF 
     359      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
     360 
    279361      ! Flux-Gradient arrays. 
    280362      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
     
    287369      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    288370 
     371      zdhdt_2(:,:) = 0._wp 
    289372      ! hbl = MAX(hbl,epsln) 
    290373      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    350433              ! Use wind speed wndm included in sbc_oce module 
    351434              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    352               dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
     435              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    353436           END DO 
    354437        END DO 
     
    362445              ! It could represent the effects of the spread of wave directions 
    363446              ! around the mean wind. The effect of this adjustment needs to be tested. 
    364               zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
    365                    &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
    366               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
     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 ! 
    367450           END DO 
    368451        END DO 
     
    375458           ! Langmuir velocity scale (zwstrl), at T-point 
    376459           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
    377            ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
    378            ! Intended as a possible test. Doesn't affect LES results for entrainment, 
    379            !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
    380            zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
    381                 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
    382                 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
    383            ! define La this way so effects of Stokes penetration depth on velocity scale are included 
    384            zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
     460           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 
     461           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 
    385462           ! Velocity scale that tends to zustar for large Langmuir numbers 
    386463           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
     
    389466           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    390467           ! Note zustke and zwstrl are not amended. 
    391            IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    392468           ! 
    393469           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
     
    406482     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 
    407483     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408      ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
    410       ibld(:,:) = 3 
    411       DO jk = 4, jpkm1 
    412          DO jj = 2, jpjm1 
    413             DO ji = 2, jpim1 
     484     ! BL must be always 4 levels deep. 
     485     ! For calculation of lateral buoyancy gradients for FK in 
     486     ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 
     487     ! previously exist for hbl also. 
     488      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
     489      ibld(:,:) = 4 
     490      DO jk = 5, jpkm1 
     491         DO jj = 1, jpj 
     492            DO ji = 1, jpi 
    414493               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    415494                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
     
    419498      END DO 
    420499 
    421       DO jj = 2, jpjm1                                 !  Vertical slab 
     500      DO jj = 2, jpjm1 
    422501         DO ji = 2, jpim1 
    423                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    424                zbeta    = rab_n(ji,jj,1,jp_sal) 
    425                zt   = 0._wp 
    426                zs   = 0._wp 
    427                zu   = 0._wp 
    428                zv   = 0._wp 
    429                ! average over depth of boundary layer 
    430                zthick=0._wp 
    431                DO jm = 2, ibld(ji,jj) 
    432                   zthick=zthick+e3t_n(ji,jj,jm) 
    433                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    434                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    435                   zu   = zu  + e3t_n(ji,jj,jm) & 
    436                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    437                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    438                   zv   = zv  + e3t_n(ji,jj,jm) & 
    439                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    440                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    441                END DO 
    442                zt_bl(ji,jj) = zt / zthick 
    443                zs_bl(ji,jj) = zs / zthick 
    444                zu_bl(ji,jj) = zu / zthick 
    445                zv_bl(ji,jj) = zv / zthick 
    446                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    447                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    448                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    449                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    450                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    451                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    452                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    453                IF ( lconv(ji,jj) ) THEN    ! Convective 
    454                       zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    455                            &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    456  
    457                       zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    458                            &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    460 !                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    461 !                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    462  
    463 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    464 !                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    465                       zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    466                ELSE                        ! Stable 
    467                       zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
    468                            &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
    469                            & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
    470                            &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
    471                       zzdhdt = zzdhdt + zwbav(ji,jj) 
    472                       IF ( zzdhdt < 0._wp ) THEN 
    473                       ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    474                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    475                       ELSE 
    476                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    477                               &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    478                       ENDIF 
    479                       zzdhdt = 2.0 * zzdhdt / zpert 
    480                ENDIF 
    481                zdhdt(ji,jj) = zzdhdt 
    482            END DO 
     502            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     503            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 
     504            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     505         END DO 
    483506      END DO 
    484  
     507      ! Averages over well-mixed and boundary layer 
     508      ! Alan: do we need zb_nl?, zb_ml? 
     509      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     510      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     511! External gradient 
     512      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     513 
     514 
     515      IF ( ln_osm_mle ) THEN 
     516         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     517         CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     518       ENDIF 
     519 
     520! Rate of change of hbl 
     521      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
    485522      ! Calculate averages over depth of boundary layer 
    486       imld = ibld           ! use imld to hold previous blayer index 
    487       ibld(:,:) = 3 
    488  
    489       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
    490       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    491       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     523         DO jj = 2, jpjm1 
     524            DO ji = 2, jpim1 
     525               zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 
     526               ! adjustment to represent limiting by ocean bottom 
     527               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
     528            END DO 
     529         END DO 
     530 
     531      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index 
     532      ibld(:,:) = 4 
    492533 
    493534      DO jk = 4, jpkm1 
     
    495536            DO ji = 2, jpim1 
    496537               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    497                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     538                  ibld(ji,jj) = jk 
    498539               ENDIF 
    499540            END DO 
     
    504545! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    505546! 
    506       DO jj = 2, jpjm1 
    507          DO ji = 2, jpim1 
    508             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     547      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     548      ! Alan: do we need zb_ml? 
     549      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    509550! 
    510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    511551! 
    512                zhbl_s = hbl(ji,jj) 
    513                jm = imld(ji,jj) 
    514                zthermal = rab_n(ji,jj,1,jp_tem) 
    515                zbeta = rab_n(ji,jj,1,jp_sal) 
    516                IF ( lconv(ji,jj) ) THEN 
    517 !unstable 
    518                   zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    519                        &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    520  
    521                   DO jk = imld(ji,jj), ibld(ji,jj) 
    522                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
    523                           & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524  
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
    526                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527  
    528                      IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
    529                   END DO 
    530                   hbl(ji,jj) = zhbl_s 
    531                   ibld(ji,jj) = jm 
    532                   hbli(ji,jj) = hbl(ji,jj) 
    533                ELSE 
    534 ! stable 
    535                   DO jk = imld(ji,jj), ibld(ji,jj) 
    536                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          & 
    537                           &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 
    538                           & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    539  
    540                      zhbl_s = zhbl_s +  (                                                                                & 
    541                           &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
    542                           &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
    543                           &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    544                           &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    545                           &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here 
    546  
    547                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    548                      IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
    549                   END DO 
    550                   hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 
    551                   ibld(ji,jj) = MAX(jm, 3 ) 
    552                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    553                ENDIF   ! IF ( lconv ) 
    554             ELSE 
    555 ! change zero or one model level. 
    556                hbl(ji,jj) = zhbl_t(ji,jj) 
    557                IF ( lconv(ji,jj) ) THEN 
    558                   hbli(ji,jj) = hbl(ji,jj) 
    559                ELSE 
    560                   hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 
    561                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    562                ENDIF 
    563             ENDIF 
    564             zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    565          END DO 
    566       END DO 
     552      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
    567553      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    568  
    569 ! Recalculate averages over boundary layer after depth updated 
    570      ! Consider later  combining this into the loop above and looking for columns 
    571      ! where the index for base of the boundary layer have changed 
    572       DO jj = 2, jpjm1                                 !  Vertical slab 
    573          DO ji = 2, jpim1 
    574                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    575                zbeta    = rab_n(ji,jj,1,jp_sal) 
    576                zt   = 0._wp 
    577                zs   = 0._wp 
    578                zu   = 0._wp 
    579                zv   = 0._wp 
    580                ! average over depth of boundary layer 
    581                zthick=0._wp 
    582                DO jm = 2, ibld(ji,jj) 
    583                   zthick=zthick+e3t_n(ji,jj,jm) 
    584                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    585                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    586                   zu   = zu  + e3t_n(ji,jj,jm) & 
    587                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    588                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    589                   zv   = zv  + e3t_n(ji,jj,jm) & 
    590                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    591                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    592                END DO 
    593                zt_bl(ji,jj) = zt / zthick 
    594                zs_bl(ji,jj) = zs / zthick 
    595                zu_bl(ji,jj) = zu / zthick 
    596                zv_bl(ji,jj) = zv / zthick 
    597                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    598                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    599                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    600                       &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    601                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    602                       &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    603                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    604                zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    605                IF ( lconv(ji,jj) ) THEN 
    606                   IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    607                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    608                            zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
    609                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    610                      ELSE                                                     ! unstable 
    611                            zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
    612                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    613                      ENDIF 
    614                      IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
    615                         zari = 0.2 
    616                         zwb_ent(ji,jj) = 0._wp 
    617                      ENDIF 
    618                      inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    619                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    620                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    621                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    622                   ELSE  ! IF (zdb_bl) 
    623                      imld(ji,jj) = ibld(ji,jj) - 1 
    624                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    625                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    626                   ENDIF 
    627                ELSE   ! IF (lconv) 
    628                   IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    629                   ! boundary layer deepening 
    630                      IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    631                   ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    632                         zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    633                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    634                         inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    635                         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    636                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    637                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    638                      ELSE 
    639                         imld(ji,jj) = ibld(ji,jj) - 1 
    640                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    641                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    642                      ENDIF ! IF (zdb_bl > 0.0) 
    643                   ELSE     ! IF(dhdt >= 0) 
    644                   ! boundary layer collapsing. 
    645                      imld(ji,jj) = ibld(ji,jj) 
    646                      zhml(ji,jj) = zhbl(ji,jj) 
    647                      zdh(ji,jj) = 0._wp 
    648                   ENDIF    ! IF (dhdt >= 0) 
    649                ENDIF       ! IF (lconv) 
    650          END DO 
    651       END DO 
    652  
    653       ! Average over the depth of the mixed layer in the convective boundary layer 
    654       ! Also calculate entrainment fluxes for temperature and salinity 
    655       DO jj = 2, jpjm1                                 !  Vertical slab 
    656          DO ji = 2, jpim1 
    657             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    658             zbeta    = rab_n(ji,jj,1,jp_sal) 
    659             IF ( lconv(ji,jj) ) THEN 
    660                zt   = 0._wp 
    661                zs   = 0._wp 
    662                zu   = 0._wp 
    663                zv   = 0._wp 
    664                ! average over depth of boundary layer 
    665                zthick=0._wp 
    666                DO jm = 2, imld(ji,jj) 
    667                   zthick=zthick+e3t_n(ji,jj,jm) 
    668                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    669                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    670                   zu   = zu  + e3t_n(ji,jj,jm) & 
    671                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    672                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    673                   zv   = zv  + e3t_n(ji,jj,jm) & 
    674                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    675                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    676                END DO 
    677                zt_ml(ji,jj) = zt / zthick 
    678                zs_ml(ji,jj) = zs / zthick 
    679                zu_ml(ji,jj) = zu / zthick 
    680                zv_ml(ji,jj) = zv / zthick 
    681                zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    682                zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    683                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    684                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    685                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    686                      &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    687                zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    688             ELSE 
    689             ! stable, if entraining calulate average below interface layer. 
    690                IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    691                   zt   = 0._wp 
    692                   zs   = 0._wp 
    693                   zu   = 0._wp 
    694                   zv   = 0._wp 
    695                   ! average over depth of boundary layer 
    696                   zthick=0._wp 
    697                   DO jm = 2, imld(ji,jj) 
    698                      zthick=zthick+e3t_n(ji,jj,jm) 
    699                      zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    700                      zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    701                      zu   = zu  + e3t_n(ji,jj,jm) & 
    702                         &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    703                         &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    704                      zv   = zv  + e3t_n(ji,jj,jm) & 
    705                         &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    706                         &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    707                   END DO 
    708                   zt_ml(ji,jj) = zt / zthick 
    709                   zs_ml(ji,jj) = zs / zthick 
    710                   zu_ml(ji,jj) = zu / zthick 
    711                   zv_ml(ji,jj) = zv / zthick 
    712                   zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    713                   zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    714                   zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    715                         &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    716                   zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    717                         &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    718                   zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    719                ENDIF 
    720             ENDIF 
    721          END DO 
    722       END DO 
    723     ! 
     554! 
     555    ! Average over the depth of the mixed layer in the convective boundary layer 
     556    ! Alan: do we need zb_ml? 
     557    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    724558    ! rotate mean currents and changes onto wind align co-ordinates 
    725559    ! 
    726  
    727       DO jj = 2, jpjm1 
    728          DO ji = 2, jpim1 
    729             ztemp = zu_ml(ji,jj) 
    730             zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
    731             zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    732             ztemp = zdu_ml(ji,jj) 
    733             zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
    734             zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    735     ! 
    736             ztemp = zu_bl(ji,jj) 
    737             zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
    738             zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    739             ztemp = zdu_bl(ji,jj) 
    740             zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
    741             zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    742          END DO 
    743       END DO 
    744  
     560    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     561    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    745562     zuw_bse = 0._wp 
    746563     zvw_bse = 0._wp 
     564     zwth_ent = 0._wp 
     565     zws_ent = 0._wp 
    747566     DO jj = 2, jpjm1 
    748567        DO ji = 2, jpim1 
    749  
    750            IF ( lconv(ji,jj) ) THEN 
    751               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    752                  zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    753                  zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     568           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 
     569              IF ( lconv(ji,jj) ) THEN 
     570             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 
     571                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 
     572                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 
     573            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 
     574                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 
     575                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 
     576                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     577                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     578                 ENDIF 
     579              ELSE 
     580                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
     581                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    754582              ENDIF 
    755            ELSE 
    756               zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    757               zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    758583           ENDIF 
    759584        END DO 
     
    764589      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    765590 
    766        DO jj = 2, jpjm1 
    767           DO ji = 2, jpim1 
    768           ! 
    769              IF ( lconv (ji,jj) ) THEN 
    770              ! Unstable conditions 
    771                 IF( zdb_bl(ji,jj) > 0._wp ) THEN 
    772                 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
    773                    ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
    774                    zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
    775                    zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    776                    DO jk = 2 , ibld(ji,jj) 
    777                       znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    778                       zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    779                       zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    780                       zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    781                    END DO 
    782                 ENDIF 
    783              ELSE 
    784              ! stable conditions 
    785              ! if pycnocline profile only defined when depth steady of increasing. 
    786                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
    787                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    788                      IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    789                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
    790                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
    791                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    792                          DO jk = 2, ibld(ji,jj) 
    793                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    794                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    795                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    796                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    797                          END DO 
    798                      ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    799                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
    800                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
    801                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    802                          DO jk = 2, ibld(ji,jj) 
    803                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    804                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    805                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    806                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    807                          END DO 
    808                       ENDIF ! IF (zhol >=0.5) 
    809                    ENDIF    ! IF (zdb_bl> 0.) 
    810                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
    811              ENDIF          ! IF (lconv) 
    812             ! 
    813           END DO 
    814        END DO 
    815 ! 
    816        DO jj = 2, jpjm1 
    817           DO ji = 2, jpim1 
    818           ! 
    819              IF ( lconv (ji,jj) ) THEN 
    820              ! Unstable conditions 
    821                  zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
    822                & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
    823                 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
    824               & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    825                 DO jk = 2 , ibld(ji,jj)-1 
    826                    znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    827                    zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    828                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    829                 END DO 
    830              ELSE 
    831              ! stable conditions 
    832                 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    833                 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    834                 DO jk = 2, ibld(ji,jj) 
    835                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    836                    IF ( znd < 1.0 ) THEN 
    837                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    838                    ELSE 
    839                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    840                    ENDIF 
    841                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
    842                 END DO 
    843              ENDIF 
    844             ! 
    845           END DO 
    846        END DO 
     591      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     592      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 
     593      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    847594       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    848595       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    849596       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    850597 
    851       ! WHERE ( lconv ) 
    852       !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 
    853       !     zvisml_sc = zdifml_sc 
    854       !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 
    855       !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 
    856       !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 
    857       !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 
    858       !  ELSEWHERE 
    859       !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    860       !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    861       !  ENDWHERE 
    862598       DO jj = 2, jpjm1 
    863599          DO ji = 2, jpim1 
     
    868604               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    869605               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    870                zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
     606               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 
    871607             ELSE 
    872608               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    873609               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    874             END IF 
    875         END DO 
    876     END DO 
     610             END IF 
     611          END DO 
     612       END DO 
    877613! 
    878614       DO jj = 2, jpjm1 
     
    889625                ! pycnocline - if present linear profile 
    890626                IF ( zdh(ji,jj) > 0._wp ) THEN 
     627                   zgamma_b = 6.0 
    891628                   DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892629                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    893630                       ! 
    894                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     631                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    895632                       ! 
    896                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     633                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    897634                   END DO 
     635                   IF ( ibld_ext == 0 ) THEN 
     636                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     637                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     638                   ELSE 
     639                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     640                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     641                   ENDIF 
    898642                ENDIF 
    899                 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    900                 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 
     643                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
     644                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 
     645                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 
    901646                ! could be taken out, take account of entrainment represents as a diffusivity 
    902647                ! should remove w from here, represents entrainment 
     
    908653                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    909654                END DO 
     655 
     656                IF ( ibld_ext == 0 ) THEN 
     657                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     658                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     659                ELSE 
     660                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     661                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     662                ENDIF 
    910663             ENDIF   ! end if ( lconv ) 
    911 ! 
     664             ! 
    912665          END DO  ! end of ji loop 
    913666       END DO     ! end of jj loop 
     
    952705       END DO     ! end of jj loop 
    953706 
    954  
    955707! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 
    956708       WHERE ( lconv ) 
    957           zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ) 
    958           zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln ) 
    959           zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
     709          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 ) 
     710          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
     711          zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
    960712       ELSEWHERE 
    961713          zsc_uw_1 = zustar**2 
    962           zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 
     714          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    963715       ENDWHERE 
    964  
     716       IF(ln_dia_osm) THEN 
     717          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     718          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     719       END IF 
    965720       DO jj = 2, jpjm1 
    966721          DO ji = 2, jpim1 
     
    1007762                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1008763                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    1009                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     764                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
    1010765                   ! non-gradient buoyancy terms 
    1011766                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
     
    1020775          END DO   ! ji loop 
    1021776       END DO      ! jj loop 
    1022  
    1023777 
    1024778       WHERE ( lconv ) 
     
    1051805       END DO           ! jj loop 
    1052806 
     807       IF(ln_dia_osm) THEN 
     808          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
     809          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     810       END IF 
    1053811! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    1054812 
     
    1089847       END DO      ! jj loop 
    1090848 
    1091  
    1092849       WHERE ( lconv ) 
    1093850          zsc_uw_1 = zustar**2 
     
    1134891          END DO 
    1135892       END DO 
     893 
     894       IF(ln_dia_osm) THEN 
     895          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
     896          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
     897          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     898          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     899          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     900          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     901       END IF 
    1136902! 
    1137903! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
     
    1165931      END DO 
    1166932 
     933       IF(ln_dia_osm) THEN 
     934          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     935          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     936       END IF 
    1167937      ! pynocline contributions 
    1168938       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
     
    1170940       DO jj = 2, jpjm1 
    1171941          DO ji = 2, jpim1 
    1172              DO jk= 2, ibld(ji,jj) 
    1173                 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1174                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1175                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1176                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1177                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
    1178                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    1179              END DO 
    1180            END DO 
     942             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     943                DO jk= 2, ibld(ji,jj) 
     944                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     945                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     946                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     947                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     948                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
     949                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     950                END DO 
     951             END IF 
     952          END DO 
    1181953       END DO 
    1182954 
     
    1185957       DO jj=2, jpjm1 
    1186958          DO ji = 2, jpim1 
    1187              IF ( lconv(ji,jj) ) THEN 
     959            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 
    1188960               DO jk = 1, imld(ji,jj) - 1 
    1189961                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1190                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1191                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
     962                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
     963                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1192964                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1193965                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
     
    1195967               DO jk = imld(ji,jj), ibld(ji,jj) 
    1196968                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    1197                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1198                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
     969                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
     970                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1199971                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1200972                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1201973                END DO 
    1202974             ENDIF 
    1203              ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
    1204              ghams(ji,jj,ibld(ji,jj)) = 0._wp 
    1205              ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
    1206              ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
     975 
     976             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     977             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     978             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     979             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1207980          END DO       ! ji loop 
    1208981       END DO          ! jj loop 
    1209982 
    1210  
     983       IF(ln_dia_osm) THEN 
     984          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     985          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     986          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 
     987          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 
     988          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
     989          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
     990          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
     991       END IF 
    1211992       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1212993       ! Need to put in code for contributions that are applied explicitly to 
     
    12321013       IF(ln_dia_osm) THEN 
    12331014          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1015          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1016          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
    12341017       END IF 
    12351018 
     
    12861069       END IF ! ln_convmix = .true. 
    12871070 
     1071 
     1072  
     1073       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
     1074          DO jj = 2 , jpjm1 
     1075             DO ji = 2, jpim1 
     1076                IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 
     1077            ! Calculate MLE flux profiles 
     1078                   ! DO jk = 1, mld_prof(ji,jj) 
     1079                   !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1080                   !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 
     1081                   !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1082                   !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 
     1083                   !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1084                   ! END DO 
     1085            ! Calculate MLE flux contribution from surface fluxes 
     1086                   DO jk = 1, ibld(ji,jj) 
     1087                     znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 
     1088                     ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 
     1089                     ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 
     1090                    END DO 
     1091                    DO jk = 1, mld_prof(ji,jj) 
     1092                      znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1093                      ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 
     1094                      ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 
     1095                    END DO 
     1096            ! Viscosity for MLEs 
     1097                    DO jk = ibld(ji,jj), mld_prof(ji,jj) 
     1098                      zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 
     1099                    END DO 
     1100            ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 
     1101                    jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 
     1102                    jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 
     1103                    DO jk = mld_prof(ji,jj),  jl 
     1104                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 
     1105                            & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 
     1106                            & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 
     1107                    END DO 
     1108                ENDIF 
     1109            END DO 
     1110          END DO 
     1111       ENDIF 
     1112 
     1113       IF(ln_dia_osm) THEN 
     1114          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1115          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1116          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
     1117       END IF 
     1118 
     1119 
    12881120       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
    1289        CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. ) 
     1121       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 
    12901122 
    12911123       ! GN 25/8: need to change tmask --> wmask 
     
    13161148           END DO 
    13171149        END DO 
     1150        ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged) 
     1151        CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 
    13181152        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged) 
    13191153        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    13201154        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   & 
    1321          &                  ghamu, 'U', 1. , ghamv, 'V', 1. ) 
    1322  
    1323        IF(ln_dia_osm) THEN 
     1155         &                  ghamu, 'U', -1. , ghamv, 'V', -1. ) 
     1156 
     1157      IF(ln_dia_osm) THEN 
    13241158         SELECT CASE (nn_osm_wave) 
    13251159         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1). 
     
    13301164         ! Stokes drift read in from sbcwave  (=2). 
    13311165         CASE(2) 
    1332             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    1333             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
     1166            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
     1167            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     1168            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period 
     1169            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height 
     1170            IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) )                  ! wave mean period from NP spectrum 
     1171            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum 
     1172            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10 
    13341173            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
    13351174                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
     
    13421181         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0> 
    13431182         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth 
    1344          IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth 
     1183         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1184         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
     1185         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
     1186         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base 
     1187         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base 
     1188         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base 
     1189         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth 
     1190         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth 
    13451191         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth 
    13461192         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points 
     
    13481194         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    13491195         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
     1196         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale 
     1197         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir # 
    13501198         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    13511199         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    13521200         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    13531201         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    1354          IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine 
     1202         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1203         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    13551204         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    1356          IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine 
    1357          IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine 
    1358          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML 
     1205         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
     1206         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
     1207         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
     1208         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
     1209         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
     1210 
     1211         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )               ! FK layer depth 
     1212         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )               ! FK target layer depth 
     1213         IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux 
     1214         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML 
     1215         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 
     1216         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )            ! FK dtdx at u-pt 
     1217         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )            ! FK dtdy at v-pt 
     1218         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )            ! FK dtdx at u-pt 
     1219         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )            ! FK dsdy at v-pt 
     1220         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )            ! FK dbdx at u-pt 
     1221         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )            ! FK dbdy at v-pt 
     1222         IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1223         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1224 
    13591225      END IF 
    1360       ! Lateral boundary conditions on p_avt  (sign unchanged) 
    1361       CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 
     1226 
     1227CONTAINS 
     1228 
     1229 
     1230   ! Alan: do we need zb? 
     1231   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1232     !!--------------------------------------------------------------------- 
     1233     !!                   ***  ROUTINE zdf_vertical_average  *** 
     1234     !! 
     1235     !! ** Purpose : Determines vertical averages from surface to jnlev. 
     1236     !! 
     1237     !! ** Method  : Averages are calculated from the surface to jnlev. 
     1238     !!              The external level used to calculate differences is ibld+ibld_ext 
     1239     !! 
     1240     !!---------------------------------------------------------------------- 
     1241 
     1242        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1243 
     1244        ! Alan: do we need zb? 
     1245        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity 
     1246        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
     1247        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
     1248        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
     1249 
     1250        INTEGER :: jk, ji, jj 
     1251        REAL(wp) :: zthick, zthermal, zbeta 
     1252 
     1253 
     1254        zt   = 0._wp 
     1255        zs   = 0._wp 
     1256        zu   = 0._wp 
     1257        zv   = 0._wp 
     1258        DO jj = 2, jpjm1                                 !  Vertical slab 
     1259         DO ji = 2, jpim1 
     1260            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1261            zbeta    = rab_n(ji,jj,1,jp_sal) 
     1262               ! average over depth of boundary layer 
     1263            zthick = epsln 
     1264            DO jk = 2, jnlev_av(ji,jj) 
     1265               zthick = zthick + e3t_n(ji,jj,jk) 
     1266               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 
     1267               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     1268               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) & 
     1269                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 
     1270                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 
     1271               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) & 
     1272                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 
     1273                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 
     1274            END DO 
     1275            zt(ji,jj) = zt(ji,jj) / zthick 
     1276            zs(ji,jj) = zs(ji,jj) / zthick 
     1277            zu(ji,jj) = zu(ji,jj) / zthick 
     1278            zv(ji,jj) = zv(ji,jj) / zthick 
     1279            ! Alan: do we need zb? 
     1280            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 
     1281            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 
     1282            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 
     1283                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 
     1284            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 
     1285                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 
     1286            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1287         END DO 
     1288        END DO 
     1289   END SUBROUTINE zdf_osm_vertical_average 
     1290 
     1291   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 
     1292     !!--------------------------------------------------------------------- 
     1293     !!                   ***  ROUTINE zdf_velocity_rotation  *** 
     1294     !! 
     1295     !! ** Purpose : Rotates frame of reference of averaged velocity components. 
     1296     !! 
     1297     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w 
     1298     !! 
     1299     !!---------------------------------------------------------------------- 
     1300 
     1301        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle 
     1302        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current 
     1303        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline 
     1304 
     1305        INTEGER :: ji, jj 
     1306        REAL(wp) :: ztemp 
     1307 
     1308        DO jj = 2, jpjm1 
     1309           DO ji = 2, jpim1 
     1310              ztemp = zu(ji,jj) 
     1311              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 
     1312              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1313              ztemp = zdu(ji,jj) 
     1314              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
     1315              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1316           END DO 
     1317        END DO 
     1318    END SUBROUTINE zdf_osm_velocity_rotation 
     1319 
     1320    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 
     1321     !!--------------------------------------------------------------------- 
     1322     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     1323     !! 
     1324     !! ** Purpose : Calculates the gradients below the OSBL 
     1325     !! 
     1326     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1327     !! 
     1328     !!---------------------------------------------------------------------- 
     1329 
     1330     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
     1331 
     1332     INTEGER :: jj, ji, jkb, jkb1 
     1333     REAL(wp) :: zthermal, zbeta 
     1334 
     1335 
     1336     DO jj = 2, jpjm1 
     1337        DO ji = 2, jpim1 
     1338           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1339              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1340              zbeta    = rab_n(ji,jj,1,jp_sal) 
     1341              jkb = ibld(ji,jj) + ibld_ext 
     1342              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     1343              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 
     1344                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1345              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 
     1346                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1347              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1348           END IF 
     1349        END DO 
     1350     END DO 
     1351    END SUBROUTINE zdf_osm_external_gradients 
     1352 
     1353    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 
     1354 
     1355     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1356 
     1357     INTEGER :: jk, jj, ji 
     1358     REAL(wp) :: ztgrad, zsgrad, zbgrad 
     1359     REAL(wp) :: zgamma_b_nd, zgamma_c, znd 
     1360     REAL(wp) :: zzeta_s=0.3, ztmp 
     1361 
     1362     DO jj = 2, jpjm1 
     1363        DO ji = 2, jpim1 
     1364           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1365              IF ( lconv(ji,jj) ) THEN  ! convective conditions 
     1366                    IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 
     1367                         & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 
     1368                         &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0 
     1369                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
     1370                       ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 
     1371                       zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 
     1372                       zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 
     1373                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
     1374                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 
     1375                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 
     1376                       DO jk = 2, ibld(ji,jj)+ibld_ext 
     1377                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 
     1378                          IF ( znd <= zzeta_s ) THEN 
     1379                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 
     1380                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1381                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 
     1382                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1383                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 
     1384                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1385                          ELSE 
     1386                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1387                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1388                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1389                          ENDIF 
     1390                       END DO 
     1391                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 
     1392              ELSE 
     1393                 ! stable conditions 
     1394                 ! if pycnocline profile only defined when depth steady of increasing. 
     1395                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     1396                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1397                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     1398                          ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 
     1399                          ztgrad = zdt_bl(ji,jj) * ztmp 
     1400                          zsgrad = zds_bl(ji,jj) * ztmp 
     1401                          zbgrad = zdb_bl(ji,jj) * ztmp 
     1402                          DO jk = 2, ibld(ji,jj) 
     1403                             znd = gdepw_n(ji,jj,jk) * ztmp 
     1404                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1405                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1406                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1407                          END DO 
     1408                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1409                          ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
     1410                          ztgrad = zdt_bl(ji,jj) * ztmp 
     1411                          zsgrad = zds_bl(ji,jj) * ztmp 
     1412                          zbgrad = zdb_bl(ji,jj) * ztmp 
     1413                          DO jk = 2, ibld(ji,jj) 
     1414                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp 
     1415                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1416                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1417                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1418                          END DO 
     1419                       ENDIF ! IF (zhol >=0.5) 
     1420                    ENDIF    ! IF (zdb_bl> 0.) 
     1421                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1422              ENDIF          ! IF (lconv) 
     1423           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1424        END DO 
     1425     END DO 
     1426 
     1427    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
     1428 
     1429    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
     1430      !!--------------------------------------------------------------------- 
     1431      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
     1432      !! 
     1433      !! ** Purpose : Calculates velocity shear in the pycnocline 
     1434      !! 
     1435      !! ** Method  : 
     1436      !! 
     1437      !!---------------------------------------------------------------------- 
     1438 
     1439      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
     1440 
     1441      INTEGER :: jk, jj, ji 
     1442      REAL(wp) :: zugrad, zvgrad, znd 
     1443      REAL(wp) :: zzeta_v = 0.45 
    13621444      ! 
    1363    END SUBROUTINE zdf_osm 
     1445      DO jj = 2, jpjm1 
     1446         DO ji = 2, jpim1 
     1447            ! 
     1448            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1449               IF ( lconv (ji,jj) ) THEN 
     1450                  ! Unstable conditions 
     1451                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1452                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1453                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1454                  !Alan is this right? 
     1455                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1456                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1457                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1458                       &      )/ (zdh(ji,jj)  + epsln ) 
     1459                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     1460                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     1461                     IF ( znd <= 0.0 ) THEN 
     1462                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     1463                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     1464                     ELSE 
     1465                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     1466                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     1467                     ENDIF 
     1468                  END DO 
     1469               ELSE 
     1470                  ! stable conditions 
     1471                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     1472                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     1473                  DO jk = 2, ibld(ji,jj) 
     1474                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1475                     IF ( znd < 1.0 ) THEN 
     1476                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     1477                     ELSE 
     1478                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     1479                     ENDIF 
     1480                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     1481                  END DO 
     1482               ENDIF 
     1483               ! 
     1484            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1485         END DO 
     1486      END DO 
     1487    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1488 
     1489    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
     1490     !!--------------------------------------------------------------------- 
     1491     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     1492     !! 
     1493     !! ** Purpose : Calculates the rate at which hbl changes. 
     1494     !! 
     1495     !! ** Method  : 
     1496     !! 
     1497     !!---------------------------------------------------------------------- 
     1498 
     1499    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl 
     1500 
     1501    INTEGER :: jj, ji 
     1502    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 
     1503    REAL(wp) :: zvel_max, zwb_min 
     1504    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1505    REAL(wp) :: zzeta_m = 0.3 
     1506    REAL(wp) :: zgamma_c = 2.0 
     1507    REAL(wp) :: zdhoh = 0.1 
     1508    REAL(wp) :: alpha_bc = 0.5 
     1509 
     1510    DO jj = 2, jpjm1 
     1511       DO ji = 2, jpim1 
     1512          IF ( lconv(ji,jj) ) THEN    ! Convective 
     1513             ! Alan is this right?  Yes, it's a bit different from the previous relationship 
     1514             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1515             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1516             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1517             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1518             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1519             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) ) ) 
     1522 
     1523             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1524                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1525                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1526                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1527             ! 
     1528             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 
     1529 
     1530             IF ( ln_osm_mle ) THEN 
     1531                  !  zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 
     1532                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL 
     1533                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1534                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     1535                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     1536                ELSE 
     1537                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     1538                ENDIF 
     1539                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1540                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     1541                   ! OSBL is deepening, entrainment > restratification 
     1542                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1543                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1544                   ELSE 
     1545                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     1546                   ENDIF 
     1547                ELSE 
     1548                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     1549                   zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     1550 
     1551 
     1552                ENDIF 
     1553 
     1554             ELSE 
     1555                ! Fox-Kemper not used. 
     1556 
     1557                  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 
     1559                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1560                ! added ajgn 23 July as temporay fix 
     1561 
     1562             ENDIF 
     1563 
     1564             zdhdt_2(ji,jj) = 0._wp 
     1565 
     1566                ! commented out ajgn 23 July as temporay fix 
     1567!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1568! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 
     1569!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1570!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 
     1571!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1572!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 
     1573!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 
     1574!                     ! Alan no idea what this should be? 
     1575!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 
     1576!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 
     1577!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 
     1578!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 
     1579!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 
     1580!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 
     1581!                 ENDIF 
     1582            ELSE                        ! Stable 
     1583                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     1584                zdhdt_2(ji,jj) = 0._wp 
     1585                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     1586                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     1587                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     1588                ELSE 
     1589                    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) ) 
     1590                ENDIF 
     1591                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 
     1592            ENDIF 
     1593         END DO 
     1594      END DO 
     1595    END SUBROUTINE zdf_osm_calculate_dhdt 
     1596 
     1597    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     1598     !!--------------------------------------------------------------------- 
     1599     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
     1600     !! 
     1601     !! ** Purpose : Increments hbl. 
     1602     !! 
     1603     !! ** Method  : If thechange in hbl exceeds one model level the change is 
     1604     !!              is calculated by moving down the grid, changing the buoyancy 
     1605     !!              jump. This is to ensure that the change in hbl does not 
     1606     !!              overshoot a stable layer. 
     1607     !! 
     1608     !!---------------------------------------------------------------------- 
     1609 
     1610 
     1611    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl. 
     1612 
     1613    INTEGER :: jk, jj, ji, jm 
     1614    REAL(wp) :: zhbl_s, zvel_max, zdb 
     1615    REAL(wp) :: zthermal, zbeta 
     1616 
     1617     DO jj = 2, jpjm1 
     1618         DO ji = 2, jpim1 
     1619           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     1620! 
     1621! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     1622! 
     1623              zhbl_s = hbl(ji,jj) 
     1624              jm = imld(ji,jj) 
     1625              zthermal = rab_n(ji,jj,1,jp_tem) 
     1626              zbeta = rab_n(ji,jj,1,jp_sal) 
     1627 
     1628 
     1629              IF ( lconv(ji,jj) ) THEN 
     1630!unstable 
     1631 
     1632                 IF( ln_osm_mle ) THEN 
     1633                    zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1634                 ELSE 
     1635 
     1636                    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) / & 
     1637                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1638 
     1639                 ENDIF 
     1640 
     1641                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1642                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
     1643                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 
     1644                         &  0.0 ) + zvel_max 
     1645 
     1646 
     1647                    IF ( ln_osm_mle ) THEN 
     1648                       zhbl_s = zhbl_s + MIN( & 
     1649                        & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1650                        & e3w_n(ji,jj,jm) ) 
     1651                    ELSE 
     1652                      zhbl_s = zhbl_s + MIN( & 
     1653                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1654                        & e3w_n(ji,jj,jm) ) 
     1655                    ENDIF 
     1656 
     1657                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1658 
     1659                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
     1660                 END DO 
     1661                 hbl(ji,jj) = zhbl_s 
     1662                 ibld(ji,jj) = jm 
     1663             ELSE 
     1664! stable 
     1665                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1666                    zdb = MAX( & 
     1667                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 
     1668                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 
     1669                         & 0.0 ) + & 
     1670             &       2.0 * zvstr(ji,jj)**2 / zhbl_s 
     1671 
     1672                    ! Alan is thuis right? I have simply changed hbli to hbl 
     1673                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 
     1674                    zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & 
     1675               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 
     1676                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 
     1677                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 
     1678 
     1679                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1680                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
     1681                 END DO 
     1682             ENDIF   ! IF ( lconv ) 
     1683             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 
     1684             ibld(ji,jj) = MAX(jm, 4 ) 
     1685           ELSE 
     1686! change zero or one model level. 
     1687             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 
     1688           ENDIF 
     1689           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     1690         END DO 
     1691      END DO 
     1692 
     1693    END SUBROUTINE zdf_osm_timestep_hbl 
     1694 
     1695    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     1696      !!--------------------------------------------------------------------- 
     1697      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     1698      !! 
     1699      !! ** Purpose : Calculates thickness of the pycnocline 
     1700      !! 
     1701      !! ** Method  : The thickness is calculated from a prognostic equation 
     1702      !!              that relaxes the pycnocine thickness to a diagnostic 
     1703      !!              value. The time change is calculated assuming the 
     1704      !!              thickness relaxes exponentially. This is done to deal 
     1705      !!              with large timesteps. 
     1706      !! 
     1707      !!---------------------------------------------------------------------- 
     1708 
     1709      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
     1710      ! 
     1711      INTEGER :: jj, ji 
     1712      INTEGER :: inhml 
     1713      REAL(wp) :: zari, ztau, zddhdt 
     1714 
     1715 
     1716      DO jj = 2, jpjm1 
     1717         DO ji = 2, jpim1 
     1718            IF ( lconv(ji,jj) ) THEN 
     1719 
     1720               IF( ln_osm_mle ) THEN 
     1721                  IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
     1722                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
     1723                     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 
     1725                           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 ) 
     1727                        ELSE                                                     ! unstable 
     1728                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1729                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1730                        ENDIF 
     1731                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1732                        zddhdt = zari * hbl(ji,jj) 
     1733                     ELSE 
     1734                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1735                        zddhdt = 0.2 * hbl(ji,jj) 
     1736                     ENDIF 
     1737                  ELSE 
     1738                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1739                     zddhdt = 0.2 * hbl(ji,jj) 
     1740                  ENDIF 
     1741               ELSE ! ln_osm_mle 
     1742                  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 ) 
     1746                     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 ) 
     1749                     ENDIF 
     1750                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1751                     zddhdt = zari * hbl(ji,jj) 
     1752                  ELSE 
     1753                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1754                     zddhdt = 0.2 * hbl(ji,jj) 
     1755                  ENDIF 
     1756 
     1757               END IF  ! ln_osm_mle 
     1758 
     1759               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1760               ! Alan: this hml is never defined or used 
     1761            ELSE   ! IF (lconv) 
     1762               ztau = hbl(ji,jj) / zvstr(ji,jj) 
     1763               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     1764                  ! boundary layer deepening 
     1765                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1766                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     1767                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     1768                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     1769                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
     1770                  ELSE 
     1771                     zddhdt = 0.2 * hbl(ji,jj) 
     1772                  ENDIF 
     1773               ELSE     ! IF(dhdt < 0) 
     1774                  zddhdt = 0.2 * hbl(ji,jj) 
     1775               ENDIF    ! IF (dhdt >= 0) 
     1776               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1777               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse 
     1778               ! Alan: this hml is never defined or used -- do we need it? 
     1779            ENDIF       ! IF (lconv) 
     1780 
     1781            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     1782            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     1783            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     1784            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     1785            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1786         END DO 
     1787      END DO 
     1788 
     1789    END SUBROUTINE zdf_osm_pycnocline_thickness 
     1790 
     1791 
     1792   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     1793      !!---------------------------------------------------------------------- 
     1794      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  *** 
     1795      !! 
     1796      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 
     1797      !! 
     1798      !! ** Method  : 
     1799      !! 
     1800      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1801      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1802 
     1803 
     1804      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 
     1805      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
     1806      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy 
     1807 
     1808      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1809      INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     1810      REAL(wp)                         :: zc 
     1811      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value 
     1812      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH 
     1813      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 
     1814      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv 
     1815!!---------------------------------------------------------------------- 
     1816      ! 
     1817      !                                      !==  MLD used for MLE  ==! 
     1818 
     1819      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point 
     1820      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
     1821      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     1822      DO jk = nlb10, jpkm1 
     1823         DO jj = 1, jpj                ! Mixed layer level: w-level  
     1824            DO ji = 1, jpi 
     1825               ikt = mbkt(ji,jj) 
     1826               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 
     1827               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     1828            END DO 
     1829         END DO 
     1830      END DO 
     1831      DO jj = 1, jpj 
     1832         DO ji = 1, jpi 
     1833            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 
     1834            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1835         END DO 
     1836      END DO 
     1837      ! ensure mld_prof .ge. ibld 
     1838      ! 
     1839      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     1840      ! 
     1841      ztm(:,:) = 0._wp 
     1842      zsm(:,:) = 0._wp 
     1843      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
     1844         DO jj = 1, jpj 
     1845            DO ji = 1, jpi 
     1846               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     1847               ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 
     1848               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 
     1849            END DO 
     1850         END DO 
     1851      END DO 
     1852      ! average temperature and salinity. 
     1853      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1854      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1855      ! calculate horizontal gradients at u & v points 
     1856 
     1857      DO jj = 2, jpjm1 
     1858         DO ji = 1, jpim1 
     1859            zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1860            zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1861            zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 
     1862            ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 
     1863            ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 
     1864         END DO 
     1865      END DO 
     1866 
     1867      DO jj = 1, jpjm1 
     1868         DO ji = 2, jpim1 
     1869            zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1870            zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1871            zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 
     1872            ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 
     1873            ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 
     1874         END DO 
     1875      END DO 
     1876 
     1877      CALL eos_rab(ztsm_midu, zmld_midu, zabu) 
     1878      CALL eos_rab(ztsm_midv, zmld_midv, zabv) 
     1879 
     1880      DO jj = 2, jpjm1 
     1881         DO ji = 1, jpim1 
     1882            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 
     1883         END DO 
     1884      END DO 
     1885      DO jj = 1, jpjm1 
     1886         DO ji = 2, jpim1 
     1887            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 
     1888         END DO 
     1889      END DO 
     1890 
     1891 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
     1892  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     1893      !!---------------------------------------------------------------------- 
     1894      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
     1895      !! 
     1896      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 
     1897      !! 
     1898      !! ** Method  : 
     1899      !! 
     1900      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1901      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1902 
     1903      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle 
     1904      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1905      INTEGER  ::   ii, ij, ik  ! local integers 
     1906      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     1907      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle 
     1908 
     1909      mld_prof(:,:) = 4 
     1910      DO jk = 5, jpkm1 
     1911        DO jj = 2, jpjm1 
     1912          DO ji = 2, jpim1 
     1913            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1914          END DO 
     1915        END DO 
     1916      END DO 
     1917      ! DO jj = 2, jpjm1 
     1918      !    DO ji = 1, jpim1 
     1919      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1920      !    END DO 
     1921      !  END DO 
     1922   ! Timestep mixed layer eddy depth. 
     1923      DO jj = 2, jpjm1 
     1924        DO ji = 2, jpim1 
     1925          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 
     1926          IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 
     1927             hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh )  ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 
     1928          ELSE 
     1929            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 
     1930             ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 
     1931             !   hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     1932             ! ELSE 
     1933             !   hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
     1934             ! ENDIF 
     1935             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1936               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     1937             ELSE 
     1938               hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
     1939             ENDIF 
     1940          ENDIF 
     1941          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
     1942          hmle(ji,jj) = MIN(hmle(ji,jj), 1.2*zmld(ji,jj))  
     1943        END DO 
     1944      END DO 
     1945 
     1946      mld_prof = 4 
     1947      DO jk = 5, jpkm1 
     1948        DO jj = 2, jpjm1 
     1949          DO ji = 2, jpim1 
     1950            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1951          END DO 
     1952        END DO 
     1953      END DO 
     1954      DO jj = 2, jpjm1 
     1955         DO ji = 2, jpim1 
     1956            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 
     1957         END DO 
     1958       END DO 
     1959   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
     1960 
     1961      DO jj = 2, jpjm1 
     1962        DO ji = 2, jpim1 
     1963          IF ( lconv(ji,jj) ) THEN 
     1964             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     1965             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 
     1966             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 
     1967             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 
     1968             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 
     1969             zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
     1970                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
     1971             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)  
     1974             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf 
     1975          ENDIF 
     1976        END DO 
     1977      END DO 
     1978END SUBROUTINE zdf_osm_mle_parameters 
     1979 
     1980END SUBROUTINE zdf_osm 
    13641981 
    13651982 
     
    13781995     INTEGER  ::   ios            ! local integer 
    13791996     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     1997     REAL z1_t2 
    13801998     !! 
    13811999     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    13822000          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
    1383           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 
     2001          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
     2002          & ,ln_osm_mle 
     2003! Namelist for Fox-Kemper parametrization. 
     2004      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 
     2006 
    13842007     !!---------------------------------------------------------------------- 
    13852008     ! 
     
    13972020        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 
    13982021        WRITE(numout,*) '~~~~~~~~~~~~' 
    1399         WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters' 
    1400         WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la 
     2022        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters' 
     2023        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la 
     2024        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    14012025        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
    14022026        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    1403         WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes 
     2027        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
    14042028        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave 
    14052029        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave 
     
    14232047     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    14242048 
    1425      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     2049 
     2050     IF( ln_osm_mle ) THEN 
     2051! Initialise Fox-Kemper parametrization 
     2052         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme 
     2053         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 
     2054903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 
     2055 
     2056         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 
     2057         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 
     2058904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 
     2059         IF(lwm) WRITE ( numond, namosm_mle ) 
     2060 
     2061         IF(lwp) THEN                     ! Namelist print 
     2062            WRITE(numout,*) 
     2063            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 
     2064            WRITE(numout,*) '~~~~~~~~~~~~~' 
     2065            WRITE(numout,*) '   Namelist namosm_mle : ' 
     2066            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle 
     2067            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce 
     2068            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm' 
     2069            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's' 
     2070            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg' 
     2071            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c 
     2072            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
     2073            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
     2074         ENDIF         ! 
     2075     ENDIF 
     2076      ! 
     2077      IF(lwp) THEN 
     2078         WRITE(numout,*) 
     2079         IF( ln_osm_mle ) THEN 
     2080            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 
     2081            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation' 
     2082            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation' 
     2083         ELSE 
     2084            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 
     2085         ENDIF 
     2086      ENDIF 
     2087      ! 
     2088      IF( ln_osm_mle ) THEN                ! MLE initialisation 
     2089         ! 
     2090         rb_c = grav * rn_osm_mle_rho_c /rau0        ! Mixed Layer buoyancy criteria 
     2091         IF(lwp) WRITE(numout,*) 
     2092         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
     2093         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 
     2094         ! 
     2095         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            ! 
     2096! 
     2097         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation 
     2098            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  ) 
     2099            ! 
     2100         ENDIF 
     2101         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
     2102         z1_t2 = 2.e-5 
     2103         do jj=1,jpj 
     2104            do ji = 1,jpi 
     2105               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 
     2106            end do 
     2107         end do 
     2108         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 
     2109         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  ) 
     2110         ! 
     2111      ENDIF 
     2112 
     2113     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 
     2114 
    14262115 
    14272116     IF( ln_zdfddm) THEN 
     
    15122201        CALL iom_set_rstw_var_active('wn') 
    15132202        CALL iom_set_rstw_var_active('hbl') 
    1514         CALL iom_set_rstw_var_active('hbli') 
     2203        CALL iom_set_rstw_var_active('dh') 
     2204        IF( ln_osm_mle ) THEN 
     2205            CALL iom_set_rstw_var_active('hmle') 
     2206        END IF 
    15152207     ENDIF 
    15162208   END SUBROUTINE zdf_osm_init 
     
    15302222     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    15312223 
    1532      INTEGER ::   id1, id2   ! iom enquiry index 
     2224     INTEGER ::   id1, id2, id3   ! iom enquiry index 
    15332225     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    15342226     INTEGER  ::   iiki, ikt ! local integer 
     
    15362228     REAL(wp) ::   zN2_c           ! local scalar 
    15372229     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    1538      INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top) 
     2230     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 
    15392231     !!---------------------------------------------------------------------- 
    15402232     ! 
     
    15512243           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    15522244        END IF 
     2245 
    15532246        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    1554         id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. ) 
     2247        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. ) 
    15552248        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return 
    15562249           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 
    1557            CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios  ) 
    1558            WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file' 
     2250           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  ) 
     2251           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file' 
     2252           IF( ln_osm_mle ) THEN 
     2253              id3 = iom_varid( numror, 'hmle'   , ldstop = .FALSE. ) 
     2254              IF( id3 > 0) THEN 
     2255                 CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios ) 
     2256                 WRITE(numout,*) ' ===>>>> :  hmle read from restart file' 
     2257              ELSE 
     2258                 WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl' 
     2259                 hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2260              END IF 
     2261           END IF 
    15592262           RETURN 
    1560         ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate 
     2263        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate 
    15612264           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 
    15622265        END IF 
     
    15682271     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return 
    15692272        IF(lwp) WRITE(numout,*) '---- osm-rst ----' 
    1570          CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
    1571          CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    1572          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
     2273         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn,   ldxios = lwxios ) 
     2274         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl,  ldxios = lwxios ) 
     2275         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh,   ldxios = lwxios ) 
     2276         IF( ln_osm_mle ) THEN 
     2277            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios ) 
     2278         END IF 
    15732279        RETURN 
    15742280     END IF 
     
    15782284     !!----------------------------------------------------------------------------- 
    15792285     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 
    1580      ALLOCATE( imld_rst(jpi,jpj) ) 
    15812286     ! w-level of the mixing and mixed layers 
    15822287     CALL eos_rab( tsn, rab_n ) 
     
    15992304     DO jj = 1, jpj 
    16002305        DO ji = 1, jpi 
    1601            iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
     2306           iiki = MAX(4,imld_rst(ji,jj)) 
     2307           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth 
     2308           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth 
    16032309        END DO 
    16042310     END DO 
    1605      hbl = MAX(hbl,epsln) 
    1606      hbli(:,:) = hbl(:,:) 
    1607      DEALLOCATE( imld_rst ) 
     2311 
    16082312     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
     2313 
     2314     IF( ln_osm_mle ) THEN 
     2315        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2316        WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 
     2317     END IF 
     2318 
     2319     wn(:,:,:) = 0._wp 
     2320     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    16092321   END SUBROUTINE osm_rst 
    16102322 
     
    16342346      ENDIF 
    16352347 
    1636       ! add non-local temperature and salinity flux 
    16372348      DO jk = 1, jpkm1 
    16382349         DO jj = 2, jpjm1 
     
    16482359      END DO 
    16492360 
    1650  
    1651       ! save the non-local tracer flux trends for diagnostic 
     2361      ! save the non-local tracer flux trends for diagnostics 
    16522362      IF( l_trdtra )   THEN 
    16532363         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    16542364         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    1655 !!bug gm jpttdzdf ==> jpttosm 
    1656          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    1657          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     2365 
     2366         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 
     2367         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 
    16582368         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    16592369      ENDIF 
     
    17232433 
    17242434   !!====================================================================== 
     2435 
    17252436END MODULE zdfosm 
  • NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF/zdfphy.F90

    r12658 r12912  
    5656   !!---------------------------------------------------------------------- 
    5757   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    58    !! $Id$ 
     58   !! $Id: zdfphy.F90 12178 2019-12-11 11:02:38Z agn $ 
    5959   !! Software governed by the CeCILL license (see ./LICENSE) 
    6060   !!---------------------------------------------------------------------- 
     
    172172      IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' ) 
    173173      IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 
    174       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' ) 
     174      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 
    175175      IF(lwp) THEN 
    176176         WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.