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

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

Changeset 14051


Ignore:
Timestamp:
2020-12-03T14:28:02+01:00 (4 years ago)
Author:
techene
Message:

#2385 branch updated with trunk 14045 (pb in ORCA2_ICE_OBS to fix later)

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
Files:
1 deleted
12 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/cfgs/SHARED/field_def_nemo-oce.xml

    r14050 r14051  
    271271 
    272272      <field_group id="OSMOSIS_T" grid_ref="grid_T_2D"> 
     273        <field id="hml"                 long_name="mixed layr depth"                         unit="m"       /> 
     274        <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
     275        <field id="dh"                  long_name="Pycnocline thickness"                     unit=" m"      /> 
     276        <field id="ibld"                long_name="index of boundary layer depth"            unit="#"       /> 
     277        <field id="imld"                long_name="index of mixed layer depth"            unit="#"       /> 
     278        <field id="zhbl"                long_name="boundary layer depth -grid"                     unit="m"       /> 
     279        <field id="zhml"                long_name="mixed layer depth - grid"                        unit="m"       /> 
     280        <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
     281        <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
     282        <field id="us_x"        long_name="i component of active Stokes drift"                      unit="m/s"     /> 
     283        <field id="us_y"        long_name="j component of active Stokes drift"                      unit="m/s"     /> 
     284        <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    273285        <field id="zwth0"               long_name="surface non-local temperature flux"       unit="deg m/s" /> 
    274286        <field id="zws0"                long_name="surface non-local salinity flux"          unit="psu m/s" /> 
    275         <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
    276         <field id="hbli"                long_name="initial boundary layer depth"             unit="m"       /> 
    277         <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    278         <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
    279287        <field id="zwstrc"              long_name="convective velocity scale"                unit="m/s"     /> 
     288        <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    280289        <field id="zwstrl"              long_name="langmuir velocity scale"                  unit="m/s"     /> 
    281         <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    282         <field id="zhbl"                long_name="boundary layer depth"                     unit="m"       /> 
    283         <field id="zhml"                long_name="mixed layer depth"                        unit="m"       /> 
     290        <field id="zvstr"               long_name="mixed velocity scale"                     unit="m/s"     /> 
     291        <field id="zla"                 long_name="langmuir number"                          unit="m/s"     /> 
    284292        <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2"                        unit="mW"      /> 
    285293        <field id="wind_wave_power"     long_name="U_s \dot  tau"                            unit="mW"      /> 
    286294        <field id="wind_power"          long_name="\rho  u*^3"                               unit="mW"      /> 
    287295 
    288         <!-- extra OSMOSIS diagnostics --> 
     296       <!-- interior BL OSMOSIS diagnostics --> 
    289297        <field id="zwthav"              long_name="av turb flux of T in ml"                  unit="deg m/s" /> 
    290298        <field id="zt_ml"               long_name="av T in ml"                               unit="deg"     /> 
     299        <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
     300        <field id="zws_ent"            long_name="entrainment turb flux of S"                unit="10^-3 m/s" /> 
    291301        <field id="zwth_ent"            long_name="entrainment turb flux of T"               unit="deg m/s" /> 
    292         <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
    293         <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
    294       </field_group> 
    295  
    296       <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" operation="instant" > 
     302        <field id="zwb_ent"            long_name="entrainment turb flux of buoyancy"         unit="m^2/s^-3" /> 
     303  
     304        <field id="zdt_bl"             long_name="temperature jump at base of BL"                 unit="deg"      /> 
     305        <field id="zds_bl"             long_name="salinity jump at base of BL"                 unit="10^-3"      /> 
     306        <field id="zdb_bl"             long_name="buoyancy jump at base of BL"                 unit="m/s^2"      /> 
     307        <field id="zdu_bl"             long_name="u jump at base of BL"                       unit="m/s"      /> 
     308        <field id="zdv_bl"             long_name="v jump at base of BL"                       unit="m/s"      /> 
     309 
     310        <!-- extra OSMOSIS diagnostics for debugging --> 
     311       <field id="zsc_uw_1_0"       long_name="zsc u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 
     312        <field id="zsc_uw_1_f"       long_name="zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     313        <field id="zsc_vw_1_f"       long_name="zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     314        <field id="zsc_uw_2_f"       long_name="2nd zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     315        <field id="zsc_vw_2_f"       long_name="2nd zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     316        <field id="zuw_bse"       long_name="base u-flux T-points"                          unit="m^2/s^2" /> 
     317        <field id="zvw_bse"       long_name="base v-flux T-points"                          unit="m^2/s^2" /> 
     318 
     319       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     320         <field id="hmle"          long_name="OBL FK-layer thickness"                                     unit="m"        /> 
     321        <field id="mld_prof"              long_name="FK-layer depth index"                  unit="#" /> 
     322        <field id="zmld"          long_name="target FK-layer thickness"                                     unit="m"        /> 
     323        <field id="zwb_fk"          long_name="FK b-flux"                                     unit="m^2 s^-3"        /> 
     324        <field id="zwb_fk_b"          long_name="layer averaged FK b-flux"                 unit="m^2 s^-3"       /> 
     325        <field id="zdiff_mle"          long_name="max FK diffusivity in MLE"       unit=" 10^-4 m^2 s^-1"       /> 
     326        <field id="zvel_mle"          long_name="FK velocity scale in MLE"       unit=" m s^-1"       /> 
     327    </field_group> 
     328 
     329      <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" > 
     330        <field id="zviscos"       long_name="BL viscosity"   unit="m^2/s" /> 
    297331        <field id="ghamt"       long_name="non-local temperature flux"                       unit="deg m/s" /> 
    298332        <field id="ghams"       long_name="non-local salinity flux"                          unit="psu m/s" /> 
    299333        <field id="zdtdz_pyc"   long_name="Pycnocline temperature gradient"                  unit=" deg/m"  /> 
    300       </field_group> 
     334        <field id="zdsdz_pyc"   long_name="Pycnocline salinity gradient"                  unit=" 10^-3/m"  /> 
     335        <field id="zdbdz_pyc"   long_name="Pycnocline buoyancy gradient"                  unit=" s^-2"  /> 
     336        <field id="zdudz_pyc"   long_name="Pycnocline u gradient"                  unit=" s^-2"  /> 
     337        <field id="zdvdz_pyc"   long_name="Pycnocline v gradient"                  unit=" s^-2"  /> 
     338 
     339        <!-- extra OSMOSIS diagnostics for debugging --> 
     340         <field id="ghamu_00"       long_name="initial non-local u-momentum flux"   unit="m^2/s^2" /> 
     341        <field id="ghamv_00"       long_name="initial non-local v-momentum flux"   unit="m^2/s^2" /> 
     342        <field id="ghamu_0"       long_name="after dstokes non-local u-momentum flux"   unit="m^2/s^2" /> 
     343        <field id="ghamu_f"       long_name="after Coriolis non-local u-momentum flux"   unit="m^2/s^2" /> 
     344        <field id="ghamv_f"       long_name="after Coriolis  non-local v-momentum flux"   unit="m^2/s^2" /> 
     345        <field id="ghamu_b"       long_name="after buoyancy added non-local u-momentum flux"   unit="m^2/s^2" /> 
     346        <field id="ghamv_b"       long_name="after buoyancy added  non-local v-momentum flux"  unit="m^2/s^2" /> 
     347        <field id="ghamu_1"       long_name="after entrainment non-local u-momentum flux"   unit="m^2/s^2" /> 
     348        <field id="ghamv_1"       long_name="after entrainment  non-local v-momentum flux"  unit="m^2/s^2" /> 
     349     </field_group> 
    301350 
    302351      <field_group id="OSMOSIS_U" grid_ref="grid_U_2D" > 
    303352        <field id="ghamu"       long_name="non-local u-momentum flux"   grid_ref="grid_U_3D" unit="m^2/s^2" /> 
    304         <field id="us_x"        long_name="i component of Stokes drift"                      unit="m/s"     /> 
    305       </field_group> 
     353       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     354       <field id="zdtdx"          long_name="FK  T x-gradient"                                     unit=" deg C m^-1"        /> 
     355        <field id="zdsdx"          long_name="FK  S x-gradient"                                     unit=" 10^-3 m^-1"        /> 
     356        <field id="dbdx_mle"          long_name="FK  B x-gradient"                                     unit=" s^-2"        /> 
     357     </field_group> 
    306358 
    307359      <field_group id="OSMOSIS_V" grid_ref="grid_V_2D" > 
    308360        <field id="ghamv"       long_name="non-local v-momentum flux"   grid_ref="grid_V_3D" unit="m^2/s^2" /> 
    309         <field id="us_y"        long_name="j component of Stokes drift"                      unit="m/s"     /> 
     361        <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     362        <field id="zdtdy"          long_name="FK T y-gradient"                                     unit=" deg C m^-1"        /> 
     363        <field id="zdsdy"          long_name="FK S y-gradient"                                     unit=" 10^-3 m^-1"        /> 
     364        <field id="dbdy_mle"          long_name="FK B y-gradient"                                     unit=" s^-2"        /> 
    310365      </field_group> 
    311366 
    … …  
    856911     <field id="strd_zdfp"     long_name="salinity   -trend: pure vert. diffusion"   unit="1e-3/s" /> 
    857912 
    858      <!-- --> 
     913     <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     914     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s" /> 
     915     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s" /> 
     916 
     917 
     918    <!-- --> 
    859919     <field id="ttrd_dmp"      long_name="temperature-trend: interior restoring"        unit="degC/s" /> 
    860920     <field id="strd_dmp"      long_name="salinity   -trend: interior restoring"        unit="1e-3/s" /> 
    … …  
    892952     <field id="strd_zdfp_e3t"     unit="1e-3/s * m"  >  strd_zdfp * e3t </field> 
    893953 
     954          <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     955     <field id="ttrd_osm_e3t"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field> 
     956     <field id="strd_osm_e3t"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field> 
     957      
    894958     <!-- --> 
    895959     <field id="ttrd_dmp_e3t"      unit="degC/s * m"  >  ttrd_dmp * e3t </field> 
    … …  
    907971     <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> 
    908972     <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> 
     973     <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> 
     974     <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> 
    909975     <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> 
    910976     <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> 
    … …  
    11141180    </field_group> 
    11151181 
     1182    <!-- TMB diagnostic output --> 
     1183    <field_group  id="1h_grid_T_tmb" grid_ref="grid_T_2D" operation="instant"> 
     1184      <field id="top_temp"           name="votemper_top"  unit="degC"  /> 
     1185      <field id="mid_temp"           name="votemper_mid"  unit="degC"  /> 
     1186      <field id="bot_temp"           name="votemper_bot"  unit="degC"  /> 
     1187      <field id="top_sal"            name="vosaline_top"  unit="psu"   /> 
     1188      <field id="mid_sal"            name="vosaline_mid"  unit="psu"   /> 
     1189      <field id="bot_sal"            name="vosaline_bot"  unit="psu"   /> 
     1190      <field id="sshnmasked"         name="sossheig"      unit="m"     />  
     1191    </field_group> 
     1192 
    11161193    <field_group  id="1h_grid_U_tmb" grid_ref="grid_U_2D" operation="instant"> 
    11171194      <field id="top_u"           name="vozocrtx_top"  unit="m/s"  /> 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/cfgs/SHARED/namelist_ref

    r14050 r14051  
    12191219&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
    12201220!----------------------------------------------------------------------- 
    1221    ln_use_osm_la = .false.      !  Use namelist  rn_osm_la 
     1221   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    12221222   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    1223    rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     1223   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor 
     1224   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4 
     1225   rn_osm_bl_thresh   = 5.e-5      !Threshold buoyancy for deepening of OSBL base 
    12241226   nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv 
    12251227   ln_dia_osm = .true.         ! output OSMOSIS-OBL variables 
    … …  
    12291231   rn_difri  =  0.005          ! max Ri# diffusivity at Ri_g = 0 (m^2/s) 
    12301232   ln_convmix  = .true.        ! Use convective instability mixing below BL 
    1231    rn_difconv = 1.             ! diffusivity when unstable below BL  (m2/s) 
     1233   rn_difconv = 1. !0.01 !1.             ! diffusivity when unstable below BL  (m2/s) 
     1234   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
    12321235   nn_osm_wave = 0             ! Method used to calculate Stokes drift 
    12331236      !                        !  = 2: Use ECMWF wave fields 
    12341237      !                        !  = 1: Pierson Moskowitz wave spectrum 
    12351238      !                        !  = 0: Constant La# = 0.3 
    1236 / 
     1239   nn_osm_SD_reduce = 0        ! Method used to get active Stokes drift from surface value 
     1240      !                        !  = 0: No reduction 
     1241                               !  = 1: use SD avged over top 10% hbl 
     1242                               !  = 2:use surface value of SD fit to slope at rn_osm_hblfrac*hbl below surface 
     1243   ln_zdfosm_ice_shelter = .true.  ! reduce surface SD and depth scale under ice 
     1244   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
     1245/ 
     1246!----------------------------------------------------------------------- 
     1247&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     1248!----------------------------------------------------------------------- 
     1249   rn_osm_mle_ce       = 0.06      ! magnitude of the MLE (typical value: 0.06 to 0.08) 
     1250   nn_osm_mle          = 0         ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 
     1251   rn_osm_mle_lf       = 5.e+3     ! typical scale of mixed layer front (meters)                      (case rn_osm_mle=0) 
     1252   rn_osm_mle_time     = 172800.   ! time scale for mixing momentum across the mixed layer (seconds)  (case rn_osm_mle=0) 
     1253   rn_osm_mle_lat      = 20.       ! reference latitude (degrees) of MLE coef.                        (case rn_mle=1) 
     1254   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
     1255   rn_osm_mle_thresh  = 0.0005     ! delta b criterion used for FK MLE criterion 
     1256   rn_osm_mle_tau     = 172800.    ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
     1257   ln_osm_hmle_limit   = .false.   ! limit hmle to rn_osm_hmle_limit*hbl 
     1258   rn_osm_hmle_limit   = 1.2 
     1259   / 
    12371260!----------------------------------------------------------------------- 
    12381261&namzdf_mfc     !   Mass Flux Convection 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/ICE/icerst.F90

    r14018 r14051  
    9999               CALL iom_swap( cxios_context ) 
    100100#else 
    101                clinfo = 'Can not use XIOS in rst_opn' 
    102                CALL ctl_stop(TRIM(clinfo)) 
     101               CALL ctl_stop( 'Can not use XIOS in rst_opn' ) 
    103102#endif 
    104103            ENDIF 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA/diawri.F90

    r13998 r14051  
    4747   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
    4848   USE zdfmxl         ! mixed layer 
     49   USE zdfosm         ! mixed layer 
    4950   ! 
    5051   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    … …  
    11621163         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
    11631164      ENDIF 
     1165      IF( ln_zdfosm ) THEN 
     1166         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)  )      ! now boundary-layer depth 
     1167         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)  )      ! now mixed-layer depth 
     1168         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask     )      ! w-level diffusion 
     1169         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask     )      ! now w-level viscosity 
     1170         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask     )      ! non-local t forcing 
     1171         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask     )      ! non-local s forcing 
     1172         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask     )      ! non-local u forcing 
     1173         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask     )      ! non-local v forcing 
     1174         IF( ln_osm_mle ) THEN 
     1175            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)  ) ! now transition-layer depth 
     1176         END IF 
     1177      ENDIF 
    11641178      ! 
    11651179      CALL iom_close( inum ) 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/iom.F90

    r14023 r14051  
    2929   USE in_out_manager  ! I/O manager 
    3030   USE lib_mpp           ! MPP library 
    31 #if defined key_iomput 
    3231   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 
    3332   USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes 
    … …  
    3736   USE phycst          ! physical constants 
    3837   USE dianam          ! build name of file 
     38#if defined key_iomput 
    3939   USE xios 
    4040# endif 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/TRA/traatf.F90

    r14023 r14051  
    117117      IF( l_trdtra )   THEN                     
    118118         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    119          ztrdt(:,:,jpk) = 0._wp 
    120          ztrds(:,:,jpk) = 0._wp 
     119         ztrdt(:,:,:) = 0._wp 
     120         ztrds(:,:,:) = 0._wp 
    121121         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
    122122            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/TRA/tramle.F90

    r14023 r14051  
    2020   USE lib_mpp        ! MPP library 
    2121   USE lbclnk         ! lateral boundary condition / mpp link 
     22 
     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 
    2226 
    2327   IMPLICIT NONE 
    … …  
    99103      !!---------------------------------------------------------------------- 
    100104      ! 
    101       !                                      !==  MLD used for MLE  ==! 
    102       !                                                ! compute from the 10m density to deal with the diurnal cycle 
    103       DO_2D( 1, 1, 1, 1 ) 
    104          inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1) 
    105       END_2D 
    106       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
    107          DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m) 
    108             IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
     105      ! 
     106      IF(ln_osm_mle.and.ln_zdfosm) THEN 
     107         ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     108         ! 
     109         ! 
     110         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     111         CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     112            DO_2D( 1, 0, 1, 0 ) 
     113               zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 
     114               zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 
     115            END_2D 
     116         CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     117            DO_2D( 1, 0, 1, 0 ) 
     118               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
     119               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
     120            END_2D 
     121         CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     122            DO_2D( 1, 0, 1, 0 ) 
     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_2D 
     126         END SELECT 
     127         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     128            DO_2D( 1, 0, 1, 0 ) 
     129               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            & 
     130                    &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     131                    &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     132               ! 
     133               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            & 
     134                    &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     135                    &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     136            END_2D 
     137            ! 
     138         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     139            DO_2D( 1, 0, 1, 0 ) 
     140               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               & 
     141                    &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     142               ! 
     143               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               & 
     144                    &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     145            END_2D 
     146         ENDIF 
     147 
     148      ELSE !do not use osn_mle 
     149         !                                      !==  MLD used for MLE  ==! 
     150         !                                                ! compute from the 10m density to deal with the diurnal cycle 
     151         DO_2D( 1, 1, 1, 1 ) 
     152            inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1) 
     153         END_2D 
     154         IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
     155           DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m) 
     156              IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
     157           END_3D 
     158         ENDIF 
     159         ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
     160         ! 
     161         ! 
     162         zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
     163         zbm (:,:) = 0._wp 
     164         zn2 (:,:) = 0._wp 
     165         DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer 
     166            zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     167            zmld(ji,jj) = zmld(ji,jj) + zc 
     168            zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 
     169            zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
    109170         END_3D 
    110       ENDIF 
    111       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
    112       ! 
    113       ! 
    114       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
    115       zbm (:,:) = 0._wp 
    116       zn2 (:,:) = 0._wp 
    117       DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer 
    118          zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
    119          zmld(ji,jj) = zmld(ji,jj) + zc 
    120          zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 
    121          zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
    122       END_3D 
    123  
    124       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    125       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
    126          DO_2D( 1, 0, 1, 0 ) 
    127             zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
    128             zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
     171    
     172         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     173         CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     174            DO_2D( 1, 0, 1, 0 ) 
     175               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
     176               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
     177            END_2D 
     178         CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     179            DO_2D( 1, 0, 1, 0 ) 
     180               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
     181               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
     182            END_2D 
     183         CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     184            DO_2D( 1, 0, 1, 0 ) 
     185               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
     186               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
     187            END_2D 
     188         END SELECT 
     189         !                                                ! convert density into buoyancy 
     190         DO_2D( 1, 1, 1, 1 ) 
     191            zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) ) 
    129192         END_2D 
    130       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
    131          DO_2D( 1, 0, 1, 0 ) 
    132             zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
    133             zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
    134          END_2D 
    135       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
    136          DO_2D( 1, 0, 1, 0 ) 
    137             zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
    138             zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
    139          END_2D 
    140       END SELECT 
    141       !                                                ! convert density into buoyancy 
    142       DO_2D( 1, 1, 1, 1 ) 
    143          zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) ) 
    144       END_2D 
    145       ! 
    146       ! 
    147       !                                      !==  Magnitude of the MLE stream function  ==! 
    148       ! 
    149       !                 di[bm]  Ds 
    150       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
    151       !                  e1u   Lf fu                      and the e2u for the "transport" 
    152       !                                                      (not *e3u as divided by e3u at the end) 
    153       ! 
    154       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
    155          DO_2D( 1, 0, 1, 0 ) 
    156             zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
    157                &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
    158                &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     193         ! 
     194         ! 
     195         !                                      !==  Magnitude of the MLE stream function  ==! 
     196         ! 
     197         !                 di[bm]  Ds 
     198         ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
     199         !                  e1u   Lf fu                      and the e2u for the "transport" 
     200         !                                                      (not *e3u as divided by e3u at the end) 
     201         ! 
     202         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     203            DO_2D( 1, 0, 1, 0 ) 
     204               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
     205                    &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     206                    &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
    159207               ! 
    160             zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
    161                &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
    162                &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
    163          END_2D 
    164          ! 
    165       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
    166          DO_2D( 1, 0, 1, 0 ) 
    167             zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
    168                &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     208               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
     209                    &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     210                    &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     211            END_2D 
     212            ! 
     213         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     214            DO_2D( 1, 0, 1, 0 ) 
     215               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
     216                    &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    169217               ! 
    170             zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
    171                &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
    172          END_2D 
    173       ENDIF 
    174       ! 
    175       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
    176          DO_2D( 1, 0, 1, 0 ) 
    177             IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
    178             IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
    179          END_2D 
    180       ENDIF 
    181       ! 
    182       !                                      !==  structure function value at uw- and vw-points  ==! 
    183       DO_2D( 1, 0, 1, 0 ) 
    184          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu 
    185          zhv(ji,jj) = 1._wp / zhv(ji,jj) 
    186       END_2D 
    187       ! 
    188       zpsi_uw(:,:,:) = 0._wp 
    189       zpsi_vw(:,:,:) = 0._wp 
    190       ! 
     218               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
     219                    &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     220            END_2D 
     221         ENDIF 
     222         ! 
     223         IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
     224            DO_2D( 1, 0, 1, 0 ) 
     225               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
     226               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
     227            END_2D 
     228         ENDIF 
     229         ! 
     230      ENDIF  ! end of ln_osm_mle conditional 
     231    !                                      !==  structure function value at uw- and vw-points  ==! 
     232    DO_2D( 1, 0, 1, 0 ) 
     233       zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu 
     234       zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall)  
     235    END_2D 
     236    ! 
     237    zpsi_uw(:,:,:) = 0._wp 
     238    zpsi_vw(:,:,:) = 0._wp 
     239    ! 
    191240      DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0 
    192241         zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) 
    … …  
    220269         ENDIF 
    221270         ! 
    222          DO_2D( 0, 0, 0, 0 ) 
    223             zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f 
    224          END_2D 
     271         IF (ln_osm_mle.and.ln_zdfosm) THEN 
     272            DO_2D( 0, 0, 0, 0 ) 
     273               zLf_NH(ji,jj) = SQRT( rb_c * hmle(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f 
     274            END_2D 
     275         ELSE 
     276            DO_2D( 0, 0, 0, 0 ) 
     277               zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f 
     278            END_2D 
     279         ENDIF 
    225280         ! 
    226281         ! divide by cross distance to give streamfunction with dimensions m^2/s 
    … …  
    239294      ! 
    240295   END SUBROUTINE tra_mle_trp 
    241  
    242296 
    243297   SUBROUTINE tra_mle_init 
    … …  
    301355            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' ) 
    302356            z1_t2 = 1._wp / ( rn_time * rn_time ) 
    303             DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )                      ! "coriolis+ time^-1" at u- & v-points 
     357            DO_2D( 0, 1, 0, 1 )                      ! "coriolis+ time^-1" at u- & v-points 
    304358               zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 
    305359               zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp 
    … …  
    307361               rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 ) 
    308362            END_2D 
    309             IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp ) 
     363            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp ) 
    310364            ! 
    311365         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/TRD/trd_oce.F90

    r10068 r14051  
    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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfosm.F90

    r14023 r14051  
    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 
     105   REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by 
     106   REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl 
     107   LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering 
    86108   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs 
    87109   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt 
    88110   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave 
     111   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value 
    89112   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la 
    90113 
    … …  
    96119   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s) 
    97120 
     121! OSMOSIS mixed layer eddy parametrization constants 
     122   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt 
     123   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient 
     124   !                                        ! parameters used in nn_osm_mle = 0 case 
     125   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front 
     126   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer 
     127   !                                        ! parameters used in nn_osm_mle = 1 case 
     128   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
     129   LOGICAL  ::   ln_osm_hmle_limit           ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
     130   REAL(wp) ::   rn_osm_hmle_limit           ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
     131   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
     132   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
     133   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 
     134   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
     135   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
     136   REAL(wp) ::   rn_osm_bl_thresh          ! Threshold buoyancy for deepening of OSBL base. 
     137   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
     138 
     139 
    98140   !                                    !!! ** General constants  ** 
    99    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number 
     141   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
     142   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    100143   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    101144   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
    … …  
    118161      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    119162      !!---------------------------------------------------------------------- 
    120      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 
    121           &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     & 
    122           &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc ) 
     163     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
     164          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     165          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     166 
     167     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 
     168          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 
     169 
     170     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
    123171     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
    124      CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     172 
    125173   END FUNCTION zdf_osm_alloc 
    126174 
    … …  
    166214      !! 
    167215      INTEGER ::   ji, jj, jk                   ! dummy loop indices 
     216 
     217      INTEGER ::   jl                   ! dummy loop indices 
     218 
    168219      INTEGER ::   ikbot, jkmax, jkm1, jkp2     ! 
    169220 
    … …  
    196247      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    197248      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     249      REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 
     250 
     251 
     252      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL 
     253      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux 
     254      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 
     255      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK 
     256 
    198257      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift 
    199258      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number 
    … …  
    201260      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 
    202261      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    203       LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 
     262      LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl 
     263      LOGICAL, DIMENSION(jpi,jpj)  :: lshear    ! Shear layers 
     264      LOGICAL, DIMENSION(jpi,jpj)  :: lpyc      ! OSBL pycnocline present 
     265      LOGICAL, DIMENSION(jpi,jpj)  :: lflux     ! surface flux extends below OSBL into MLE layer. 
     266      LOGICAL, DIMENSION(jpi,jpj)  :: lmle      ! MLE layer increases in hickness. 
    204267 
    205268      ! mixed-layer variables 
    … …  
    207270      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 
    208271      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 
     272      INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level 
     273      INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 
    209274 
    210275      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients 
    … …  
    213278      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    214279      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid 
     280 
     281      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 
     282      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid 
     283 
    215284      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    216285      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
    217       REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    218       REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
    219       REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 
    220       REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 
    221       REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
    222       REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
     286      REAL(wp), DIMENSION(jpi,jpj) :: zddhdt                                    ! correction to dhdt due to internal structure. 
     287      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients 
     288      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext              ! external temperature/salinity and buoyancy gradients 
     289      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
     290 
     291      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl  ! averages over the depth of the blayer 
     292      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml  ! averages over the depth of the mixed layer 
     293      REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle  ! averages over the depth of the MLE layer 
     294      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 
     295      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 
     296      REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer 
     297!      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
     298      REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
     299      REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
    223300      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline 
    224301      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline 
    … …  
    226303      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline 
    227304      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline 
    228  
     305      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
    229306      ! Flux-gradient relationship variables 
     307      REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 
    230308 
    231309      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale. 
    232310 
    233       REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity 
     311      REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.   
    234312      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. 
     313      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ 
    235314      REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. 
    236315      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep 
    … …  
    243322      ! Temporary variables 
    244323      INTEGER :: inhml 
    245       INTEGER :: i_lconv_alloc 
    246324      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 
    247325      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables 
    248326      REAL(wp) :: zthick, zz0, zz1 ! temporary variables 
    249327      REAL(wp) :: zvel_max, zhbl_s ! temporary variables 
    250       REAL(wp) :: zfac             ! temporary variable 
     328      REAL(wp) :: zfac, ztmp       ! temporary variable 
    251329      REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift 
    252330      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity 
    253331      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
     332      REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 
     333      REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. 
     334      REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc 
     335      REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max 
     336      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
     337      REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 
     338      REAL(wp) :: zzeta_s = 0._wp 
     339      REAL(wp) :: zzeta_v = 0.46 
     340      REAL(wp) :: zabsstke 
     341      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 
     342      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 
    254343 
    255344      ! For debugging 
    … …  
    257346      !!-------------------------------------------------------------------- 
    258347      ! 
    259       ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc ) 
    260       IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv') 
    261  
    262348      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    263349      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
    … …  
    267353      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 
    268354      zhol(:,:)   = 0._wp 
    269       lconv(:,:)  = .FALSE. 
     355      lconv(:,:)  = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE. 
    270356      ! mixed layer 
    271357      ! no initialization of zhbl or zhml (or zdh?) 
    272358      zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp 
    273       zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    274       zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
    275       zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    276       zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     359      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp 
     360      zv_bl(:,:)   = 0._wp ; zb_bl(:,:)  = 0._wp 
     361      zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     362      zt_mle(:,:)   = 0._wp ; zs_mle(:,:)    = 0._wp ; zu_mle(:,:)   = 0._wp 
     363      zb_mle(:,:) = 0._wp 
     364      zv_ml(:,:)   = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
     365      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
    277366      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp 
    278       zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp 
    279       zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp 
     367      zdb_ml(:,:)  = 0._wp 
     368      zdt_mle(:,:)  = 0._wp ; zds_mle(:,:)  = 0._wp ; zdu_mle(:,:)   = 0._wp 
     369      zdv_mle(:,:)  = 0._wp ; zdb_mle(:,:)  = 0._wp 
     370      zwth_ent = 0._wp ; zws_ent = 0._wp 
    280371      ! 
    281372      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 
    282373      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    283374      ! 
     375      zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 
     376 
     377      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
     378         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 
     379         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 
     380         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 
     381         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp 
     382      ENDIF 
     383      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
     384 
    284385      ! Flux-Gradient arrays. 
    285       zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
    286       zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp 
    287386      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp 
    288387      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp 
    … …  
    292391      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    293392 
     393      zddhdt(:,:) = 0._wp 
    294394      ! hbl = MAX(hbl,epsln) 
    295395      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    326426        zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    327427        ! Surface upward velocity fluxes 
    328         zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
    329         zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
     428        zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
     429        zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
    330430        ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    331431        zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
    … …  
    340440           zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
    341441           zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
     442           ! Linearly 
    342443           zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 
    343            ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 
     444           dstokes(ji,jj) = rn_osm_dstokes 
    344445        END_2D 
    345446     ! Assume Pierson-Moskovitz wind-wave spectrum 
    … …  
    347448        DO_2D( 0, 0, 0, 0 ) 
    348449           ! Use wind speed wndm included in sbc_oce module 
    349            zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    350            dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
     450           zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     451           dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    351452        END_2D 
    352453     ! Use ECMWF wave fields as output from SBCWAVE 
    353454     CASE(2) 
    354455        zfac =  2.0_wp * rpi / 16.0_wp 
     456 
    355457        DO_2D( 0, 0, 0, 0 ) 
    356            ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 
    357            !    The coefficient 0.8 gives La=0.3  in this situation. 
    358            ! It could represent the effects of the spread of wave directions 
    359            ! around the mean wind. The effect of this adjustment needs to be tested. 
    360            zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
    361                 &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
    362            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 ! 
     458           IF (hsw(ji,jj) > 1.e-4) THEN 
     459              ! Use  wave fields 
     460              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     461              zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     462              dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) 
     463           ELSE 
     464              ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 
     465              ! .. so default to Pierson-Moskowitz 
     466              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     467              dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
     468           END IF 
     469        END_2D 
     470     END SELECT 
     471 
     472     IF (ln_zdfosm_ice_shelter) THEN 
     473        ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 
     474        DO_2D( 0, 0, 0, 0 ) 
     475           zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     476           dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     477        END_2D 
     478     END IF 
     479 
     480     SELECT CASE (nn_osm_SD_reduce) 
     481     ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van  Roekel (2012) or Grant (2020). 
     482     CASE(0) 
     483        ! The Langmur number from the ECMWF model (or from PM)  appears to give La<0.3 for wind-driven seas. 
     484        !    The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3  in this situation. 
     485        ! It could represent the effects of the spread of wave directions 
     486        ! around the mean wind. The effect of this adjustment needs to be tested. 
     487        IF(nn_osm_wave > 0) THEN 
     488           zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) 
     489        END IF 
     490     CASE(1) 
     491        ! van  Roekel (2012): consider average SD over top 10% of boundary layer 
     492        ! assumes approximate depth profile of SD from Breivik (2016) 
     493        zsqrtpi = SQRT(rpi) 
     494        z_two_thirds = 2.0_wp / 3.0_wp 
     495 
     496        DO_2D( 0, 0, 0, 0 ) 
     497           zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     498           z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     499           zsqrt_depth = SQRT(z2k_times_thickness) 
     500           zexp_depth  = EXP(-z2k_times_thickness) 
     501           zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth  & 
     502                &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & 
     503                &              + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness 
     504 
     505        END_2D 
     506     CASE(2) 
     507        ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer 
     508        ! assumes approximate depth profile of SD from Breivik (2016) 
     509        zsqrtpi = SQRT(rpi) 
     510 
     511        DO_2D( 0, 0, 0, 0 ) 
     512           zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     513           z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     514 
     515           IF(z2k_times_thickness < 50._wp) THEN 
     516              zsqrt_depth = SQRT(z2k_times_thickness) 
     517              zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) 
     518           ELSE 
     519              ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness 
     520              ! See Abramowitz and Stegun, Eq. 7.1.23 
     521              ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness)  + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 
     522              zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp 
     523           END IF 
     524           zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 
     525           dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 
     526           zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 
    363527        END_2D 
    364528     END SELECT 
    … …  
    369533        ! Langmuir velocity scale (zwstrl), at T-point 
    370534        zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
    371         ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
    372         ! Intended as a possible test. Doesn't affect LES results for entrainment, 
    373         !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
    374         zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
    375              & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
    376              & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
    377         ! define La this way so effects of Stokes penetration depth on velocity scale are included 
    378         zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
     535        zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 
     536        IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 
    379537        ! Velocity scale that tends to zustar for large Langmuir numbers 
    380538        zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
    … …  
    383541        ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    384542        ! Note zustke and zwstrl are not amended. 
    385         IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    386543        ! 
    387544        ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
    … …  
    389546           zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 
    390547           zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 
    391            lconv(ji,jj) = .TRUE. 
    392         ELSE 
     548         ELSE 
    393549           zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln ) 
    394            lconv(ji,jj) = .FALSE. 
    395550        ENDIF 
    396551     END_2D 
    … …  
    399554     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 
    400555     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    401      ! BL must be always 2 levels deep. 
    402       hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 
    403       ibld(:,:) = 3 
    404       DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 
     556     ! BL must be always 4 levels deep. 
     557     ! For calculation of lateral buoyancy gradients for FK in 
     558     ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 
     559     ! previously exist for hbl also. 
     560 
     561     ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 
     562     ! ########################################################################## 
     563      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 
     564      ibld(:,:) = 4 
     565      DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 
    405566         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    406567            ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
    407568         ENDIF 
    408569      END_3D 
     570     ! ########################################################################## 
    409571 
    410572      DO_2D( 0, 0, 0, 0 ) 
    411             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    412             zbeta    = rab_n(ji,jj,1,jp_sal) 
    413             zt   = 0._wp 
    414             zs   = 0._wp 
    415             zu   = 0._wp 
    416             zv   = 0._wp 
    417             ! average over depth of boundary layer 
    418             zthick=0._wp 
    419             DO jm = 2, ibld(ji,jj) 
    420                zthick=zthick+e3t(ji,jj,jm,Kmm) 
    421                zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
    422                zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
    423                zu   = zu  + e3t(ji,jj,jm,Kmm) & 
    424                   &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    425                   &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    426                zv   = zv  + e3t(ji,jj,jm,Kmm) & 
    427                   &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    428                   &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    429             END DO 
    430             zt_bl(ji,jj) = zt / zthick 
    431             zs_bl(ji,jj) = zs / zthick 
    432             zu_bl(ji,jj) = zu / zthick 
    433             zv_bl(ji,jj) = zv / zthick 
    434             zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
    435             zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
    436             zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    437                   &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    438             zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    439                   &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    440             zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    441             IF ( lconv(ji,jj) ) THEN    ! Convective 
    442                    zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    443                         &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    444  
    445                    zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    446                         &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    447 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    448 !                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    449 !                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    450  
    451 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    452 !                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    453                    zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    454             ELSE                        ! Stable 
    455                    zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
    456                         &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
    457                         & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
    458                         &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
    459                    zzdhdt = zzdhdt + zwbav(ji,jj) 
    460                    IF ( zzdhdt < 0._wp ) THEN 
    461                    ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    462                       zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    463                    ELSE 
    464                       zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    465                            &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    466                    ENDIF 
    467                    zzdhdt = 2.0 * zzdhdt / zpert 
    468             ENDIF 
    469             zdhdt(ji,jj) = zzdhdt 
     573         zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
     574         imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 )) 
     575         zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
     576         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    470577      END_2D 
    471  
    472       ! Calculate averages over depth of boundary layer 
    473       imld = ibld           ! use imld to hold previous blayer index 
    474       ibld(:,:) = 3 
    475  
    476       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 
    477       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 
    478       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     578      ! Averages over well-mixed and boundary layer 
     579      jp_ext(:,:) = 2 
     580      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     581!      jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 
     582      CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     583! Velocity components in frame aligned with surface stress. 
     584      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     585      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
     586! Determine the state of the OSBL, stable/unstable, shear/no shear 
     587      CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
     588 
     589      IF ( ln_osm_mle ) THEN 
     590! Fox-Kemper Scheme 
     591         mld_prof = 4 
     592         DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 
     593         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     594         END_3D 
     595         jp_ext_mle(:,:) = 2 
     596        CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) 
     597 
     598         DO_2D( 0, 0, 0, 0 ) 
     599           zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
     600         END_2D 
     601 
     602!! External gradient 
     603         CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     604         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
     605         CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 
     606         CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
     607         CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
     608      ELSE    ! ln_osm_mle 
     609! FK not selected, Boundary Layer only. 
     610         lpyc(:,:) = .TRUE. 
     611         lflux(:,:) = .FALSE. 
     612         lmle(:,:) = .FALSE. 
     613         DO_2D( 0, 0, 0, 0 ) 
     614          IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 
     615         END_2D 
     616      ENDIF   ! ln_osm_mle 
     617 
     618! Test if pycnocline well resolved 
     619      DO_2D( 0, 0, 0, 0 ) 
     620       IF (lconv(ji,jj) ) THEN 
     621          ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) 
     622          IF ( ztmp > 6 ) THEN 
     623   ! pycnocline well resolved 
     624            jp_ext(ji,jj) = 1 
     625          ELSE 
     626   ! pycnocline poorly resolved 
     627            jp_ext(ji,jj) = 0 
     628          ENDIF 
     629       ELSE 
     630   ! Stable conditions 
     631         jp_ext(ji,jj) = 0 
     632       ENDIF 
     633      END_2D 
     634 
     635      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     636!      jp_ext = ibld-imld+1 
     637      CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     638! Rate of change of hbl 
     639      CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     640      DO_2D( 0, 0, 0, 0 ) 
     641       zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it 
     642            ! adjustment to represent limiting by ocean bottom 
     643       IF ( zhbl_t(ji,jj) >= gdepw(ji, jj, mbkt(ji,jj) + 1, Kmm ) ) THEN 
     644          zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm) - depth_tol)! ht(:,:)) 
     645          lpyc(ji,jj) = .FALSE. 
     646       ENDIF 
     647      END_2D 
     648 
     649      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index 
     650      ibld(:,:) = 4 
    479651 
    480652      DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 
    481653         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    482             ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     654            ibld(ji,jj) = jk 
    483655         ENDIF 
    484656      END_3D 
    … …  
    487659! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    488660! 
     661      CALL zdf_osm_timestep_hbl( zdhdt ) 
     662! is external level in bounds? 
     663 
     664      CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     665! 
     666! 
     667! Check to see if lpyc needs to be changed  
     668 
     669      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
     670 
    489671      DO_2D( 0, 0, 0, 0 ) 
    490          IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     672       IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE.  
     673      END_2D 
     674 
     675      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    491676! 
    492 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    493 ! 
    494             zhbl_s = hbl(ji,jj) 
    495             jm = imld(ji,jj) 
    496             zthermal = rab_n(ji,jj,1,jp_tem) 
    497             zbeta = rab_n(ji,jj,1,jp_sal) 
    498             IF ( lconv(ji,jj) ) THEN 
    499 !unstable 
    500                zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    501                     &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    502  
    503                DO jk = imld(ji,jj), ibld(ji,jj) 
    504                   zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 
    505                        & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 
    506  
    507                   zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ),   & 
    508                      &                     e3w(ji,jj,jk,Kmm) ) 
    509                   zhbl_s = MIN(zhbl_s, ht(ji,jj)) 
    510  
    511                   IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
    512                END DO 
    513                hbl(ji,jj) = zhbl_s 
    514                ibld(ji,jj) = jm 
    515                hbli(ji,jj) = hbl(ji,jj) 
    516             ELSE 
    517 ! stable 
    518                DO jk = imld(ji,jj), ibld(ji,jj) 
    519                   zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          & 
    520                        &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 
    521                        & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    522  
    523                   zhbl_s = zhbl_s +  (                                                                                & 
    524                        &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
    525                        &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
    526                        &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    527                        &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    528                        &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here 
    529  
    530                   zhbl_s = MIN(zhbl_s, ht(ji,jj)) 
    531                   IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
    532                END DO 
    533                hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 
    534                ibld(ji,jj) = MAX(jm, 3 ) 
    535                IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    536             ENDIF   ! IF ( lconv ) 
    537          ELSE 
    538 ! change zero or one model level. 
    539             hbl(ji,jj) = zhbl_t(ji,jj) 
    540             IF ( lconv(ji,jj) ) THEN 
    541                hbli(ji,jj) = hbl(ji,jj) 
    542             ELSE 
    543                hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 
    544                IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    545             ENDIF 
    546          ENDIF 
    547          zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    548       END_2D 
    549       dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    550  
    551 ! Recalculate averages over boundary layer after depth updated 
    552      ! Consider later  combining this into the loop above and looking for columns 
    553      ! where the index for base of the boundary layer have changed 
    554       DO_2D( 0, 0, 0, 0 ) 
    555             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    556             zbeta    = rab_n(ji,jj,1,jp_sal) 
    557             zt   = 0._wp 
    558             zs   = 0._wp 
    559             zu   = 0._wp 
    560             zv   = 0._wp 
    561             ! average over depth of boundary layer 
    562             zthick=0._wp 
    563             DO jm = 2, ibld(ji,jj) 
    564                zthick=zthick+e3t(ji,jj,jm,Kmm) 
    565                zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
    566                zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
    567                zu   = zu  + e3t(ji,jj,jm,Kmm) & 
    568                   &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    569                   &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    570                zv   = zv  + e3t(ji,jj,jm,Kmm) & 
    571                   &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    572                   &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    573             END DO 
    574             zt_bl(ji,jj) = zt / zthick 
    575             zs_bl(ji,jj) = zs / zthick 
    576             zu_bl(ji,jj) = zu / zthick 
    577             zv_bl(ji,jj) = zv / zthick 
    578             zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
    579             zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
    580             zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    581                    &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    582             zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    583                    &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    584             zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    585             zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    586             IF ( lconv(ji,jj) ) THEN 
    587                IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    588                   IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    589                         zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
    590                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    591                   ELSE                                                     ! unstable 
    592                         zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
    593                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    594                   ENDIF 
    595                   IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
    596                      zari = 0.2 
    597                      zwb_ent(ji,jj) = 0._wp 
    598                   ENDIF 
    599                   inhml = MAX( INT( zari * zhbl(ji,jj)   & 
    600                      &              / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 
    601                   imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    602                   zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    603                   zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    604                ELSE  ! IF (zdb_bl) 
    605                   imld(ji,jj) = ibld(ji,jj) - 1 
    606                   zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    607                   zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    608                ENDIF 
    609             ELSE   ! IF (lconv) 
    610                IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    611                ! boundary layer deepening 
    612                   IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    613                ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    614                      zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    615                        & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    616                      inhml = MAX( INT( zari * zhbl(ji,jj)   & 
    617                         &             / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 
    618                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    619                      zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    620                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    621                   ELSE 
    622                      imld(ji,jj) = ibld(ji,jj) - 1 
    623                      zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    624                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    625                   ENDIF ! IF (zdb_bl > 0.0) 
    626                ELSE     ! IF(dhdt >= 0) 
    627                ! boundary layer collapsing. 
    628                   imld(ji,jj) = ibld(ji,jj) 
    629                   zhml(ji,jj) = zhbl(ji,jj) 
    630                   zdh(ji,jj) = 0._wp 
    631                ENDIF    ! IF (dhdt >= 0) 
    632             ENDIF       ! IF (lconv) 
    633       END_2D 
    634  
    635       ! Average over the depth of the mixed layer in the convective boundary layer 
    636       ! Also calculate entrainment fluxes for temperature and salinity 
    637       DO_2D( 0, 0, 0, 0 ) 
    638          zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    639          zbeta    = rab_n(ji,jj,1,jp_sal) 
    640          IF ( lconv(ji,jj) ) THEN 
    641             zt   = 0._wp 
    642             zs   = 0._wp 
    643             zu   = 0._wp 
    644             zv   = 0._wp 
    645             ! average over depth of boundary layer 
    646             zthick=0._wp 
    647             DO jm = 2, imld(ji,jj) 
    648                zthick=zthick+e3t(ji,jj,jm,Kmm) 
    649                zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
    650                zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
    651                zu   = zu  + e3t(ji,jj,jm,Kmm) & 
    652                   &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    653                   &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    654                zv   = zv  + e3t(ji,jj,jm,Kmm) & 
    655                   &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    656                   &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    657             END DO 
    658             zt_ml(ji,jj) = zt / zthick 
    659             zs_ml(ji,jj) = zs / zthick 
    660             zu_ml(ji,jj) = zu / zthick 
    661             zv_ml(ji,jj) = zv / zthick 
    662             zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
    663             zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
    664             zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    665                   &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    666             zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    667                   &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    668             zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    669          ELSE 
    670          ! stable, if entraining calulate average below interface layer. 
    671             IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    672                zt   = 0._wp 
    673                zs   = 0._wp 
    674                zu   = 0._wp 
    675                zv   = 0._wp 
    676                ! average over depth of boundary layer 
    677                zthick=0._wp 
    678                DO jm = 2, imld(ji,jj) 
    679                   zthick=zthick+e3t(ji,jj,jm,Kmm) 
    680                   zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
    681                   zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
    682                   zu   = zu  + e3t(ji,jj,jm,Kmm) & 
    683                      &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    684                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    685                   zv   = zv  + e3t(ji,jj,jm,Kmm) & 
    686                      &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    687                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    688                END DO 
    689                zt_ml(ji,jj) = zt / zthick 
    690                zs_ml(ji,jj) = zs / zthick 
    691                zu_ml(ji,jj) = zu / zthick 
    692                zv_ml(ji,jj) = zv / zthick 
    693                zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
    694                zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
    695                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    696                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    697                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    698                      &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    699                zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    700             ENDIF 
    701          ENDIF 
    702       END_2D 
    703     ! 
     677    ! Average over the depth of the mixed layer in the convective boundary layer 
     678!      jp_ext = ibld - imld +1 
     679      CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    704680    ! rotate mean currents and changes onto wind align co-ordinates 
    705681    ! 
    706  
    707       DO_2D( 0, 0, 0, 0 ) 
    708          ztemp = zu_ml(ji,jj) 
    709          zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
    710          zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    711          ztemp = zdu_ml(ji,jj) 
    712          zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
    713          zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    714  ! 
    715          ztemp = zu_bl(ji,jj) 
    716          zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
    717          zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    718          ztemp = zdu_bl(ji,jj) 
    719          zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
    720          zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    721       END_2D 
    722  
    723      zuw_bse = 0._wp 
    724      zvw_bse = 0._wp 
    725      DO_2D( 0, 0, 0, 0 ) 
    726  
    727         IF ( lconv(ji,jj) ) THEN 
    728            IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    729               zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    730               zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    731            ENDIF 
    732         ELSE 
    733            zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    734            zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    735         ENDIF 
    736      END_2D 
    737  
     682     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     683     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    738684      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    739685      !  Pycnocline gradients for scalars and velocity 
    740686      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    741687 
    742        DO_2D( 0, 0, 0, 0 ) 
    743        ! 
    744           IF ( lconv (ji,jj) ) THEN 
    745           ! Unstable conditions 
    746              IF( zdb_bl(ji,jj) > 0._wp ) THEN 
    747              ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
    748                 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
    749                 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
    750                 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    751                 DO jk = 2 , ibld(ji,jj) 
    752                    znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    753                    zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    754                    zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    755                    zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    756                 END DO 
    757              ENDIF 
    758           ELSE 
    759           ! stable conditions 
    760           ! if pycnocline profile only defined when depth steady of increasing. 
    761              IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
    762                 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    763                   IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    764                       ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
    765                       zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
    766                       zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    767                       DO jk = 2, ibld(ji,jj) 
    768                          znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    769                          zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    770                          zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    771                          zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    772                       END DO 
    773                   ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    774                       ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
    775                       zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
    776                       zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    777                       DO jk = 2, ibld(ji,jj) 
    778                          znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    779                          zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    780                          zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    781                          zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    782                       END DO 
    783                    ENDIF ! IF (zhol >=0.5) 
    784                 ENDIF    ! IF (zdb_bl> 0.) 
    785              ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
    786           ENDIF          ! IF (lconv) 
    787          ! 
    788        END_2D 
    789 ! 
    790        DO_2D( 0, 0, 0, 0 ) 
    791        ! 
    792           IF ( lconv (ji,jj) ) THEN 
    793           ! Unstable conditions 
    794               zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
    795             & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
    796              zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
    797            & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    798              DO jk = 2 , ibld(ji,jj)-1 
    799                 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    800                 zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    801                 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    802              END DO 
    803           ELSE 
    804           ! stable conditions 
    805              zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    806              zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    807              DO jk = 2, ibld(ji,jj) 
    808                 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    809                 IF ( znd < 1.0 ) THEN 
    810                    zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    811                 ELSE 
    812                    zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    813                 ENDIF 
    814                 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
    815              END DO 
    816           ENDIF 
    817          ! 
    818        END_2D 
     688      CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     689      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 
     690      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    819691       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    820692       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    821693       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    822  
    823       ! WHERE ( lconv ) 
    824       !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 
    825       !     zvisml_sc = zdifml_sc 
    826       !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 
    827       !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 
    828       !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 
    829       !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 
    830       !  ELSEWHERE 
    831       !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    832       !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    833       !  ENDWHERE 
    834        DO_2D( 0, 0, 0, 0 ) 
    835           IF ( lconv(ji,jj) ) THEN 
    836             zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    837             zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 
    838             zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    839             zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    840             zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    841             zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
    842           ELSE 
    843             zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    844             zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    845          END IF 
    846        END_2D 
    847 ! 
    848        DO_2D( 0, 0, 0, 0 ) 
    849           IF ( lconv(ji,jj) ) THEN 
    850              DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
    851                  zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    852                  ! 
    853                  zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5 
    854                  ! 
    855                  zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) & 
    856                       &            *                                      ( 1.0 -               0.5 * zznd_ml**2 ) 
    857              END DO 
    858              ! pycnocline - if present linear profile 
    859              IF ( zdh(ji,jj) > 0._wp ) THEN 
    860                 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    861                     zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    862                     ! 
    863                     zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    864                     ! 
    865                     zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    866                 END DO 
    867              ENDIF 
    868              ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 
    869              zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 
    870              ! could be taken out, take account of entrainment represents as a diffusivity 
    871              ! should remove w from here, represents entrainment 
    872           ELSE 
    873           ! stable conditions 
    874              DO jk = 2, ibld(ji,jj) 
    875                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    876                 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
    877                 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    878              END DO 
    879           ENDIF   ! end if ( lconv ) 
    880 ! 
    881        END_2D 
     694       CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 
    882695 
    883696       ! 
    … …  
    918731       END_2D 
    919732 
    920  
    921733! 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) 
    922734       WHERE ( lconv ) 
    923           zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ) 
    924           zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln ) 
    925           zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
     735          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 ) 
     736          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
     737          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 ) 
    926738       ELSEWHERE 
    927739          zsc_uw_1 = zustar**2 
    928           zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 
     740          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    929741       ENDWHERE 
    930  
     742       IF(ln_dia_osm) THEN 
     743          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     744          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     745       END IF 
    931746       DO_2D( 0, 0, 0, 0 ) 
    932747          IF ( lconv(ji,jj) ) THEN 
    … …  
    970785                zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    971786                     &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    972                 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     787                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 
    973788                ! non-gradient buoyancy terms 
    974789                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 ) 
    975790                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    976791             END DO 
    977           ELSE 
     792              
     793             IF ( lpyc(ji,jj) ) THEN 
     794               ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     795               ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) 
     796               zwth_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)                   
     797               zws_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 
     798! Cubic profile used for buoyancy term 
     799               za_cubic = 0.755 * ztau_sc_u(ji,jj) 
     800               zb_cubic = 0.25 * ztau_sc_u(ji,jj) 
     801               DO jk = 2, ibld(ji,jj) 
     802                 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
     803                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
     804 
     805                 ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
     806               END DO 
     807! 
     808               zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 
     809               zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) ) 
     810! 
     811               zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
     812! 
     813               zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
     814! 
     815               zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )  
     816               DO jk = 2, ibld(ji,jj) 
     817                 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
     818                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     819! 
     820                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     821               END DO 
     822            ENDIF ! End of pycnocline                   
     823          ELSE ! lconv test - stable conditions 
    978824             DO jk = 2, ibld(ji,jj) 
    979825                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
    … …  
    982828          ENDIF 
    983829       END_2D 
    984  
    985830 
    986831       WHERE ( lconv ) 
    … …  
    1011856       END_2D 
    1012857 
     858       DO_2D( 0, 0, 0, 0 ) 
     859        IF ( lpyc(ji,jj) ) THEN 
     860          IF ( j_ddh(ji,jj) == 0 ) THEN 
     861! Place holding code. Parametrization needs checking for these conditions. 
     862            zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     863            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
     864            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
     865          ELSE 
     866            zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     867            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
     868            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
     869          ENDIF 
     870          zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 
     871          zc_cubic = zuw_bse - zd_cubic 
     872! need ztau_sc_u to be available. Change to array.  
     873          DO jk = imld(ji,jj), ibld(ji,jj) 
     874             zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
     875             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     876          END DO 
     877          zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 
     878          zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 
     879          zc_cubic = zvw_bse - zd_cubic 
     880          DO jk = imld(ji,jj), ibld(ji,jj) 
     881            zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj) 
     882            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     883          END DO 
     884        ENDIF  ! lpyc 
     885       END_2D 
     886 
     887       IF(ln_dia_osm) THEN 
     888          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
     889          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     890       END IF 
    1013891! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    1014892 
    1015        WHERE ( lconv ) 
    1016           zsc_wth_1 = zwth0 
    1017           zsc_ws_1 = zws0 
    1018        ELSEWHERE 
    1019           zsc_wth_1 = 2.0 * zwthav 
    1020           zsc_ws_1 = zws0 
    1021        ENDWHERE 
     893       DO_2D( 1, 0, 1, 0 ) 
     894        
     895         IF ( lconv(ji,jj) ) THEN 
     896           zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 
     897           zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 
     898           IF ( lpyc(ji,jj) ) THEN 
     899! Pycnocline scales 
     900              zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 
     901              zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 
     902            ENDIF 
     903         ELSE 
     904           zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 
     905           zsc_ws_1(ji,jj) = zws0(ji,jj) 
     906         ENDIF 
     907       END_2D 
    1022908 
    1023909       DO_2D( 0, 0, 0, 0 ) 
    … …  
    1035921                    &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) ) 
    1036922            END DO 
     923! 
     924            IF ( lpyc(ji,jj) ) THEN 
     925! pycnocline 
     926              DO jk = imld(ji,jj), ibld(ji,jj) 
     927                zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
     928                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )  
     929                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )  
     930              END DO 
     931           ENDIF 
    1037932         ELSE 
    1038             DO jk = 2, ibld(ji,jj) 
    1039                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    1040                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    1041                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    1042             &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
    1043                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    1044             &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
    1045             END DO 
     933            IF( zdhdt(ji,jj) > 0. ) THEN 
     934              DO jk = 2, ibld(ji,jj) 
     935                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     936                 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     937                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     938              &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
     939                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     940               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
     941              END DO 
     942            ENDIF 
    1046943         ENDIF 
    1047944       END_2D 
    1048  
    1049945 
    1050946       WHERE ( lconv ) 
    … …  
    1090986          ENDIF 
    1091987       END_2D 
     988 
     989       IF(ln_dia_osm) THEN 
     990          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
     991          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
     992          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     993          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     994          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     995          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     996       END IF 
    1092997! 
    1093998! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
    1094999 
     1000 
     1001 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 
     1002 
    10951003      DO_2D( 0, 0, 0, 0 ) 
    1096          IF ( lconv(ji,jj) ) THEN 
     1004         IF ( .not. lconv(ji,jj) ) THEN 
    10971005            DO jk = 2, ibld(ji,jj) 
    1098                znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    1099                IF ( znd >= 0.0 ) THEN 
    1100                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    1101                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    1102                ELSE 
    1103                   ghamu(ji,jj,jk) = 0._wp 
    1104                   ghamv(ji,jj,jk) = 0._wp 
    1105                ENDIF 
    1106             END DO 
    1107          ELSE 
    1108             DO jk = 2, ibld(ji,jj) 
    1109                znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1006               znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 
    11101007               IF ( znd >= 0.0 ) THEN 
    11111008                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    … …  
    11201017 
    11211018      ! pynocline contributions 
    1122        ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
    1123        zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 
    11241019       DO_2D( 0, 0, 0, 0 ) 
    1125           DO jk= 2, ibld(ji,jj) 
    1126              znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    1127              ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1128              ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1129              ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1130              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) 
    1131              ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    1132           END DO 
     1020         IF ( .not. lconv(ji,jj) ) THEN 
     1021          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     1022             DO jk= 2, ibld(ji,jj) 
     1023                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1024                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     1025                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     1026                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     1027                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     1028             END DO 
     1029          END IF 
     1030         END IF 
    11331031       END_2D 
    1134  
    1135 ! Entrainment contribution. 
     1032      IF(ln_dia_osm) THEN 
     1033          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     1034          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     1035       END IF 
    11361036 
    11371037       DO_2D( 0, 0, 0, 0 ) 
    1138           IF ( lconv(ji,jj) ) THEN 
    1139             DO jk = 1, imld(ji,jj) - 1 
    1140                znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    1141                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1142                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1143                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1144                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
    1145             END DO 
    1146             DO jk = imld(ji,jj), ibld(ji,jj) 
    1147                znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    1148                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1149                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1150                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1151                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1152              END DO 
    1153           ENDIF 
    1154           ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
    1155           ghams(ji,jj,ibld(ji,jj)) = 0._wp 
    1156           ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
    1157           ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
     1038          ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     1039          ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     1040          ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     1041          ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    11581042       END_2D 
    11591043 
    1160  
     1044       IF(ln_dia_osm) THEN 
     1045          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     1046          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     1047          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
     1048          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
     1049          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
     1050       END IF 
    11611051       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    11621052       ! Need to put in code for contributions that are applied explicitly to 
    … …  
    11801070       IF(ln_dia_osm) THEN 
    11811071          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1072          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1073          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
    11821074       END IF 
    11831075 
    … …  
    12221114       END IF ! ln_convmix = .true. 
    12231115 
     1116 
     1117 
     1118       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
     1119          DO_2D( 0, 0, 0, 0 ) 
     1120              IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 
     1121             ! Calculate MLE flux contribution from surface fluxes 
     1122                DO jk = 1, ibld(ji,jj) 
     1123                  znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) 
     1124                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 
     1125                  ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 
     1126                 END DO 
     1127                 DO jk = 1, mld_prof(ji,jj) 
     1128                   znd = gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 
     1129                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 
     1130                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 
     1131                 END DO 
     1132         ! Viscosity for MLEs 
     1133                 DO jk = 1, mld_prof(ji,jj) 
     1134                   znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 
     1135                   zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 
     1136                 END DO 
     1137              ELSE 
     1138! Surface transports limited to OSBL.                  
     1139         ! Viscosity for MLEs 
     1140                 DO jk = 1, mld_prof(ji,jj) 
     1141                   znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 
     1142                   zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 
     1143                 END DO 
     1144              ENDIF 
     1145          END_2D 
     1146       ENDIF 
     1147 
     1148       IF(ln_dia_osm) THEN 
     1149          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1150          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1151          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
     1152       END IF 
     1153 
     1154 
    12241155       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
    1225        CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 
     1156       !CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 
    12261157 
    12271158       ! GN 25/8: need to change tmask --> wmask 
    … …  
    12441175            ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk) 
    12451176       END_3D 
     1177        ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged) 
     1178        CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 
    12461179        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged) 
    1247         ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    1248         CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp,   & 
    1249          &                  ghamu, 'U', 1.0_wp , ghamv, 'V', 1.0_wp ) 
    1250  
    1251        IF(ln_dia_osm) THEN 
     1180        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign changed) 
     1181        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W',  1.0_wp , ghams, 'W', 1.0_wp,   & 
     1182         &                            ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) 
     1183 
     1184      IF(ln_dia_osm) THEN 
    12521185         SELECT CASE (nn_osm_wave) 
    12531186         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1). 
    … …  
    12571190            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    12581191         ! Stokes drift read in from sbcwave  (=2). 
    1259          CASE(2) 
    1260             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    1261             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
     1192         CASE(2:3) 
     1193            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
     1194            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     1195            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period 
     1196            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height 
     1197            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 
     1198            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum 
     1199            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10 
    12621200            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 
    12631201                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
    … …  
    12701208         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0> 
    12711209         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth 
    1272          IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth 
     1210         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1211         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
     1212         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
     1213         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base 
     1214         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base 
     1215         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base 
     1216         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth 
     1217         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth 
    12731218         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth 
    12741219         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points 
    … …  
    12761221         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    12771222         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
     1223         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale 
     1224         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir # 
    12781225         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    12791226         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    12801227         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    12811228         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    1282          IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine 
     1229         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1230         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    12831231         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    1284          IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine 
    1285          IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine 
    1286          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML 
     1232         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
     1233         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
     1234         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
     1235         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
     1236         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
     1237 
     1238         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )               ! FK layer depth 
     1239         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )               ! FK target layer depth 
     1240         IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux 
     1241         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML 
     1242         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 
     1243         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )            ! FK dtdx at u-pt 
     1244         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )            ! FK dtdy at v-pt 
     1245         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )            ! FK dtdx at u-pt 
     1246         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )            ! FK dsdy at v-pt 
     1247         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )            ! FK dbdx at u-pt 
     1248         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )            ! FK dbdy at v-pt 
     1249         IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1250         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1251 
    12871252      END IF 
    1288       ! Lateral boundary conditions on p_avt  (sign unchanged) 
    1289       CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1.0_wp ) 
     1253 
     1254CONTAINS 
     1255! subroutine code changed, needs syntax checking. 
     1256  SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 
     1257 
     1258!!--------------------------------------------------------------------- 
     1259     !!                   ***  ROUTINE zdf_osm_diffusivity_viscosity  *** 
     1260     !! 
     1261     !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. 
     1262     !! 
     1263     !! ** Method  :  
     1264     !! 
     1265     !! !!---------------------------------------------------------------------- 
     1266     REAL(wp), DIMENSION(:,:,:) :: zdiffut 
     1267     REAL(wp), DIMENSION(:,:,:) :: zviscos 
     1268! local 
     1269 
     1270! Scales used to calculate eddy diffusivity and viscosity profiles 
     1271      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 
     1272      REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 
     1273      REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 
     1274      REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 
     1275! 
     1276      REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 
     1277       
     1278      REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 
     1279      REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 
     1280      REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 
     1281       
     1282      DO_2D( 0, 0, 0, 0 ) 
     1283          IF ( lconv(ji,jj) ) THEN 
     1284           
     1285            zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 
     1286            zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1287            zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 
     1288 
     1289            zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 
     1290            zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 
     1291 
     1292            IF ( lpyc(ji,jj) ) THEN 
     1293              zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 
     1294 
     1295              IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     1296                zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 
     1297              ENDIF 
     1298             
     1299              zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 
     1300              zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 
     1301              zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 
     1302               
     1303              zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 
     1304              zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 
     1305              IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     1306                zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 
     1307              ENDIF 
     1308              
     1309              zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 
     1310              zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 
     1311              zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 
     1312 
     1313              zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third 
     1314              zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 
     1315            ELSE 
     1316              zbeta_d_sc(ji,jj) = 1.0 
     1317              zbeta_v_sc(ji,jj) = 1.0 
     1318            ENDIF 
     1319          ELSE 
     1320            zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1321            zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1322          END IF 
     1323      END_2D 
     1324! 
     1325       DO_2D( 0, 0, 0, 0 ) 
     1326          IF ( lconv(ji,jj) ) THEN 
     1327             DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
     1328                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     1329                 ! 
     1330                 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
     1331                 ! 
     1332                 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 
     1333   &            *                                      ( 1.0 - 0.5 * zznd_ml**2 ) 
     1334             END DO 
     1335! pycnocline 
     1336             IF ( lpyc(ji,jj) ) THEN 
     1337! Diffusivity profile in the pycnocline given by cubic polynomial. 
     1338                za_cubic = 0.5 
     1339                zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
     1340                zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & 
     1341                     & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 
     1342                zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic  - zb_cubic ) 
     1343                zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1344                DO jk = imld(ji,jj) , ibld(ji,jj) 
     1345                  zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1346                      ! 
     1347                  zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 +   zd_cubic * zznd_pyc**3 ) 
     1348 
     1349                  zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) 
     1350                END DO 
     1351 ! viscosity profiles. 
     1352                za_cubic = 0.5 
     1353                zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
     1354                zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 
     1355                zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) 
     1356                zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1357                DO jk = imld(ji,jj) , ibld(ji,jj) 
     1358                   zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1359                    zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 
     1360                    zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) 
     1361                END DO 
     1362                IF ( zdhdt(ji,jj) > 0._wp ) THEN 
     1363                 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 
     1364                 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 
     1365                ELSE 
     1366                  zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     1367                  zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     1368                ENDIF 
     1369             ENDIF 
     1370          ELSE 
     1371          ! stable conditions 
     1372             DO jk = 2, ibld(ji,jj) 
     1373                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1374                zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
     1375                zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
     1376             END DO 
     1377 
     1378             IF ( zdhdt(ji,jj) > 0._wp ) THEN 
     1379                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) 
     1380                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) 
     1381             ENDIF 
     1382          ENDIF   ! end if ( lconv ) 
     1383          ! 
     1384       END_2D 
     1385        
     1386  END SUBROUTINE zdf_osm_diffusivity_viscosity 
     1387   
     1388  SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
     1389 
     1390!!--------------------------------------------------------------------- 
     1391     !!                   ***  ROUTINE zdf_osm_osbl_state  *** 
     1392     !! 
     1393     !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number 
     1394     !! 
     1395     !! ** Method  :  
     1396     !! 
     1397     !! !!---------------------------------------------------------------------- 
     1398 
     1399     INTEGER, DIMENSION(jpi,jpj) :: j_ddh  ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low. 
     1400      
     1401     LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear 
     1402 
     1403     REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 
     1404     REAL(wp), DIMENSION(jpi,jpj) :: zshear  ! production of TKE due to shear across the pycnocline 
     1405     REAL(wp), DIMENSION(jpi,jpj) :: zri_i  ! Interfacial Richardson Number 
     1406 
     1407! Local Variables 
     1408 
     1409     INTEGER :: jj, ji 
     1410      
     1411     REAL(wp), DIMENSION(jpi,jpj) :: zekman 
     1412     REAL(wp) :: zri_p, zri_b   ! Richardson numbers 
     1413     REAL(wp) :: zshear_u, zshear_v, zwb_shr 
     1414     REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1415 
     1416     REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 
     1417     REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 
     1418     REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04      
     1419     REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 
     1420     REAL, PARAMETER :: rn_ri_p_thresh = 27.0 
     1421     REAL, PARAMETER :: zrot=0._wp  ! dummy rotation rate of surface stress. 
     1422      
     1423! Determins stability and set flag lconv 
     1424     DO_2D( 0, 0, 0, 0 ) 
     1425      IF ( zhol(ji,jj) < 0._wp ) THEN 
     1426         lconv(ji,jj) = .TRUE. 
     1427       ELSE 
     1428          lconv(ji,jj) = .FALSE. 
     1429       ENDIF 
     1430     END_2D 
     1431  
     1432     zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 
     1433      
     1434     WHERE ( lconv ) 
     1435       zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 
     1436     END WHERE 
     1437 
     1438     zshear(:,:) = 0._wp 
     1439     j_ddh(:,:) = 1      
     1440  
     1441     DO_2D( 0, 0, 0, 0 ) 
     1442      IF ( lconv(ji,jj) ) THEN 
     1443         IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1444           zri_p = MAX (  SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) )  *  ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 
     1445                & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp ) 
     1446          
     1447           zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 
     1448                      
     1449           zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 
     1450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1451! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  ! 
     1452! full code available                                          ! 
     1453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1454           IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 
     1455! Growing shear layer 
     1456             j_ddh(ji,jj) = 0 
     1457             lshear(ji,jj) = .TRUE. 
     1458           ELSE 
     1459             j_ddh(ji,jj) = 1 
     1460             IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 
     1461! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 
     1462               lshear(ji,jj) = .TRUE. 
     1463             ELSE 
     1464! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 
     1465               zshear(ji,jj) = 0.5 * zshear(ji,jj) 
     1466               lshear(ji,jj) = .FALSE. 
     1467             ENDIF  
     1468           ENDIF                 
     1469         ELSE                ! zdb_bl test, note zshear set to zero 
     1470           j_ddh(ji,jj) = 2 
     1471           lshear(ji,jj) = .FALSE. 
     1472         ENDIF 
     1473       ENDIF 
     1474     END_2D 
     1475  
     1476! Calculate entrainment buoyancy flux due to surface fluxes. 
     1477 
     1478     DO_2D( 0, 0, 0, 0 ) 
     1479      IF ( lconv(ji,jj) ) THEN 
     1480        zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1481        zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1482        zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1483        zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
     1484        IF (nn_osm_SD_reduce > 0 ) THEN 
     1485        ! Effective Stokes drift already reduced from surface value 
     1486           zr_stokes = 1.0_wp 
     1487        ELSE 
     1488         ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 
     1489          ! requires further reduction where BL is deep 
     1490           zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1491         &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1492        END IF 
     1493        zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1494               &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1495               &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1496               &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1497          ! 
     1498      ENDIF 
     1499     END_2D 
     1500 
     1501     zwb_min(:,:) = 0._wp 
     1502 
     1503     DO_2D( 0, 0, 0, 0 ) 
     1504      IF ( lshear(ji,jj) ) THEN 
     1505        IF ( lconv(ji,jj) ) THEN 
     1506! Unstable OSBL 
     1507           zwb_shr = -za_wb_s * zshear(ji,jj) 
     1508           IF ( j_ddh(ji,jj) == 0 ) THEN 
     1509 
     1510! Developing shear layer, additional shear production possible. 
     1511 
     1512             zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) /  zhbl(ji,jj), 0._wp ) 
     1513             zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 
     1514             zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 
     1515              
     1516             zwb_shr = -za_wb_s * zshear(ji,jj) 
     1517              
     1518           ENDIF                 
     1519           zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
     1520           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1521        ELSE    ! IF ( lconv ) THEN - ENDIF 
     1522! Stable OSBL  - shear production not coded for first attempt.            
     1523        ENDIF  ! lconv 
     1524      ELSE  ! lshear 
     1525        IF ( lconv(ji,jj) ) THEN 
     1526! Unstable OSBL 
     1527           zwb_shr = -za_wb_s * zshear(ji,jj) 
     1528           zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
     1529           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1530        ENDIF  ! lconv 
     1531      ENDIF    ! lshear 
     1532     END_2D 
     1533   END SUBROUTINE zdf_osm_osbl_state 
     1534      
     1535      
     1536   SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1537     !!--------------------------------------------------------------------- 
     1538     !!                   ***  ROUTINE zdf_vertical_average  *** 
     1539     !! 
     1540     !! ** Purpose : Determines vertical averages from surface to jnlev. 
     1541     !! 
     1542     !! ** Method  : Averages are calculated from the surface to jnlev. 
     1543     !!              The external level used to calculate differences is ibld+ibld_ext 
     1544     !! 
     1545     !!---------------------------------------------------------------------- 
     1546 
     1547        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1548        INTEGER, DIMENSION(jpi,jpj) :: jp_ext 
     1549 
     1550        ! Alan: do we need zb? 
     1551        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb        ! Average temperature and salinity 
     1552        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
     1553        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
     1554        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
     1555 
     1556        INTEGER :: jk, ji, jj, ibld_ext 
     1557        REAL(wp) :: zthick, zthermal, zbeta 
     1558 
     1559 
     1560        zt   = 0._wp 
     1561        zs   = 0._wp 
     1562        zu   = 0._wp 
     1563        zv   = 0._wp 
     1564        DO_2D( 0, 0, 0, 0 ) 
     1565         zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1566         zbeta    = rab_n(ji,jj,1,jp_sal) 
     1567            ! average over depth of boundary layer 
     1568         zthick = epsln 
     1569         DO jk = 2, jnlev_av(ji,jj) 
     1570            zthick = zthick + e3t(ji,jj,jk,Kmm) 
     1571            zt(ji,jj)   = zt(ji,jj)  + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 
     1572            zs(ji,jj)   = zs(ji,jj)  + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
     1573            zu(ji,jj)   = zu(ji,jj)  + e3t(ji,jj,jk,Kmm) & 
     1574                  &            * ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) & 
     1575                  &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 
     1576            zv(ji,jj)   = zv(ji,jj)  + e3t(ji,jj,jk,Kmm) & 
     1577                  &            * ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) & 
     1578                  &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 
     1579         END DO 
     1580         zt(ji,jj) = zt(ji,jj) / zthick 
     1581         zs(ji,jj) = zs(ji,jj) / zthick 
     1582         zu(ji,jj) = zu(ji,jj) / zthick 
     1583         zv(ji,jj) = zv(ji,jj) / zthick 
     1584         zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 
     1585         ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 
     1586         IF ( ibld_ext < mbkt(ji,jj) ) THEN 
     1587           zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm) 
     1588           zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm) 
     1589           zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) & 
     1590                  &    / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) 
     1591           zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) & 
     1592                  &   / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 
     1593           zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1594         ELSE 
     1595           zdt(ji,jj) = 0._wp 
     1596           zds(ji,jj) = 0._wp 
     1597           zdu(ji,jj) = 0._wp 
     1598           zdv(ji,jj) = 0._wp 
     1599           zdb(ji,jj) = 0._wp 
     1600         ENDIF 
     1601        END_2D 
     1602   END SUBROUTINE zdf_osm_vertical_average 
     1603 
     1604   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 
     1605     !!--------------------------------------------------------------------- 
     1606     !!                   ***  ROUTINE zdf_velocity_rotation  *** 
     1607     !! 
     1608     !! ** Purpose : Rotates frame of reference of averaged velocity components. 
     1609     !! 
     1610     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w 
     1611     !! 
     1612     !!---------------------------------------------------------------------- 
     1613 
     1614        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle 
     1615        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current 
     1616        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline 
     1617 
     1618        INTEGER :: ji, jj 
     1619        REAL(wp) :: ztemp 
     1620 
     1621        DO_2D( 0, 0, 0, 0 ) 
     1622           ztemp = zu(ji,jj) 
     1623           zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 
     1624           zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1625           ztemp = zdu(ji,jj) 
     1626           zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
     1627           zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1628        END_2D 
     1629    END SUBROUTINE zdf_osm_velocity_rotation 
     1630 
     1631    SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
     1632     !!--------------------------------------------------------------------- 
     1633     !!                   ***  ROUTINE zdf_osm_osbl_state_fk  *** 
     1634     !! 
     1635     !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme. 
     1636     !!  lpyc :: determines whether pycnocline flux-grad relationship needs to be determined 
     1637     !!  lflux :: determines whether effects of surface flux extend below the base of the OSBL 
     1638     !!  lmle  :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.  
     1639     !! 
     1640     !! ** Method  :  
     1641     !! 
     1642     !!  
     1643     !!---------------------------------------------------------------------- 
     1644       
     1645! Outputs 
     1646      LOGICAL,  DIMENSION(jpi,jpj)  :: lpyc, lflux, lmle 
     1647      REAL(wp), DIMENSION(jpi,jpj)  :: zwb_fk 
     1648! 
     1649      REAL(wp), DIMENSION(jpi,jpj)  :: znd_param 
     1650      REAL(wp)                      :: zbuoy, ztmp, zpe_mle_layer 
     1651      REAL(wp)                      :: zpe_mle_ref, zwb_ent, zdbdz_mle_int 
     1652       
     1653      znd_param(:,:) = 0._wp 
     1654 
     1655        DO_2D( 0, 0, 0, 0 ) 
     1656          ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     1657          zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 
     1658        END_2D 
     1659        DO_2D( 0, 0, 0, 0 ) 
     1660                 ! 
     1661         IF ( lconv(ji,jj) ) THEN 
     1662           IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 
     1663             zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1664             zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1665             zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1666             zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1667! Calculate potential energies of actual profile and reference profile. 
     1668             zpe_mle_layer = 0._wp 
     1669             zpe_mle_ref = 0._wp 
     1670             DO jk = ibld(ji,jj), mld_prof(ji,jj) 
     1671               zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) 
     1672               zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
     1673               zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
     1674             END DO 
     1675! Non-dimensional parameter to diagnose the presence of thermocline 
     1676                 
     1677             znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) 
     1678           ENDIF 
     1679         ENDIF 
     1680        END_2D 
     1681 
     1682! Diagnosis 
     1683        DO_2D( 0, 0, 0, 0 ) 
     1684          IF ( lconv(ji,jj) ) THEN 
     1685            zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1686               &                  - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1687               &         + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 
     1688               &         - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1689            IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 
     1690              IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 
     1691! MLE layer growing 
     1692                IF ( znd_param (ji,jj) > 100. ) THEN 
     1693! Thermocline present 
     1694                  lflux(ji,jj) = .FALSE. 
     1695                  lmle(ji,jj) =.FALSE. 
     1696                ELSE 
     1697! Thermocline not present 
     1698                  lflux(ji,jj) = .TRUE. 
     1699                  lmle(ji,jj) = .TRUE. 
     1700                ENDIF  ! znd_param > 100 
     1701! 
     1702                IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
     1703                  lpyc(ji,jj) = .FALSE. 
     1704                ELSE 
     1705                   lpyc = .TRUE. 
     1706                ENDIF 
     1707              ELSE 
     1708! MLE layer restricted to OSBL or just below. 
     1709                IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
     1710! Weak stratification MLE layer can grow. 
     1711                  lpyc(ji,jj) = .FALSE. 
     1712                  lflux(ji,jj) = .TRUE. 
     1713                  lmle(ji,jj) = .TRUE. 
     1714                ELSE 
     1715! Strong stratification 
     1716                  lpyc(ji,jj) = .TRUE. 
     1717                  lflux(ji,jj) = .FALSE. 
     1718                  lmle(ji,jj) = .FALSE. 
     1719                ENDIF ! zdb_bl < rn_mle_thresh_bl and  
     1720              ENDIF  ! zhmle > 1.2 zhbl 
     1721            ELSE 
     1722              lpyc(ji,jj) = .TRUE. 
     1723              lflux(ji,jj) = .FALSE. 
     1724              lmle(ji,jj) = .FALSE. 
     1725              IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 
     1726            ENDIF !  -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5  
     1727          ELSE 
     1728! Stable Boundary Layer 
     1729            lpyc(ji,jj) = .FALSE. 
     1730            lflux(ji,jj) = .FALSE. 
     1731            lmle(ji,jj) = .FALSE. 
     1732          ENDIF  ! lconv 
     1733        END_2D 
     1734    END SUBROUTINE zdf_osm_osbl_state_fk 
     1735 
     1736    SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 
     1737     !!--------------------------------------------------------------------- 
     1738     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     1739     !! 
     1740     !! ** Purpose : Calculates the gradients below the OSBL 
     1741     !! 
     1742     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1743     !! 
     1744     !!---------------------------------------------------------------------- 
     1745 
     1746     INTEGER, DIMENSION(jpi,jpj)  :: jbase 
     1747     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
     1748 
     1749     INTEGER :: jj, ji, jkb, jkb1 
     1750     REAL(wp) :: zthermal, zbeta 
     1751 
     1752 
     1753     DO_2D( 0, 0, 0, 0 ) 
     1754        IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
     1755           zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1756           zbeta    = rab_n(ji,jj,1,jp_sal) 
     1757           jkb = jbase(ji,jj) 
     1758           jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     1759           zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & 
     1760                &   / e3t(ji,jj,ibld(ji,jj),Kmm) 
     1761           zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & 
     1762                &   / e3t(ji,jj,ibld(ji,jj),Kmm) 
     1763           zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1764        ELSE 
     1765           zdtdz(ji,jj) = 0._wp 
     1766           zdsdz(ji,jj) = 0._wp 
     1767           zdbdz(ji,jj) = 0._wp 
     1768        END IF 
     1769     END_2D 
     1770    END SUBROUTINE zdf_osm_external_gradients 
     1771 
     1772    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 
     1773 
     1774     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1775     REAL(wp), DIMENSION(jpi,jpj) :: zalpha 
     1776 
     1777     INTEGER :: jk, jj, ji 
     1778     REAL(wp) :: ztgrad, zsgrad, zbgrad 
     1779     REAL(wp) :: zgamma_b_nd, znd 
     1780     REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 
     1781     REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 
     1782 
     1783     DO_2D( 0, 0, 0, 0 ) 
     1784        IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     1785           IF ( lconv(ji,jj) ) THEN  ! convective conditions 
     1786             IF ( lpyc(ji,jj) ) THEN 
     1787                zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
     1788                zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 
     1789                zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 
     1790 
     1791                ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
     1792!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1793! Commented lines in this section are not needed in new code, once tested ! 
     1794! can be removed                                                          ! 
     1795!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1796!                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
     1797!                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
     1798                zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 
     1799                zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
     1800                DO jk = 2, ibld(ji,jj)+ibld_ext 
     1801                  znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 
     1802                  IF ( znd <= zzeta_m ) THEN 
     1803!                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 
     1804!                &        EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1805!                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 
     1806!                                  & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1807                     zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 
     1808                               & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1809                  ELSE 
     1810!                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1811!                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1812                      zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1813                  ENDIF 
     1814               END DO 
     1815            ENDIF ! if no pycnocline pycnocline gradients set to zero 
     1816           ELSE 
     1817              ! stable conditions 
     1818              ! if pycnocline profile only defined when depth steady of increasing. 
     1819              IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady. 
     1820                 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1821                    IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     1822                       ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 
     1823                       ztgrad = zdt_bl(ji,jj) * ztmp 
     1824                       zsgrad = zds_bl(ji,jj) * ztmp 
     1825                       zbgrad = zdb_bl(ji,jj) * ztmp 
     1826                       DO jk = 2, ibld(ji,jj) 
     1827                          znd = gdepw(ji,jj,jk,Kmm) * ztmp 
     1828                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1829                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1830                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1831                       END DO 
     1832                    ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1833                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
     1834                       ztgrad = zdt_bl(ji,jj) * ztmp 
     1835                       zsgrad = zds_bl(ji,jj) * ztmp 
     1836                       zbgrad = zdb_bl(ji,jj) * ztmp 
     1837                       DO jk = 2, ibld(ji,jj) 
     1838                          znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 
     1839                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1840                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1841                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1842                       END DO 
     1843                    ENDIF ! IF (zhol >=0.5) 
     1844                 ENDIF    ! IF (zdb_bl> 0.) 
     1845              ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1846           ENDIF          ! IF (lconv) 
     1847        ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 
     1848     END_2D 
     1849 
     1850    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
     1851 
     1852    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
     1853      !!--------------------------------------------------------------------- 
     1854      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
     1855      !! 
     1856      !! ** Purpose : Calculates velocity shear in the pycnocline 
     1857      !! 
     1858      !! ** Method  : 
     1859      !! 
     1860      !!---------------------------------------------------------------------- 
     1861 
     1862      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
     1863 
     1864      INTEGER :: jk, jj, ji 
     1865      REAL(wp) :: zugrad, zvgrad, znd 
     1866      REAL(wp) :: zzeta_v = 0.45 
    12901867      ! 
    1291    END SUBROUTINE zdf_osm 
    1292  
    1293  
    1294    SUBROUTINE zdf_osm_init( Kmm )  
     1868      DO_2D( 0, 0, 0, 0 ) 
     1869         ! 
     1870         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     1871            IF ( lconv (ji,jj) ) THEN 
     1872               ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
     1873!                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1874!                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1875!                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1876               !Alan is this right? 
     1877!                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1878!                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1879!                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1880!                       &      )/ (zdh(ji,jj)  + epsln ) 
     1881!                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     1882!                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     1883!                     IF ( znd <= 0.0 ) THEN 
     1884!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     1885!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     1886!                     ELSE 
     1887!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     1888!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     1889!                     ENDIF 
     1890!                  END DO 
     1891            ELSE 
     1892               ! stable conditions 
     1893               zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     1894               zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     1895               DO jk = 2, ibld(ji,jj) 
     1896                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1897                  IF ( znd < 1.0 ) THEN 
     1898                     zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     1899                  ELSE 
     1900                     zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     1901                  ENDIF 
     1902                  zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     1903               END DO 
     1904            ENDIF 
     1905            ! 
     1906         END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1907      END_2D 
     1908    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1909 
     1910   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     1911     !!--------------------------------------------------------------------- 
     1912     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     1913     !! 
     1914     !! ** Purpose : Calculates the rate at which hbl changes. 
     1915     !! 
     1916     !! ** Method  : 
     1917     !! 
     1918     !!---------------------------------------------------------------------- 
     1919 
     1920    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt        ! Rate of change of hbl 
     1921 
     1922    INTEGER :: jj, ji 
     1923    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 
     1924    REAL(wp) :: zvel_max!, zwb_min 
     1925    REAL(wp) :: zzeta_m = 0.3 
     1926    REAL(wp) :: zgamma_c = 2.0 
     1927    REAL(wp) :: zdhoh = 0.1 
     1928    REAL(wp) :: alpha_bc = 0.5 
     1929    REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 
     1930  
     1931  DO_2D( 0, 0, 0, 0 ) 
     1932     
     1933    IF ( lshear(ji,jj) ) THEN 
     1934       IF ( lconv(ji,jj) ) THEN    ! Convective 
     1935 
     1936          IF ( ln_osm_mle ) THEN 
     1937 
     1938             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1939       ! Fox-Kemper buoyancy flux average over OSBL 
     1940                zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     1941                     (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     1942             ELSE 
     1943                zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     1944             ENDIF 
     1945             zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1946             IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     1947                ! OSBL is deepening, entrainment > restratification 
     1948                IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
     1949! *** Used for shear Needs to be changed to work stabily 
     1950!                zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 
     1951!                zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 
     1952!                zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 
     1953!                za_1 = 1.0 / zgamma_b**2 - 0.017 
     1954!                za_2 = 1.0 / zgamma_b**3 - 0.0025 
     1955!                zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 
     1956                   zpsi = 0._wp 
     1957                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1958                   zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) 
     1959                   IF ( j_ddh(ji,jj) == 1 ) THEN 
     1960                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 
     1961                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     1962                     ELSE 
     1963                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     1964                     ENDIF 
     1965! Relaxation to dh_ref = zari * hbl 
     1966                     zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     1967                      
     1968                   ELSE  ! j_ddh == 0 
     1969! Growing shear layer 
     1970                     zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     1971                   ENDIF ! j_ddh 
     1972                     zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 
     1973                ELSE    ! zdb_bl >0 
     1974                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     1975                ENDIF 
     1976             ELSE   ! zwb_min + 2*zwb_fk_b < 0 
     1977                ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     1978                zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     1979 
     1980 
     1981             ENDIF 
     1982 
     1983          ELSE 
     1984             ! Fox-Kemper not used. 
     1985 
     1986               zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     1987               &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 
     1988               zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1989             ! added ajgn 23 July as temporay fix 
     1990 
     1991          ENDIF  ! ln_osm_mle 
     1992 
     1993         ELSE    ! lconv - Stable 
     1994             zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     1995             IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     1996                ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     1997                 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     1998             ELSE 
     1999                 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     2000             ENDIF 
     2001             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
     2002         ENDIF  ! lconv 
     2003    ELSE ! lshear 
     2004      IF ( lconv(ji,jj) ) THEN    ! Convective 
     2005 
     2006          IF ( ln_osm_mle ) THEN 
     2007 
     2008             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     2009       ! Fox-Kemper buoyancy flux average over OSBL 
     2010                zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     2011                     (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     2012             ELSE 
     2013                zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     2014             ENDIF 
     2015             zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2016             IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     2017                ! OSBL is deepening, entrainment > restratification 
     2018                IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
     2019                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     2020                ELSE 
     2021                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     2022                ENDIF 
     2023             ELSE 
     2024                ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     2025                zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     2026 
     2027 
     2028             ENDIF 
     2029 
     2030          ELSE 
     2031             ! Fox-Kemper not used. 
     2032 
     2033               zvel_max = -zwb_ent(ji,jj) / & 
     2034               &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 
     2035               zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     2036             ! added ajgn 23 July as temporay fix 
     2037 
     2038          ENDIF  ! ln_osm_mle 
     2039 
     2040         ELSE                        ! Stable 
     2041             zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     2042             IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     2043                ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     2044                 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 
     2045             ELSE 
     2046                 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     2047             ENDIF 
     2048             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
     2049         ENDIF  ! lconv 
     2050      ENDIF ! lshear 
     2051  END_2D 
     2052    END SUBROUTINE zdf_osm_calculate_dhdt 
     2053 
     2054    SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 
     2055     !!--------------------------------------------------------------------- 
     2056     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
     2057     !! 
     2058     !! ** Purpose : Increments hbl. 
     2059     !! 
     2060     !! ** Method  : If thechange in hbl exceeds one model level the change is 
     2061     !!              is calculated by moving down the grid, changing the buoyancy 
     2062     !!              jump. This is to ensure that the change in hbl does not 
     2063     !!              overshoot a stable layer. 
     2064     !! 
     2065     !!---------------------------------------------------------------------- 
     2066 
     2067 
     2068    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt   ! rates of change of hbl. 
     2069 
     2070    INTEGER :: jk, jj, ji, jm 
     2071    REAL(wp) :: zhbl_s, zvel_max, zdb 
     2072    REAL(wp) :: zthermal, zbeta 
     2073 
     2074     DO_2D( 0, 0, 0, 0 ) 
     2075        IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     2076! 
     2077! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     2078! 
     2079           zhbl_s = hbl(ji,jj) 
     2080           jm = imld(ji,jj) 
     2081           zthermal = rab_n(ji,jj,1,jp_tem) 
     2082           zbeta = rab_n(ji,jj,1,jp_sal) 
     2083 
     2084 
     2085           IF ( lconv(ji,jj) ) THEN 
     2086!unstable 
     2087 
     2088              IF( ln_osm_mle ) THEN 
     2089                 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2090              ELSE 
     2091 
     2092                 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     2093                   &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     2094 
     2095              ENDIF 
     2096 
     2097              DO jk = imld(ji,jj), ibld(ji,jj) 
     2098                 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 
     2099                      & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & 
     2100                      &  0.0 ) + zvel_max 
     2101 
     2102 
     2103                 IF ( ln_osm_mle ) THEN 
     2104                    zhbl_s = zhbl_s + MIN( & 
     2105                     & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     2106                     & e3w(ji,jj,jm,Kmm) ) 
     2107                 ELSE 
     2108                   zhbl_s = zhbl_s + MIN( & 
     2109                     & rn_Dt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     2110                     & e3w(ji,jj,jm,Kmm) ) 
     2111                 ENDIF 
     2112 
     2113!                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2114                 IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN 
     2115                   zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2116                   lpyc(ji,jj) = .FALSE. 
     2117                 ENDIF 
     2118                 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
     2119              END DO 
     2120              hbl(ji,jj) = zhbl_s 
     2121              ibld(ji,jj) = jm 
     2122          ELSE 
     2123! stable 
     2124              DO jk = imld(ji,jj), ibld(ji,jj) 
     2125                 zdb = MAX( & 
     2126                      & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )& 
     2127                      &           - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& 
     2128                      & 0.0 ) + & 
     2129          &       2.0 * zvstr(ji,jj)**2 / zhbl_s 
     2130 
     2131                 ! Alan is thuis right? I have simply changed hbli to hbl 
     2132                 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 
     2133                 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) ) ) * & 
     2134            &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 
     2135                 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 
     2136                 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) 
     2137 
     2138!                    zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2139                 IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 
     2140                   zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2141                   lpyc(ji,jj) = .FALSE. 
     2142                 ENDIF 
     2143                 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
     2144              END DO 
     2145          ENDIF   ! IF ( lconv ) 
     2146          hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4,Kmm) ) 
     2147          ibld(ji,jj) = MAX(jm, 4 ) 
     2148        ELSE 
     2149! change zero or one model level. 
     2150          hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) 
     2151        ENDIF 
     2152        zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
     2153     END_2D 
     2154 
     2155    END SUBROUTINE zdf_osm_timestep_hbl 
     2156 
     2157    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     2158      !!--------------------------------------------------------------------- 
     2159      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     2160      !! 
     2161      !! ** Purpose : Calculates thickness of the pycnocline 
     2162      !! 
     2163      !! ** Method  : The thickness is calculated from a prognostic equation 
     2164      !!              that relaxes the pycnocine thickness to a diagnostic 
     2165      !!              value. The time change is calculated assuming the 
     2166      !!              thickness relaxes exponentially. This is done to deal 
     2167      !!              with large timesteps. 
     2168      !! 
     2169      !!---------------------------------------------------------------------- 
     2170 
     2171      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
     2172       ! 
     2173      INTEGER :: jj, ji 
     2174      INTEGER :: inhml 
     2175      REAL(wp) :: zari, ztau, zdh_ref 
     2176      REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 
     2177 
     2178    DO_2D( 0, 0, 0, 0 ) 
     2179 
     2180      IF ( lshear(ji,jj) ) THEN 
     2181         IF ( lconv(ji,jj) ) THEN 
     2182           IF ( j_ddh(ji,jj) == 0 ) THEN 
     2183! ddhdt for pycnocline determined in osm_calculate_dhdt 
     2184             dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt 
     2185           ELSE 
     2186! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt  
     2187             IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 
     2188               zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2189             ELSE 
     2190               zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2191             ENDIF 
     2192             ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt ) 
     2193             dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) ) 
     2194             IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 
     2195           ENDIF 
     2196             
     2197         ELSE ! lconv 
     2198! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL  
     2199 
     2200            ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 
     2201            IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     2202               ! boundary layer deepening 
     2203               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     2204                  ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     2205                  zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     2206                       & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 ) 
     2207                  zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 
     2208               ELSE 
     2209                  zdh_ref = 0.2 * hbl(ji,jj) 
     2210               ENDIF 
     2211            ELSE     ! IF(dhdt < 0) 
     2212               zdh_ref = 0.2 * hbl(ji,jj) 
     2213            ENDIF    ! IF (dhdt >= 0) 
     2214            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 
     2215            IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse 
     2216            ! Alan: this hml is never defined or used -- do we need it? 
     2217         ENDIF 
     2218           
     2219      ELSE   ! lshear   
     2220! for lshear = .FALSE. calculate ddhdt here 
     2221 
     2222          IF ( lconv(ji,jj) ) THEN 
     2223 
     2224            IF( ln_osm_mle ) THEN 
     2225               IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
     2226                  ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
     2227                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 
     2228                     IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
     2229                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2230                     ELSE                                                     ! unstable 
     2231                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2232                     ENDIF 
     2233                     ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
     2234                     zdh_ref = zari * hbl(ji,jj) 
     2235                  ELSE 
     2236                     ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
     2237                     zdh_ref = 0.2 * hbl(ji,jj) 
     2238                  ENDIF 
     2239               ELSE 
     2240                  ztau = 0.2 * hbl(ji,jj) /  MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
     2241                  zdh_ref = 0.2 * hbl(ji,jj) 
     2242               ENDIF 
     2243            ELSE ! ln_osm_mle 
     2244               IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 
     2245                  IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
     2246                     zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2247                  ELSE                                                     ! unstable 
     2248                     zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2249                  ENDIF 
     2250                  ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
     2251                  zdh_ref = zari * hbl(ji,jj) 
     2252               ELSE 
     2253                  ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
     2254                  zdh_ref = 0.2 * hbl(ji,jj) 
     2255               ENDIF 
     2256 
     2257            END IF  ! ln_osm_mle 
     2258 
     2259            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 
     2260!               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 
     2261            IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 
     2262            ! Alan: this hml is never defined or used 
     2263         ELSE   ! IF (lconv) 
     2264            ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 
     2265            IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     2266               ! boundary layer deepening 
     2267               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     2268                  ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     2269                  zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     2270                       & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 ) 
     2271                  zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 
     2272               ELSE 
     2273                  zdh_ref = 0.2 * hbl(ji,jj) 
     2274               ENDIF 
     2275            ELSE     ! IF(dhdt < 0) 
     2276               zdh_ref = 0.2 * hbl(ji,jj) 
     2277            ENDIF    ! IF (dhdt >= 0) 
     2278            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 
     2279            IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse 
     2280         ENDIF       ! IF (lconv) 
     2281      ENDIF  ! lshear 
     2282  
     2283      hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     2284      inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj),Kmm), 1.e-3) ) , 1 ) 
     2285      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     2286      zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
     2287      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     2288    END_2D 
     2289 
     2290    END SUBROUTINE zdf_osm_pycnocline_thickness 
     2291 
     2292 
     2293   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
     2294      !!---------------------------------------------------------------------- 
     2295      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  *** 
     2296      !! 
     2297      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 
     2298      !! 
     2299      !! ** Method  : 
     2300      !! 
     2301      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     2302      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     2303 
     2304 
     2305      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 
     2306      REAL(wp), DIMENSION(jpi,jpj)     :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 
     2307      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
     2308      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy 
     2309 
     2310      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     2311      INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     2312      REAL(wp)                         :: zc 
     2313      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value 
     2314      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH 
     2315      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 
     2316      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv 
     2317!!---------------------------------------------------------------------- 
     2318      ! 
     2319      !                                      !==  MLD used for MLE  ==! 
     2320 
     2321      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point 
     2322      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
     2323      zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! convert density criteria into N^2 criteria 
     2324      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 
     2325         ikt = mbkt(ji,jj) 
     2326         zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 
     2327         IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     2328      END_3D 
     2329      DO_2D( 1, 1, 1, 1 ) 
     2330         mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 
     2331         zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
     2332      END_2D 
     2333      ! ensure mld_prof .ge. ibld 
     2334      ! 
     2335      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     2336      ! 
     2337      ztm(:,:) = 0._wp 
     2338      zsm(:,:) = 0._wp 
     2339      DO_3D( 1, 1, 1, 1, 1, ikmax ) 
     2340         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     2341         ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) 
     2342         zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) 
     2343      END_3D 
     2344      ! average temperature and salinity. 
     2345      ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 
     2346      zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 
     2347      ! calculate horizontal gradients at u & v points 
     2348 
     2349      DO_2D( 0, 0, 1, 0 ) 
     2350         zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     2351         zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     2352         zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 
     2353         ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 
     2354         ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 
     2355      END_2D 
     2356 
     2357      DO_2D( 1, 0, 0, 0 ) 
     2358         zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     2359         zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     2360         zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 
     2361         ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 
     2362         ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 
     2363      END_2D 
     2364 
     2365      CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm) 
     2366      CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm) 
     2367 
     2368      DO_2D( 0, 0, 1, 0 ) 
     2369         dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 
     2370      END_2D 
     2371      DO_2D( 1, 0, 0, 0 ) 
     2372         dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 
     2373      END_2D 
     2374 
     2375      DO_2D( 0, 0, 0, 0 ) 
     2376        ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     2377        zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
     2378             & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
     2379      END_2D 
     2380       
     2381 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
     2382  SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
     2383      !!---------------------------------------------------------------------- 
     2384      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
     2385      !! 
     2386      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 
     2387      !! 
     2388      !! ** Method  : 
     2389      !! 
     2390      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     2391      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     2392 
     2393      INTEGER, DIMENSION(jpi,jpj)      :: mld_prof 
     2394      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 
     2395      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     2396      INTEGER  ::   ii, ij, ik, jkb, jkb1  ! local integers 
     2397      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     2398      REAL(wp) ::  ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle 
     2399 
     2400   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
     2401 
     2402      DO_2D( 0, 0, 0, 0 ) 
     2403       IF ( lconv(ji,jj) ) THEN 
     2404          ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     2405   ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
     2406          zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
     2407          zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 
     2408       ENDIF 
     2409      END_2D 
     2410   ! Timestep mixed layer eddy depth. 
     2411      DO_2D( 0, 0, 0, 0 ) 
     2412        IF ( lmle(ji,jj) ) THEN  ! MLE layer growing. 
     2413! Buoyancy gradient at base of MLE layer. 
     2414           zthermal = rab_n(ji,jj,1,jp_tem) 
     2415           zbeta    = rab_n(ji,jj,1,jp_sal) 
     2416           jkb = mld_prof(ji,jj) 
     2417           jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     2418!               
     2419           zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) ) 
     2420           zdb_mle = zb_bl(ji,jj) - zbuoy  
     2421! Timestep hmle.  
     2422           hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / zdb_mle 
     2423        ELSE 
     2424           IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN 
     2425              hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau 
     2426           ELSE 
     2427              hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau 
     2428           ENDIF 
     2429        ENDIF 
     2430        hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) 
     2431       IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 
     2432      END_2D 
     2433 
     2434      mld_prof = 4 
     2435      DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 
     2436      IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     2437      END_3D 
     2438      DO_2D( 0, 0, 0, 0 ) 
     2439         zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm) 
     2440      END_2D 
     2441END SUBROUTINE zdf_osm_mle_parameters 
     2442 
     2443END SUBROUTINE zdf_osm 
     2444 
     2445 
     2446   SUBROUTINE zdf_osm_init( Kmm ) 
    12952447     !!---------------------------------------------------------------------- 
    12962448     !!                  ***  ROUTINE zdf_osm_init  *** 
    … …  
    13042456     !! ** input   :   Namlist namosm 
    13052457     !!---------------------------------------------------------------------- 
    1306      INTEGER, INTENT(in)    :: Kmm ! time level index (middle) 
    1307      ! 
     2458     INTEGER, INTENT(in)   ::   Kmm       ! time level 
    13082459     INTEGER  ::   ios            ! local integer 
    13092460     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     2461     REAL z1_t2 
    13102462     !! 
    13112463     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    1312           & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
    1313           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 
     2464          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 
     2465          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
     2466          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 
     2467! Namelist for Fox-Kemper parametrization. 
     2468      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
     2469           & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit 
     2470 
    13142471     !!---------------------------------------------------------------------- 
    13152472     ! 
    … …  
    13252482        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 
    13262483        WRITE(numout,*) '~~~~~~~~~~~~' 
    1327         WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters' 
    1328         WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la 
     2484        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters' 
     2485        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la 
     2486        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    13292487        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
     2488        WRITE(numout,*) '     Stokes drift reduction factor                 rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd 
    13302489        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    1331         WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes 
     2490        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
    13322491        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave 
    13332492        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave 
    … …  
    13392498        CASE(2) 
    13402499           WRITE(numout,*) '     calculated from ECMWF wave fields' 
     2500         END SELECT 
     2501        WRITE(numout,*) '     Stokes drift reduction                        nn_osm_SD_reduce', nn_osm_SD_reduce 
     2502        WRITE(numout,*) '     fraction of hbl to average SD over/fit' 
     2503        WRITE(numout,*) '     exponential with nn_osm_SD_reduce = 1 or 2    rn_osm_hblfrac =  ', rn_osm_hblfrac 
     2504        SELECT CASE (nn_osm_SD_reduce) 
     2505        CASE(0) 
     2506           WRITE(numout,*) '     No reduction' 
     2507        CASE(1) 
     2508           WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL' 
     2509        CASE(2) 
     2510           WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL' 
    13412511        END SELECT 
     2512        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 
    13422513        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm 
     2514        WRITE(numout,*) '         Threshold used to define BL              rn_osm_bl_thresh  = ', rn_osm_bl_thresh, 'm^2/s' 
    13432515        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix 
    13442516        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty 
    … …  
    13592531     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    13602532 
    1361      call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 
     2533 
     2534     IF( ln_osm_mle ) THEN 
     2535! Initialise Fox-Kemper parametrization 
     2536         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 
     2537903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 
     2538 
     2539         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 
     2540904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 
     2541         IF(lwm) WRITE ( numond, namosm_mle ) 
     2542 
     2543         IF(lwp) THEN                     ! Namelist print 
     2544            WRITE(numout,*) 
     2545            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 
     2546            WRITE(numout,*) '~~~~~~~~~~~~~' 
     2547            WRITE(numout,*) '   Namelist namosm_mle : ' 
     2548            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle 
     2549            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce 
     2550            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm' 
     2551            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's' 
     2552            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg' 
     2553            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c 
     2554            WRITE(numout,*) '         Threshold used to define MLE for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
     2555            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
     2556            WRITE(numout,*) '         switch to limit hmle                                      ln_osm_hmle_limit  = ', ln_osm_hmle_limit 
     2557            WRITE(numout,*) '         fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T.  rn_osm_hmle_limit = ', rn_osm_hmle_limit 
     2558         ENDIF         ! 
     2559     ENDIF 
     2560      ! 
     2561      IF(lwp) THEN 
     2562         WRITE(numout,*) 
     2563         IF( ln_osm_mle ) THEN 
     2564            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 
     2565            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation' 
     2566            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation' 
     2567         ELSE 
     2568            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 
     2569         ENDIF 
     2570      ENDIF 
     2571      ! 
     2572      IF( ln_osm_mle ) THEN                ! MLE initialisation 
     2573         ! 
     2574         rb_c = grav * rn_osm_mle_rho_c /rho0        ! Mixed Layer buoyancy criteria 
     2575         IF(lwp) WRITE(numout,*) 
     2576         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
     2577         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 
     2578         ! 
     2579         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            ! 
     2580! 
     2581         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation 
     2582            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  ) 
     2583            ! 
     2584         ENDIF 
     2585         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
     2586         z1_t2 = 2.e-5 
     2587         DO_2D( 1, 1, 1, 1 ) 
     2588            r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 
     2589         END_2D 
     2590         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 
     2591         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  ) 
     2592         ! 
     2593      ENDIF 
     2594 
     2595     call osm_rst( nit000, Kmm,  'READ' ) !* read or initialize hbl, dh, hmle 
     2596 
    13622597 
    13632598     IF( ln_zdfddm) THEN 
    … …  
    14542689     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    14552690 
    1456      INTEGER ::   id1, id2   ! iom enquiry index 
     2691     INTEGER ::   id1, id2, id3   ! iom enquiry index 
    14572692     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    14582693     INTEGER  ::   iiki, ikt ! local integer 
    … …  
    14602695     REAL(wp) ::   zN2_c           ! local scalar 
    14612696     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    1462      INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top) 
     2697     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 
    14632698     !!---------------------------------------------------------------------- 
    14642699     ! 
    … …  
    14702705        IF( id1 > 0 ) THEN                       ! 'wn' exists; read 
    14712706           CALL iom_get( numror, jpdom_auto, 'wn', ww ) 
    1472            WRITE(numout,*) ' ===>>>> :  ww read from restart file' 
     2707           WRITE(numout,*) ' ===>>>> :  wn read from restart file' 
    14732708        ELSE 
    14742709           ww(:,:,:) = 0._wp 
    1475            WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially' 
     2710           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    14762711        END IF 
     2712 
    14772713        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    1478         id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. ) 
     2714        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. ) 
    14792715        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return 
    14802716           CALL iom_get( numror, jpdom_auto, 'hbl' , hbl  ) 
    1481            CALL iom_get( numror, jpdom_auto, 'hbli', hbli  ) 
    1482            WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file' 
     2717           CALL iom_get( numror, jpdom_auto, 'dh', dh ) 
     2718           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file' 
     2719           IF( ln_osm_mle ) THEN 
     2720              id3 = iom_varid( numror, 'hmle'   , ldstop = .FALSE. ) 
     2721              IF( id3 > 0) THEN 
     2722                 CALL iom_get( numror, jpdom_auto, 'hmle' , hmle ) 
     2723                 WRITE(numout,*) ' ===>>>> :  hmle read from restart file' 
     2724              ELSE 
     2725                 WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl' 
     2726                 hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2727              END IF 
     2728           END IF 
    14832729           RETURN 
    1484         ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate 
     2730        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate 
    14852731           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 
    14862732        END IF 
    … …  
    14902736     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return 
    14912737     !!----------------------------------------------------------------------------- 
    1492      IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return 
    1493         IF( ntile /= 0 .AND. ntile /= nijtile ) RETURN        ! Do only on the last tile 
    1494  
     2738     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbl into the restart file, then return 
    14952739        IF(lwp) WRITE(numout,*) '---- osm-rst ----' 
    1496          CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww   ) 
    1497          CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl  ) 
    1498          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli ) 
     2740         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  ) 
     2741         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl ) 
     2742         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh  ) 
     2743         IF( ln_osm_mle ) THEN 
     2744            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle ) 
     2745         END IF 
    14992746        RETURN 
    15002747     END IF 
    … …  
    15042751     !!----------------------------------------------------------------------------- 
    15052752     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 
    1506      ALLOCATE( imld_rst(jpi,jpj) ) 
    15072753     ! w-level of the mixing and mixed layers 
    15082754     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) 
    … …  
    15132759     ! 
    15142760     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    1515      DO_3D( 1, 1, 1, 1, 1, jpkm1 )  ! Mixed layer level: w-level 
     2761     DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
    15162762        ikt = mbkt(ji,jj) 
    15172763        hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 
    … …  
    15202766     ! 
    15212767     DO_2D( 1, 1, 1, 1 ) 
    1522         iiki = imld_rst(ji,jj) 
    1523         hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth 
     2768        iiki = MAX(4,imld_rst(ji,jj)) 
     2769        hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm  )    ! Turbocline depth 
     2770        dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm  )     ! Turbocline depth 
    15242771     END_2D 
    1525      hbl = MAX(hbl,epsln) 
    1526      hbli(:,:) = hbl(:,:) 
    1527      DEALLOCATE( imld_rst ) 
     2772 
    15282773     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
     2774 
     2775     IF( ln_osm_mle ) THEN 
     2776        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2777        WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 
     2778     END IF 
     2779 
     2780     ww(:,:,:) = 0._wp 
     2781     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    15292782   END SUBROUTINE osm_rst 
    15302783 
    … …  
    15592812      ENDIF 
    15602813 
    1561       ! add non-local temperature and salinity flux 
    15622814      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    15632815         pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      & 
    … …  
    15692821      END_3D 
    15702822 
    1571  
    1572       ! save the non-local tracer flux trends for diagnostic 
     2823      ! save the non-local tracer flux trends for diagnostics 
    15732824      IF( l_trdtra )   THEN 
    15742825         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
    15752826         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
    1576 !!bug gm jpttdzdf ==> jpttosm 
    1577          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    1578          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     2827 
     2828         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) 
     2829         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) 
    15792830         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    15802831      ENDIF 
    … …  
    16422893 
    16432894   !!====================================================================== 
     2895 
    16442896END MODULE zdfosm 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfphy.F90

    r14050 r14051  
    179179      IF( ln_zdfmfc .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfosm' ) 
    180180      IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 
    181       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' ) 
     181      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 
    182182      IF( lk_top    .AND. ln_zdfmfc )   CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' ) 
    183183      IF(lwp) THEN 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/TOP/PISCES/SED/sedrst.F90

    r14018 r14051  
    9494            CALL iom_init( cw_sedrst_cxt, kdid = numrsw, ld_closedef = .FALSE. ) 
    9595#else 
    96                clinfo = 'Can not use XIOS in trc_rst_opn' 
    97                CALL ctl_stop(TRIM(clinfo)) 
     96            CALL ctl_stop( 'Can not use XIOS in trc_rst_opn' ) 
    9897#endif 
    99             ENDIF 
     98         ENDIF 
    10099 
    101100         lrst_sed = .TRUE. 

  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/TOP/trcrst.F90

    r14018 r14051  
    105105            CALL iom_init( cw_toprst_cxt, kdid = numrtw, ld_closedef = .FALSE. ) 
    106106#else 
    107                clinfo = 'Can not use XIOS in trc_rst_opn' 
    108                CALL ctl_stop(TRIM(clinfo)) 
     107            CALL ctl_stop( 'Can not use XIOS in trc_rst_opn' ) 
    109108#endif 
    110             ENDIF 
     109         ENDIF 
    111110         lrst_trc = .TRUE. 
    112111      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.