Changeset 12245


Ignore:
Timestamp:
2019-12-13T14:31:46+01:00 (10 months ago)
Author:
cguiavarch
Message:

Merge changes from /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser:12178-12226

Location:
NEMO/branches/UKMO/NEMO_4.0_FKOSM
Files:
10 edited
2 copied

Legend:

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

    r9799 r12245  
    1  bld::tool::fppkeys   key_c1d key_mpp_mpi key_iomput 
     1 bld::tool::fppkeys   key_c1d key_iomput 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/domain_def_nemo.xml

    r9930 r12245  
    1313       <zoom_domain ibegin="1" jbegin="1" ni="1" nj="1"/> 
    1414     </domain> 
    15      <!--   Eq section --> 
     15     <domain id="danger_T" domain_ref="grid_T" > 
     16       <zoom_domain id="danger_T"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     17     </domain> 
     18     <domain id="danger_U" domain_ref="grid_U" > 
     19       <zoom_domain id="danger_U"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     20     </domain> 
     21      <domain id="danger_V" domain_ref="grid_V" > 
     22       <zoom_domain id="danger_V"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     23     </domain> 
     24     <domain id="danger_W" domain_ref="grid_W" > 
     25       <zoom_domain id="danger_W"  ibegin="45" jbegin="196" ni="5" nj="5"/> 
     26    </domain> 
     27    <!--   Eq section --> 
    1628     <domain id="EqT" domain_ref="grid_T" > <zoom_domain id="EqT"/> </domain> 
    1729     <!--   TAO : see example above   --> 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml

    r10823 r12245  
    197197 
    198198      <field_group id="OSMOSIS_T" grid_ref="grid_T_2D"> 
     199        <field id="hml"                 long_name="mixed layr depth"                         unit="m"       /> 
     200        <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
     201        <field id="dh"                  long_name="Pycnocline thickness"                     unit=" m"      /> 
     202        <field id="ibld"                long_name="index of boundary layer depth"            unit="#"       /> 
     203        <field id="imld"                long_name="index of mixed layer depth"            unit="#"       /> 
     204        <field id="zhbl"                long_name="boundary layer depth -grid"                     unit="m"       /> 
     205        <field id="zhml"                long_name="mixed layer depth - grid"                        unit="m"       /> 
     206        <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
     207        <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
     208        <field id="us_x"        long_name="i component of active Stokes drift"                      unit="m/s"     /> 
     209        <field id="us_y"        long_name="j component of active Stokes drift"                      unit="m/s"     /> 
     210        <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    199211        <field id="zwth0"               long_name="surface non-local temperature flux"       unit="deg m/s" /> 
    200212        <field id="zws0"                long_name="surface non-local salinity flux"          unit="psu m/s" /> 
    201         <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
    202         <field id="hbli"                long_name="initial boundary layer depth"             unit="m"       /> 
    203         <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    204         <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
    205213        <field id="zwstrc"              long_name="convective velocity scale"                unit="m/s"     /> 
     214        <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    206215        <field id="zwstrl"              long_name="langmuir velocity scale"                  unit="m/s"     /> 
    207         <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    208         <field id="zhbl"                long_name="boundary layer depth"                     unit="m"       /> 
    209         <field id="zhml"                long_name="mixed layer depth"                        unit="m"       /> 
     216        <field id="zvstr"               long_name="mixed velocity scale"                     unit="m/s"     /> 
     217        <field id="zla"                 long_name="langmuir number"                          unit="m/s"     /> 
    210218        <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2"                        unit="mW"      /> 
    211219        <field id="wind_wave_power"     long_name="U_s \dot  tau"                            unit="mW"      /> 
    212220        <field id="wind_power"          long_name="\rho  u*^3"                               unit="mW"      /> 
    213221 
    214         <!-- extra OSMOSIS diagnostics --> 
     222       <!-- interior BL OSMOSIS diagnostics --> 
    215223        <field id="zwthav"              long_name="av turb flux of T in ml"                  unit="deg m/s" /> 
    216224        <field id="zt_ml"               long_name="av T in ml"                               unit="deg"     /> 
     225        <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
     226        <field id="zws_ent"            long_name="entrainment turb flux of S"                unit="10^-3 m/s" /> 
    217227        <field id="zwth_ent"            long_name="entrainment turb flux of T"               unit="deg m/s" /> 
    218         <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
    219         <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
    220       </field_group> 
    221  
    222       <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" operation="instant" > 
     228        <field id="zwb_ent"            long_name="entrainment turb flux of buoyancy"         unit="m^2/s^-3" /> 
     229  
     230        <field id="zdt_bl"             long_name="temperature jump at base of BL"                 unit="deg"      /> 
     231        <field id="zds_bl"             long_name="salinity jump at base of BL"                 unit="10^-3"      /> 
     232        <field id="zdb_bl"             long_name="buoyancy jump at base of BL"                 unit="m/s^2"      /> 
     233        <field id="zdu_bl"             long_name="u jump at base of BL"                       unit="m/s"      /> 
     234        <field id="zdv_bl"             long_name="v jump at base of BL"                       unit="m/s"      /> 
     235 
     236        <!-- extra OSMOSIS diagnostics for debugging --> 
     237       <field id="zsc_uw_1_0"       long_name="zsc u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 
     238        <field id="zsc_uw_1_f"       long_name="zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     239        <field id="zsc_vw_1_f"       long_name="zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     240        <field id="zsc_uw_2_f"       long_name="2nd zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     241        <field id="zsc_vw_2_f"       long_name="2nd zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     242        <field id="zuw_bse"       long_name="base u-flux T-points"                          unit="m^2/s^2" /> 
     243        <field id="zvw_bse"       long_name="base v-flux T-points"                          unit="m^2/s^2" /> 
     244 
     245       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     246         <field id="hmle"          long_name="OBL FK-layer thickness"                                     unit="m"        /> 
     247        <field id="mld_prof"              long_name="FK-layer depth index"                  unit="#" /> 
     248        <field id="zmld"          long_name="target FK-layer thickness"                                     unit="m"        /> 
     249        <field id="zwb_fk"          long_name="FK b-flux"                                     unit="m^2 s^-3"        /> 
     250        <field id="zwb_fk_b"          long_name="layer averaged FK b-flux"                 unit="m^2 s^-3"       /> 
     251        <field id="zdiff_mle"          long_name="max FK diffusivity in MLE"       unit=" 10^-4 m^2 s^-1"       /> 
     252        <field id="zvel_mle"          long_name="FK velocity scale in MLE"       unit=" m s^-1"       /> 
     253    </field_group> 
     254 
     255      <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" > 
     256        <field id="zviscos"       long_name="BL viscosity"   unit="m^2/s" /> 
    223257        <field id="ghamt"       long_name="non-local temperature flux"                       unit="deg m/s" /> 
    224258        <field id="ghams"       long_name="non-local salinity flux"                          unit="psu m/s" /> 
    225259        <field id="zdtdz_pyc"   long_name="Pycnocline temperature gradient"                  unit=" deg/m"  /> 
    226       </field_group> 
     260        <field id="zdsdz_pyc"   long_name="Pycnocline salinity gradient"                  unit=" 10^-3/m"  /> 
     261        <field id="zdbdz_pyc"   long_name="Pycnocline buoyancy gradient"                  unit=" s^-2"  /> 
     262        <field id="zdudz_pyc"   long_name="Pycnocline u gradient"                  unit=" s^-2"  /> 
     263        <field id="zdvdz_pyc"   long_name="Pycnocline v gradient"                  unit=" s^-2"  /> 
     264 
     265        <!-- extra OSMOSIS diagnostics for debugging --> 
     266         <field id="ghamu_00"       long_name="initial non-local u-momentum flux"   unit="m^2/s^2" /> 
     267        <field id="ghamv_00"       long_name="initial non-local v-momentum flux"   unit="m^2/s^2" /> 
     268        <field id="ghamu_0"       long_name="after dstokes non-local u-momentum flux"   unit="m^2/s^2" /> 
     269        <field id="ghamu_f"       long_name="after Coriolis non-local u-momentum flux"   unit="m^2/s^2" /> 
     270        <field id="ghamv_f"       long_name="after Coriolis  non-local v-momentum flux"   unit="m^2/s^2" /> 
     271        <field id="ghamu_b"       long_name="after buoyancy added non-local u-momentum flux"   unit="m^2/s^2" /> 
     272        <field id="ghamv_b"       long_name="after buoyancy added  non-local v-momentum flux"  unit="m^2/s^2" /> 
     273        <field id="ghamu_1"       long_name="after entrainment non-local u-momentum flux"   unit="m^2/s^2" /> 
     274        <field id="ghamv_1"       long_name="after entrainment  non-local v-momentum flux"  unit="m^2/s^2" /> 
     275     </field_group> 
    227276 
    228277      <field_group id="OSMOSIS_U" grid_ref="grid_U_2D" > 
    229278        <field id="ghamu"       long_name="non-local u-momentum flux"   grid_ref="grid_U_3D" unit="m^2/s^2" /> 
    230         <field id="us_x"        long_name="i component of Stokes drift"                      unit="m/s"     /> 
    231       </field_group> 
     279       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     280       <field id="zdtdx"          long_name="FK  T x-gradient"                                     unit=" deg C m^-1"        /> 
     281        <field id="zdsdx"          long_name="FK  S x-gradient"                                     unit=" 10^-3 m^-1"        /> 
     282        <field id="dbdx_mle"          long_name="FK  B x-gradient"                                     unit=" s^-2"        /> 
     283     </field_group> 
    232284 
    233285      <field_group id="OSMOSIS_V" grid_ref="grid_V_2D" > 
    234286        <field id="ghamv"       long_name="non-local v-momentum flux"   grid_ref="grid_V_3D" unit="m^2/s^2" /> 
    235         <field id="us_y"        long_name="j component of Stokes drift"                      unit="m/s"     /> 
     287        <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     288        <field id="zdtdy"          long_name="FK T y-gradient"                                     unit=" deg C m^-1"        /> 
     289        <field id="zdsdy"          long_name="FK S y-gradient"                                     unit=" 10^-3 m^-1"        /> 
     290        <field id="dbdy_mle"          long_name="FK B y-gradient"                                     unit=" s^-2"        /> 
    236291      </field_group> 
    237292       
     
    679734     <field id="strd_zdfp"     long_name="salinity   -trend: pure vert. diffusion"   unit="1e-3/s" /> 
    680735 
    681      <!-- --> 
     736     <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     737     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s" /> 
     738     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s" /> 
     739 
     740 
     741    <!-- --> 
    682742     <field id="ttrd_dmp"      long_name="temperature-trend: interior restoring"        unit="degC/s" /> 
    683743     <field id="strd_dmp"      long_name="salinity   -trend: interior restoring"        unit="1e-3/s" /> 
     
    715775     <field id="strd_zdfp_e3t"     unit="1e-3/s * m"  >  strd_zdfp * e3t </field> 
    716776 
     777          <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     778     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field> 
     779     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field> 
     780      
    717781     <!-- --> 
    718782     <field id="ttrd_dmp_e3t"      unit="degC/s * m"  >  ttrd_dmp * e3t </field> 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/grid_def_nemo.xml

    r10226 r12245  
    4848         <axis id="depthw" /> 
    4949       </grid> 
    50         <!--  --> 
     50       <!--  --> 
    5151       <grid id="grid_1point" > 
    5252         <domain domain_ref="1point"/> 
     
    5757         <axis id="nfloat" /> 
    5858       </grid> 
     59      <grid id="danger_T_2D" > 
     60         <domain id="danger_T" /> 
     61       </grid> 
     62        <!--  --> 
     63       <grid id="danger_T_3D" > 
     64         <domain id="danger_T" /> 
     65         <axis id="deptht" /> 
     66       </grid> 
     67       <!--  --> 
     68      <grid id="danger_U_2D" > 
     69         <domain id="danger_U" /> 
     70       </grid> 
     71        <!--  --> 
     72       <grid id="danger_U_3D" > 
     73         <domain id="danger_U" /> 
     74         <axis id="depthu" /> 
     75       </grid> 
     76      <!--  --> 
     77      <grid id="danger_V_2D" > 
     78         <domain id="danger_V" /> 
     79       </grid> 
     80        <!--  --> 
     81       <grid id="danger_V_3D" > 
     82         <domain id="danger_V" /> 
     83         <axis id="depthv" /> 
     84       </grid> 
     85      <!--  --> 
     86      <grid id="danger_W_2D" > 
     87         <domain id="danger_W" /> 
     88       </grid> 
     89        <!--  --> 
     90       <grid id="danger_W_3D" > 
     91         <domain id="danger_W" /> 
     92         <axis id="depthw" /> 
     93       </grid> 
     94        <!--  --> 
     95     <grid id="danger_T_2D" > 
     96         <domain id="danger_T" /> 
     97       </grid> 
     98        <!--  --> 
     99       <grid id="danger_T_3D" > 
     100         <domain id="danger_T" /> 
     101         <axis axis_ref="deptht" /> 
     102       </grid> 
     103       <!--  --> 
     104      <grid id="danger_U_2D" > 
     105         <domain id="danger_U" /> 
     106       </grid> 
     107        <!--  --> 
     108       <grid id="danger_U_3D" > 
     109         <domain id="danger_U" /> 
     110         <axis axis_ref="depthu" /> 
     111       </grid> 
     112      <!--  --> 
     113      <grid id="danger_V_2D" > 
     114         <domain id="danger_V" /> 
     115       </grid> 
     116        <!--  --> 
     117       <grid id="danger_V_3D" > 
     118         <domain id="danger_V" /> 
     119         <axis axis_ref="depthv" /> 
     120       </grid> 
     121      <!--  --> 
     122      <grid id="danger_W_2D" > 
     123         <domain id="danger_W" /> 
     124       </grid> 
     125        <!--  --> 
     126       <grid id="danger_W_3D" > 
     127         <domain id="danger_W" /> 
     128         <axis axis_ref="depthw" /> 
     129       </grid> 
     130    </grid_definition> 
    59131 
    60     </grid_definition> 
    61      
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/namelist_ref

    r10888 r12245  
    10701070&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
    10711071!----------------------------------------------------------------------- 
    1072    ln_use_osm_la = .false.      !  Use namelist  rn_osm_la 
     1072   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    10731073   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    10741074   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     
    10851085      !                        !  = 1: Pierson Moskowitz wave spectrum 
    10861086      !                        !  = 0: Constant La# = 0.3 
     1087   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
     1088/ 
     1089!----------------------------------------------------------------------- 
     1090&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     1091!----------------------------------------------------------------------- 
     1092   rn_osm_mle_ce       = 0.06      ! magnitude of the MLE (typical value: 0.06 to 0.08) 
     1093   nn_osm_mle          = 0         ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 
     1094   rn_osm_mle_lf       = 5.e+3     ! typical scale of mixed layer front (meters)                      (case rn_osm_mle=0) 
     1095   rn_osm_mle_time     = 172800.   ! time scale for mixing momentum across the mixed layer (seconds)  (case rn_osm_mle=0) 
     1096   rn_osm_mle_lat      = 20.       ! reference latitude (degrees) of MLE coef.                        (case rn_mle=1) 
     1097   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
     1098   rn_osm_mle_thresh = 0.0005      ! delta b criterion used for FK MLD criterion 
     1099   rn_osm_mle_tau     = 172800.   ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
    10871100/ 
    10881101!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/ref_cfgs.txt

    r10888 r12245  
    99ORCA2_ICE_PISCES OCE TOP ICE NST 
    1010SPITZ12 OCE ICE 
     11eORCA1 OCE ICE 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRA/tramle.F90

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

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

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

    r10888 r12245  
    2525   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 
    2626   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 
    27    !!            (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out). 
     27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 
    2828   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 
    2929   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 
    3030   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 
     31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, 
     32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers 
     33   !!                  (b) The stable OSBL depth parametrization changed. 
     34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. 
     35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 
    3136   !!---------------------------------------------------------------------- 
    3237 
     
    4045   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD) 
    4146   !!   dyn_osm       : compute and add to u & v trensd the non-local flux 
     47   !! 
     48   !! Subroutines in revised code. 
    4249   !!---------------------------------------------------------------------- 
    4350   USE oce            ! ocean dynamics and active tracers 
     
    6976   PUBLIC   tra_osm       ! routine called by step.F90 
    7077   PUBLIC   trc_osm       ! routine called by trcstp.F90 
    71    PUBLIC   dyn_osm       ! routine called by 'step.F90' 
     78   PUBLIC   dyn_osm       ! routine called by step.F90 
     79 
     80   PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90 
    7281 
    7382   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux 
     
    7786   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt 
    7887   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth 
    79    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbli     !: intial boundary layer depth for stable blayer 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline 
     89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth 
    8090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift. 
     91 
     92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts 
     93   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML 
     96   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer. 
    8197 
    8298   !                      !!** Namelist  namzdf_osm  ** 
    8399   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la 
     100 
     101   LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 
     102 
    84103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    85104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
     
    96115   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s) 
    97116 
     117! OSMOSIS mixed layer eddy parametrization constants 
     118   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt 
     119   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient 
     120   !                                        ! parameters used in nn_osm_mle = 0 case 
     121   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front 
     122   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer 
     123   !                                        ! parameters used in nn_osm_mle = 1 case 
     124   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
     125   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
     126   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
     127   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 
     128   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
     129   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
     130   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
     131 
     132 
    98133   !                                    !!! ** General constants  ** 
    99    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number 
     134   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
     135   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    100136   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    101137   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
     
    114150      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    115151      !!---------------------------------------------------------------------- 
    116      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 
    117           &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     & 
    118           &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc ) 
     152     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
     153          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     154          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     155 
     156     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 
     157          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 
     158 
     159!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ? 
     160!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     161!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     162!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
     163!     IF ( ln_osm_mle ) THEN 
     164!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 
     165!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 
     166!     ENDIF 
     167 
    119168     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
    120169     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     
    161210      !! 
    162211      INTEGER ::   ji, jj, jk                   ! dummy loop indices 
     212 
     213      INTEGER ::   jl                   ! dummy loop indices 
     214 
    163215      INTEGER ::   ikbot, jkmax, jkm1, jkp2     ! 
    164216 
     
    166218      REAL(wp) ::   zbeta, zthermal                                  ! 
    167219      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 
    168       REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       ! 
     220      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      ! 
     221 
    169222      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density 
    170223      INTEGER  ::   jm                          ! dummy loop indices 
     
    191244      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    192245      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     246 
     247 
     248      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL 
     249      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux 
     250      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 
     251      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK 
     252 
    193253      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift 
    194254      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number 
     
    196256      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 
    197257      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    198       LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 
     258      LOGICAL, DIMENSION(jpi,jpj)  :: lconv    ! unstable/stable bl 
    199259 
    200260      ! mixed-layer variables 
     
    208268      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    209269      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid 
     270 
     271      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 
     272      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid 
     273 
    210274      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    211275      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
     276      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure. 
     277      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients 
     278 
     279      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
     280 
    212281      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    213282      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
     
    238307      ! Temporary variables 
    239308      INTEGER :: inhml 
    240       INTEGER :: i_lconv_alloc 
    241309      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 
    242310      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables 
     
    248316      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    249317 
     318      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
     319      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 
     320      REAL(wp) :: zzeta_s = 0._wp 
     321      REAL(wp) :: zzeta_v = 0.46 
     322      REAL(wp) :: zabsstke 
     323 
    250324      ! For debugging 
    251325      INTEGER :: ikt 
    252326      !!-------------------------------------------------------------------- 
    253327      ! 
    254       ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc ) 
    255       IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv') 
    256  
    257328      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    258329      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
     
    268339      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    269340      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     341 
    270342      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    271343      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     
    277349      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    278350      ! 
     351      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 
     352 
     353      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
     354         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 
     355         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 
     356         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 
     357         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp 
     358      ENDIF 
     359      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
     360 
    279361      ! Flux-Gradient arrays. 
    280362      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
     
    287369      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    288370 
     371      zdhdt_2(:,:) = 0._wp 
    289372      ! hbl = MAX(hbl,epsln) 
    290373      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    350433              ! Use wind speed wndm included in sbc_oce module 
    351434              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    352               dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
     435              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    353436           END DO 
    354437        END DO 
     
    362445              ! It could represent the effects of the spread of wave directions 
    363446              ! around the mean wind. The effect of this adjustment needs to be tested. 
    364               zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
    365                    &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
    366               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
     447              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     448              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     449              dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
    367450           END DO 
    368451        END DO 
     
    375458           ! Langmuir velocity scale (zwstrl), at T-point 
    376459           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
    377            ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
    378            ! Intended as a possible test. Doesn't affect LES results for entrainment, 
    379            !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
    380            zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
    381                 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
    382                 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
    383            ! define La this way so effects of Stokes penetration depth on velocity scale are included 
    384            zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
     460           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 
     461           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 
    385462           ! Velocity scale that tends to zustar for large Langmuir numbers 
    386463           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
     
    389466           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    390467           ! Note zustke and zwstrl are not amended. 
    391            IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    392468           ! 
    393469           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
     
    406482     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 
    407483     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408      ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
    410       ibld(:,:) = 3 
    411       DO jk = 4, jpkm1 
     484     ! BL must be always 4 levels deep. 
     485      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
     486      ibld(:,:) = 4 
     487      DO jk = 5, jpkm1 
    412488         DO jj = 2, jpjm1 
    413489            DO ji = 2, jpim1 
     
    419495      END DO 
    420496 
    421       DO jj = 2, jpjm1                                 !  Vertical slab 
     497      DO jj = 2, jpjm1 
    422498         DO ji = 2, jpim1 
    423                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    424                zbeta    = rab_n(ji,jj,1,jp_sal) 
    425                zt   = 0._wp 
    426                zs   = 0._wp 
    427                zu   = 0._wp 
    428                zv   = 0._wp 
    429                ! average over depth of boundary layer 
    430                zthick=0._wp 
    431                DO jm = 2, ibld(ji,jj) 
    432                   zthick=zthick+e3t_n(ji,jj,jm) 
    433                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    434                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    435                   zu   = zu  + e3t_n(ji,jj,jm) & 
    436                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    437                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    438                   zv   = zv  + e3t_n(ji,jj,jm) & 
    439                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    440                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    441                END DO 
    442                zt_bl(ji,jj) = zt / zthick 
    443                zs_bl(ji,jj) = zs / zthick 
    444                zu_bl(ji,jj) = zu / zthick 
    445                zv_bl(ji,jj) = zv / zthick 
    446                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    447                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    448                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    449                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    450                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    451                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    452                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    453                IF ( lconv(ji,jj) ) THEN    ! Convective 
    454                       zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    455                            &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    456  
    457                       zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    458                            &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    460 !                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    461 !                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    462  
    463 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    464 !                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    465                       zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    466                ELSE                        ! Stable 
    467                       zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
    468                            &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
    469                            & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
    470                            &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
    471                       zzdhdt = zzdhdt + zwbav(ji,jj) 
    472                       IF ( zzdhdt < 0._wp ) THEN 
    473                       ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    474                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    475                       ELSE 
    476                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    477                               &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    478                       ENDIF 
    479                       zzdhdt = 2.0 * zzdhdt / zpert 
    480                ENDIF 
    481                zdhdt(ji,jj) = zzdhdt 
    482            END DO 
     499            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     500            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 
     501            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     502         END DO 
    483503      END DO 
    484  
     504      ! Averages over well-mixed and boundary layer 
     505      ! Alan: do we need zb_nl?, zb_ml? 
     506      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     507      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     508! External gradient 
     509      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     510 
     511 
     512      IF ( ln_osm_mle ) THEN 
     513         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     514         CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     515       ENDIF 
     516 
     517! Rate of change of hbl 
     518      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
    485519      ! Calculate averages over depth of boundary layer 
     520         DO jj = 2, jpjm1 
     521            DO ji = 2, jpim1 
     522               zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 
     523               ! adjustment to represent limiting by ocean bottom 
     524               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
     525            END DO 
     526         END DO 
     527 
    486528      imld = ibld           ! use imld to hold previous blayer index 
    487       ibld(:,:) = 3 
    488  
    489       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
    490       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    491       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     529      ibld(:,:) = 4 
    492530 
    493531      DO jk = 4, jpkm1 
     
    495533            DO ji = 2, jpim1 
    496534               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    497                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     535                  ibld(ji,jj) = jk 
    498536               ENDIF 
    499537            END DO 
     
    504542! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    505543! 
    506       DO jj = 2, jpjm1 
    507          DO ji = 2, jpim1 
    508             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     544      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     545      ! Alan: do we need zb_ml? 
     546      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    509547! 
    510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    511548! 
    512                zhbl_s = hbl(ji,jj) 
    513                jm = imld(ji,jj) 
    514                zthermal = rab_n(ji,jj,1,jp_tem) 
    515                zbeta = rab_n(ji,jj,1,jp_sal) 
    516                IF ( lconv(ji,jj) ) THEN 
    517 !unstable 
    518                   zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    519                        &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    520  
    521                   DO jk = imld(ji,jj), ibld(ji,jj) 
    522                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
    523                           & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524  
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
    526                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527  
    528                      IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
    529                   END DO 
    530                   hbl(ji,jj) = zhbl_s 
    531                   ibld(ji,jj) = jm 
    532                   hbli(ji,jj) = hbl(ji,jj) 
    533                ELSE 
    534 ! stable 
    535                   DO jk = imld(ji,jj), ibld(ji,jj) 
    536                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          & 
    537                           &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 
    538                           & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    539  
    540                      zhbl_s = zhbl_s +  (                                                                                & 
    541                           &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
    542                           &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
    543                           &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    544                           &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    545                           &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here 
    546  
    547                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    548                      IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
    549                   END DO 
    550                   hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 
    551                   ibld(ji,jj) = MAX(jm, 3 ) 
    552                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    553                ENDIF   ! IF ( lconv ) 
    554             ELSE 
    555 ! change zero or one model level. 
    556                hbl(ji,jj) = zhbl_t(ji,jj) 
    557                IF ( lconv(ji,jj) ) THEN 
    558                   hbli(ji,jj) = hbl(ji,jj) 
    559                ELSE 
    560                   hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 
    561                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    562                ENDIF 
    563             ENDIF 
    564             zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    565          END DO 
    566       END DO 
     549      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
    567550      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    568  
    569 ! Recalculate averages over boundary layer after depth updated 
    570      ! Consider later  combining this into the loop above and looking for columns 
    571      ! where the index for base of the boundary layer have changed 
    572       DO jj = 2, jpjm1                                 !  Vertical slab 
    573          DO ji = 2, jpim1 
    574                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    575                zbeta    = rab_n(ji,jj,1,jp_sal) 
    576                zt   = 0._wp 
    577                zs   = 0._wp 
    578                zu   = 0._wp 
    579                zv   = 0._wp 
    580                ! average over depth of boundary layer 
    581                zthick=0._wp 
    582                DO jm = 2, ibld(ji,jj) 
    583                   zthick=zthick+e3t_n(ji,jj,jm) 
    584                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    585                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    586                   zu   = zu  + e3t_n(ji,jj,jm) & 
    587                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    588                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    589                   zv   = zv  + e3t_n(ji,jj,jm) & 
    590                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    591                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    592                END DO 
    593                zt_bl(ji,jj) = zt / zthick 
    594                zs_bl(ji,jj) = zs / zthick 
    595                zu_bl(ji,jj) = zu / zthick 
    596                zv_bl(ji,jj) = zv / zthick 
    597                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    598                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    599                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    600                       &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    601                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    602                       &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    603                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    604                zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    605                IF ( lconv(ji,jj) ) THEN 
    606                   IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    607                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    608                            zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
    609                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    610                      ELSE                                                     ! unstable 
    611                            zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
    612                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    613                      ENDIF 
    614                      IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
    615                         zari = 0.2 
    616                         zwb_ent(ji,jj) = 0._wp 
    617                      ENDIF 
    618                      inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    619                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    620                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    621                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    622                   ELSE  ! IF (zdb_bl) 
    623                      imld(ji,jj) = ibld(ji,jj) - 1 
    624                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    625                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    626                   ENDIF 
    627                ELSE   ! IF (lconv) 
    628                   IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    629                   ! boundary layer deepening 
    630                      IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    631                   ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    632                         zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    633                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    634                         inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    635                         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    636                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    637                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    638                      ELSE 
    639                         imld(ji,jj) = ibld(ji,jj) - 1 
    640                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    641                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    642                      ENDIF ! IF (zdb_bl > 0.0) 
    643                   ELSE     ! IF(dhdt >= 0) 
    644                   ! boundary layer collapsing. 
    645                      imld(ji,jj) = ibld(ji,jj) 
    646                      zhml(ji,jj) = zhbl(ji,jj) 
    647                      zdh(ji,jj) = 0._wp 
    648                   ENDIF    ! IF (dhdt >= 0) 
    649                ENDIF       ! IF (lconv) 
    650          END DO 
    651       END DO 
    652  
    653       ! Average over the depth of the mixed layer in the convective boundary layer 
    654       ! Also calculate entrainment fluxes for temperature and salinity 
    655       DO jj = 2, jpjm1                                 !  Vertical slab 
    656          DO ji = 2, jpim1 
    657             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    658             zbeta    = rab_n(ji,jj,1,jp_sal) 
    659             IF ( lconv(ji,jj) ) THEN 
    660                zt   = 0._wp 
    661                zs   = 0._wp 
    662                zu   = 0._wp 
    663                zv   = 0._wp 
    664                ! average over depth of boundary layer 
    665                zthick=0._wp 
    666                DO jm = 2, imld(ji,jj) 
    667                   zthick=zthick+e3t_n(ji,jj,jm) 
    668                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    669                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    670                   zu   = zu  + e3t_n(ji,jj,jm) & 
    671                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    672                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    673                   zv   = zv  + e3t_n(ji,jj,jm) & 
    674                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    675                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    676                END DO 
    677                zt_ml(ji,jj) = zt / zthick 
    678                zs_ml(ji,jj) = zs / zthick 
    679                zu_ml(ji,jj) = zu / zthick 
    680                zv_ml(ji,jj) = zv / zthick 
    681                zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    682                zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    683                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    684                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    685                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    686                      &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    687                zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    688             ELSE 
    689             ! stable, if entraining calulate average below interface layer. 
    690                IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    691                   zt   = 0._wp 
    692                   zs   = 0._wp 
    693                   zu   = 0._wp 
    694                   zv   = 0._wp 
    695                   ! average over depth of boundary layer 
    696                   zthick=0._wp 
    697                   DO jm = 2, imld(ji,jj) 
    698                      zthick=zthick+e3t_n(ji,jj,jm) 
    699                      zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    700                      zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    701                      zu   = zu  + e3t_n(ji,jj,jm) & 
    702                         &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    703                         &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    704                      zv   = zv  + e3t_n(ji,jj,jm) & 
    705                         &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    706                         &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    707                   END DO 
    708                   zt_ml(ji,jj) = zt / zthick 
    709                   zs_ml(ji,jj) = zs / zthick 
    710                   zu_ml(ji,jj) = zu / zthick 
    711                   zv_ml(ji,jj) = zv / zthick 
    712                   zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    713                   zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    714                   zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    715                         &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    716                   zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    717                         &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    718                   zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    719                ENDIF 
    720             ENDIF 
    721          END DO 
    722       END DO 
    723     ! 
     551! 
     552    ! Average over the depth of the mixed layer in the convective boundary layer 
     553    ! Alan: do we need zb_ml? 
     554    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    724555    ! rotate mean currents and changes onto wind align co-ordinates 
    725556    ! 
    726  
    727       DO jj = 2, jpjm1 
    728          DO ji = 2, jpim1 
    729             ztemp = zu_ml(ji,jj) 
    730             zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
    731             zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    732             ztemp = zdu_ml(ji,jj) 
    733             zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
    734             zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    735     ! 
    736             ztemp = zu_bl(ji,jj) 
    737             zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
    738             zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    739             ztemp = zdu_bl(ji,jj) 
    740             zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
    741             zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    742          END DO 
    743       END DO 
    744  
     557    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     558    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    745559     zuw_bse = 0._wp 
    746560     zvw_bse = 0._wp 
     561     zwth_ent = 0._wp 
     562     zws_ent = 0._wp 
    747563     DO jj = 2, jpjm1 
    748564        DO ji = 2, jpim1 
    749  
    750            IF ( lconv(ji,jj) ) THEN 
    751               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    752                  zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    753                  zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     565           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 
     566              IF ( lconv(ji,jj) ) THEN 
     567             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 
     568                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 
     569                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 
     570            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 
     571                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 
     572                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 
     573                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     574                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     575                 ENDIF 
     576              ELSE 
     577                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
     578                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    754579              ENDIF 
    755            ELSE 
    756               zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    757               zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    758580           ENDIF 
    759581        END DO 
     
    764586      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    765587 
    766        DO jj = 2, jpjm1 
    767           DO ji = 2, jpim1 
    768           ! 
    769              IF ( lconv (ji,jj) ) THEN 
    770              ! Unstable conditions 
    771                 IF( zdb_bl(ji,jj) > 0._wp ) THEN 
    772                 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
    773                    ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
    774                    zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
    775                    zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    776                    DO jk = 2 , ibld(ji,jj) 
    777                       znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    778                       zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    779                       zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    780                       zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    781                    END DO 
    782                 ENDIF 
    783              ELSE 
    784              ! stable conditions 
    785              ! if pycnocline profile only defined when depth steady of increasing. 
    786                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
    787                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    788                      IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    789                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
    790                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
    791                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    792                          DO jk = 2, ibld(ji,jj) 
    793                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    794                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    795                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    796                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    797                          END DO 
    798                      ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    799                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
    800                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
    801                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    802                          DO jk = 2, ibld(ji,jj) 
    803                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    804                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    805                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    806                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    807                          END DO 
    808                       ENDIF ! IF (zhol >=0.5) 
    809                    ENDIF    ! IF (zdb_bl> 0.) 
    810                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
    811              ENDIF          ! IF (lconv) 
    812             ! 
    813           END DO 
    814        END DO 
    815 ! 
    816        DO jj = 2, jpjm1 
    817           DO ji = 2, jpim1 
    818           ! 
    819              IF ( lconv (ji,jj) ) THEN 
    820              ! Unstable conditions 
    821                  zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
    822                & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
    823                 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
    824               & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    825                 DO jk = 2 , ibld(ji,jj)-1 
    826                    znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    827                    zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    828                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    829                 END DO 
    830              ELSE 
    831              ! stable conditions 
    832                 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    833                 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    834                 DO jk = 2, ibld(ji,jj) 
    835                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    836                    IF ( znd < 1.0 ) THEN 
    837                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    838                    ELSE 
    839                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    840                    ENDIF 
    841                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
    842                 END DO 
    843              ENDIF 
    844             ! 
    845           END DO 
    846        END DO 
     588      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     589      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 
     590      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    847591       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    848592       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    849593       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    850594 
    851       ! WHERE ( lconv ) 
    852       !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 
    853       !     zvisml_sc = zdifml_sc 
    854       !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 
    855       !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 
    856       !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 
    857       !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 
    858       !  ELSEWHERE 
    859       !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    860       !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    861       !  ENDWHERE 
    862595       DO jj = 2, jpjm1 
    863596          DO ji = 2, jpim1 
     
    868601               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    869602               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    870                zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
     603               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 
    871604             ELSE 
    872605               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    873606               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    874             END IF 
    875         END DO 
    876     END DO 
     607             END IF 
     608          END DO 
     609       END DO 
    877610! 
    878611       DO jj = 2, jpjm1 
     
    889622                ! pycnocline - if present linear profile 
    890623                IF ( zdh(ji,jj) > 0._wp ) THEN 
     624                   zgamma_b = 6.0 
    891625                   DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892626                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    893627                       ! 
    894                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     628                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    895629                       ! 
    896                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     630                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    897631                   END DO 
     632                   IF ( ibld_ext == 0 ) THEN 
     633                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     634                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     635                   ELSE 
     636                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     637                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     638                   ENDIF 
    898639                ENDIF 
    899                 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    900                 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 
     640                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
     641                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 
     642                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 
    901643                ! could be taken out, take account of entrainment represents as a diffusivity 
    902644                ! should remove w from here, represents entrainment 
     
    908650                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    909651                END DO 
     652 
     653                IF ( ibld_ext == 0 ) THEN 
     654                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     655                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     656                ELSE 
     657                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     658                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     659                ENDIF 
    910660             ENDIF   ! end if ( lconv ) 
    911 ! 
     661             ! 
    912662          END DO  ! end of ji loop 
    913663       END DO     ! end of jj loop 
     
    952702       END DO     ! end of jj loop 
    953703 
    954  
    955704! 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) 
    956705       WHERE ( lconv ) 
    957           zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ) 
    958           zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln ) 
    959           zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
     706          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 ) 
     707          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
     708          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 ) 
    960709       ELSEWHERE 
    961710          zsc_uw_1 = zustar**2 
    962           zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 
     711          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    963712       ENDWHERE 
    964  
     713       IF(ln_dia_osm) THEN 
     714          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     715          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     716       END IF 
    965717       DO jj = 2, jpjm1 
    966718          DO ji = 2, jpim1 
     
    1007759                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1008760                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    1009                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     761                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
    1010762                   ! non-gradient buoyancy terms 
    1011763                   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 ) 
     
    1020772          END DO   ! ji loop 
    1021773       END DO      ! jj loop 
    1022  
    1023774 
    1024775       WHERE ( lconv ) 
     
    1051802       END DO           ! jj loop 
    1052803 
     804       IF(ln_dia_osm) THEN 
     805          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
     806          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     807       END IF 
    1053808! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    1054809 
     
    1089844       END DO      ! jj loop 
    1090845 
    1091  
    1092846       WHERE ( lconv ) 
    1093847          zsc_uw_1 = zustar**2 
     
    1134888          END DO 
    1135889       END DO 
     890 
     891       IF(ln_dia_osm) THEN 
     892          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
     893          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
     894          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     895          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     896          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     897          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     898       END IF 
    1136899! 
    1137900! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
     
    1165928      END DO 
    1166929 
     930       IF(ln_dia_osm) THEN 
     931          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     932          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     933       END IF 
    1167934      ! pynocline contributions 
    1168935       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
     
    1170937       DO jj = 2, jpjm1 
    1171938          DO ji = 2, jpim1 
    1172              DO jk= 2, ibld(ji,jj) 
    1173                 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1174                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1175                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1176                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1177                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
    1178                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    1179              END DO 
    1180            END DO 
     939             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     940                DO jk= 2, ibld(ji,jj) 
     941                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     942                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     943                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     944                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     945                   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) 
     946                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     947                END DO 
     948             END IF 
     949          END DO 
    1181950       END DO 
    1182951 
     
    1185954       DO jj=2, jpjm1 
    1186955          DO ji = 2, jpim1 
    1187              IF ( lconv(ji,jj) ) THEN 
     956            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 
    1188957               DO jk = 1, imld(ji,jj) - 1 
    1189958                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1190                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1191                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
     959                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
     960                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1192961                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1193962                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
     
    1195964               DO jk = imld(ji,jj), ibld(ji,jj) 
    1196965                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    1197                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1198                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
     966                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
     967                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1199968                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1200969                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1201970                END DO 
    1202971             ENDIF 
    1203              ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
    1204              ghams(ji,jj,ibld(ji,jj)) = 0._wp 
    1205              ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
    1206              ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
     972 
     973             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     974             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     975             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     976             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1207977          END DO       ! ji loop 
    1208978       END DO          ! jj loop 
    1209979 
    1210  
     980       IF(ln_dia_osm) THEN 
     981          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     982          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     983          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 
     984          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 
     985          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
     986          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
     987          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
     988       END IF 
    1211989       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1212990       ! Need to put in code for contributions that are applied explicitly to 
     
    12291007          END DO 
    12301008       END DO 
    1231  
    1232        IF(ln_dia_osm) THEN 
    1233           IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
    1234        END IF 
    12351009 
    12361010! KPP-style Ri# mixing 
     
    12861060       END IF ! ln_convmix = .true. 
    12871061 
     1062 
     1063  
     1064       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
     1065          DO jj = 2 , jpjm1 
     1066             DO ji = 2, jpim1 
     1067                IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 
     1068            ! Calculate MLE flux profiles 
     1069                   ! DO jk = 1, mld_prof(ji,jj) 
     1070                   !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1071                   !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 
     1072                   !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1073                   !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 
     1074                   !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1075                   ! END DO 
     1076            ! Calculate MLE flux contribution from surface fluxes 
     1077                   DO jk = 1, ibld(ji,jj) 
     1078                     znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 
     1079                     ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 
     1080                     ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 
     1081                    END DO 
     1082                    DO jk = 1, mld_prof(ji,jj) 
     1083                      znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1084                      ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 
     1085                      ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 
     1086                    END DO 
     1087            ! Viscosity for MLEs 
     1088                    DO jk = ibld(ji,jj), mld_prof(ji,jj) 
     1089                      zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 
     1090                    END DO 
     1091            ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 
     1092                    jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 
     1093                    jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 
     1094                    DO jk = mld_prof(ji,jj),  jl 
     1095                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 
     1096                            & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 
     1097                            & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 
     1098                    END DO 
     1099                ENDIF 
     1100            END DO 
     1101          END DO 
     1102       ENDIF 
     1103 
     1104       IF(ln_dia_osm) THEN 
     1105          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1106          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1107          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
     1108       END IF 
     1109 
     1110 
    12881111       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
    1289        CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. ) 
     1112       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 
    12901113 
    12911114       ! GN 25/8: need to change tmask --> wmask 
     
    13191142        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    13201143        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   & 
    1321          &                  ghamu, 'U', 1. , ghamv, 'V', 1. ) 
    1322  
    1323        IF(ln_dia_osm) THEN 
     1144         &                  ghamu, 'U', -1. , ghamv, 'V', -1. ) 
     1145 
     1146      IF(ln_dia_osm) THEN 
    13241147         SELECT CASE (nn_osm_wave) 
    13251148         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1). 
     
    13301153         ! Stokes drift read in from sbcwave  (=2). 
    13311154         CASE(2) 
    1332             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    1333             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
     1155            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
     1156            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     1157            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period 
     1158            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height 
     1159            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 
     1160            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum 
     1161            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10 
    13341162            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
    13351163                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
     
    13421170         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0> 
    13431171         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth 
    1344          IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth 
     1172         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1173         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
     1174         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
     1175         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base 
     1176         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base 
     1177         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base 
     1178         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth 
     1179         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth 
    13451180         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth 
    13461181         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points 
     
    13481183         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    13491184         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
     1185         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale 
     1186         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir # 
    13501187         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    13511188         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    13521189         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    13531190         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    1354          IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine 
     1191         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1192         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    13551193         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    1356          IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine 
    1357          IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine 
    1358          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML 
     1194         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
     1195         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
     1196         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
     1197         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
     1198         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
     1199 
     1200         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )               ! FK layer depth 
     1201         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )               ! FK target layer depth 
     1202         IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux 
     1203         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML 
     1204         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 
     1205         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )            ! FK dtdx at u-pt 
     1206         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )            ! FK dtdy at v-pt 
     1207         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )            ! FK dtdx at u-pt 
     1208         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )            ! FK dsdy at v-pt 
     1209         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )            ! FK dbdx at u-pt 
     1210         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )            ! FK dbdy at v-pt 
     1211         IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1212         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1213 
    13591214      END IF 
    1360       ! Lateral boundary conditions on p_avt  (sign unchanged) 
    1361       CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 
     1215 
     1216CONTAINS 
     1217 
     1218 
     1219   ! Alan: do we need zb? 
     1220   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1221     !!--------------------------------------------------------------------- 
     1222     !!                   ***  ROUTINE zdf_vertical_average  *** 
     1223     !! 
     1224     !! ** Purpose : Determines vertical averages from surface to jnlev. 
     1225     !! 
     1226     !! ** Method  : Averages are calculated from the surface to jnlev. 
     1227     !!              The external level used to calculate differences is ibld+ibld_ext 
     1228     !! 
     1229     !!---------------------------------------------------------------------- 
     1230 
     1231        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1232 
     1233        ! Alan: do we need zb? 
     1234        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity 
     1235        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
     1236        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
     1237        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
     1238 
     1239        INTEGER :: jk, ji, jj 
     1240        REAL(wp) :: zthick, zthermal, zbeta 
     1241 
     1242 
     1243        zt   = 0._wp 
     1244        zs   = 0._wp 
     1245        zu   = 0._wp 
     1246        zv   = 0._wp 
     1247        DO jj = 2, jpjm1                                 !  Vertical slab 
     1248         DO ji = 2, jpim1 
     1249            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1250            zbeta    = rab_n(ji,jj,1,jp_sal) 
     1251               ! average over depth of boundary layer 
     1252            zthick = epsln 
     1253            DO jk = 2, jnlev_av(ji,jj) 
     1254               zthick = zthick + e3t_n(ji,jj,jk) 
     1255               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 
     1256               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     1257               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) & 
     1258                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 
     1259                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 
     1260               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) & 
     1261                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 
     1262                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 
     1263            END DO 
     1264            zt(ji,jj) = zt(ji,jj) / zthick 
     1265            zs(ji,jj) = zs(ji,jj) / zthick 
     1266            zu(ji,jj) = zu(ji,jj) / zthick 
     1267            zv(ji,jj) = zv(ji,jj) / zthick 
     1268            ! Alan: do we need zb? 
     1269            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 
     1270            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 
     1271            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 
     1272                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 
     1273            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 
     1274                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 
     1275            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1276         END DO 
     1277        END DO 
     1278   END SUBROUTINE zdf_osm_vertical_average 
     1279 
     1280   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 
     1281     !!--------------------------------------------------------------------- 
     1282     !!                   ***  ROUTINE zdf_velocity_rotation  *** 
     1283     !! 
     1284     !! ** Purpose : Rotates frame of reference of averaged velocity components. 
     1285     !! 
     1286     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w 
     1287     !! 
     1288     !!---------------------------------------------------------------------- 
     1289 
     1290        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle 
     1291        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current 
     1292        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline 
     1293 
     1294        INTEGER :: ji, jj 
     1295        REAL(wp) :: ztemp 
     1296 
     1297        DO jj = 2, jpjm1 
     1298           DO ji = 2, jpim1 
     1299              ztemp = zu(ji,jj) 
     1300              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 
     1301              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1302              ztemp = zdu(ji,jj) 
     1303              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
     1304              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1305           END DO 
     1306        END DO 
     1307    END SUBROUTINE zdf_osm_velocity_rotation 
     1308 
     1309    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 
     1310     !!--------------------------------------------------------------------- 
     1311     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     1312     !! 
     1313     !! ** Purpose : Calculates the gradients below the OSBL 
     1314     !! 
     1315     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1316     !! 
     1317     !!---------------------------------------------------------------------- 
     1318 
     1319     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
     1320 
     1321     INTEGER :: jj, ji, jkb, jkb1 
     1322     REAL(wp) :: zthermal, zbeta 
     1323 
     1324 
     1325     DO jj = 2, jpjm1 
     1326        DO ji = 2, jpim1 
     1327           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1328              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1329              zbeta    = rab_n(ji,jj,1,jp_sal) 
     1330              jkb = ibld(ji,jj) + ibld_ext 
     1331              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     1332              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 
     1333                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1334              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 
     1335                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1336              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1337           END IF 
     1338        END DO 
     1339     END DO 
     1340    END SUBROUTINE zdf_osm_external_gradients 
     1341 
     1342    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 
     1343 
     1344     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1345 
     1346     INTEGER :: jk, jj, ji 
     1347     REAL(wp) :: ztgrad, zsgrad, zbgrad 
     1348     REAL(wp) :: zgamma_b_nd, zgamma_c, znd 
     1349     REAL(wp) :: zzeta_s=0.3 
     1350 
     1351     DO jj = 2, jpjm1 
     1352        DO ji = 2, jpim1 
     1353           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1354              IF ( lconv(ji,jj) ) THEN  ! convective conditions 
     1355                    IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 
     1356                         & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 
     1357                         &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0 
     1358                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj) 
     1359                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj) 
     1360                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj) 
     1361                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1362                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 
     1363                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 
     1364                       DO jk = 2, ibld(ji,jj)+ibld_ext 
     1365                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     1366                          IF ( znd <= zzeta_s ) THEN 
     1367                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * & 
     1368                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1369                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * & 
     1370                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1371                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * & 
     1372                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1373                          ELSE 
     1374                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1375                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1376                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1377                          ENDIF 
     1378                       END DO 
     1379                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 
     1380              ELSE 
     1381                 ! stable conditions 
     1382                 ! if pycnocline profile only defined when depth steady of increasing. 
     1383                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     1384                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1385                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     1386                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
     1387                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
     1388                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
     1389                          DO jk = 2, ibld(ji,jj) 
     1390                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1391                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1392                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1393                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1394                          END DO 
     1395                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1396                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
     1397                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
     1398                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
     1399                          DO jk = 2, ibld(ji,jj) 
     1400                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     1401                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1402                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1403                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1404                          END DO 
     1405                       ENDIF ! IF (zhol >=0.5) 
     1406                    ENDIF    ! IF (zdb_bl> 0.) 
     1407                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1408              ENDIF          ! IF (lconv) 
     1409           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1410        END DO 
     1411     END DO 
     1412 
     1413    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
     1414 
     1415    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
     1416      !!--------------------------------------------------------------------- 
     1417      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
     1418      !! 
     1419      !! ** Purpose : Calculates velocity shear in the pycnocline 
     1420      !! 
     1421      !! ** Method  : 
     1422      !! 
     1423      !!---------------------------------------------------------------------- 
     1424 
     1425      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
     1426 
     1427      INTEGER :: jk, jj, ji 
     1428      REAL(wp) :: zugrad, zvgrad, znd 
     1429      REAL(wp) :: zzeta_v = 0.45 
    13621430      ! 
    1363    END SUBROUTINE zdf_osm 
     1431      DO jj = 2, jpjm1 
     1432         DO ji = 2, jpim1 
     1433            ! 
     1434            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1435               IF ( lconv (ji,jj) ) THEN 
     1436                  ! Unstable conditions 
     1437                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1438                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1439                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1440                  !Alan is this right? 
     1441                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1442                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1443                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1444                       &      )/ (zdh(ji,jj)  + epsln ) 
     1445                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     1446                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     1447                     IF ( znd <= 0.0 ) THEN 
     1448                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     1449                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     1450                     ELSE 
     1451                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     1452                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     1453                     ENDIF 
     1454                  END DO 
     1455               ELSE 
     1456                  ! stable conditions 
     1457                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     1458                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     1459                  DO jk = 2, ibld(ji,jj) 
     1460                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1461                     IF ( znd < 1.0 ) THEN 
     1462                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     1463                     ELSE 
     1464                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     1465                     ENDIF 
     1466                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     1467                  END DO 
     1468               ENDIF 
     1469               ! 
     1470            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1471         END DO 
     1472      END DO 
     1473    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1474 
     1475    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
     1476     !!--------------------------------------------------------------------- 
     1477     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     1478     !! 
     1479     !! ** Purpose : Calculates the rate at which hbl changes. 
     1480     !! 
     1481     !! ** Method  : 
     1482     !! 
     1483     !!---------------------------------------------------------------------- 
     1484 
     1485    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl 
     1486 
     1487    INTEGER :: jj, ji 
     1488    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 
     1489    REAL(wp) :: zvel_max, zwb_min 
     1490    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1491    REAL(wp) :: zzeta_m = 0.3 
     1492    REAL(wp) :: zgamma_c = 2.0 
     1493    REAL(wp) :: zdhoh = 0.1 
     1494    REAL(wp) :: alpha_bc = 0.5 
     1495 
     1496    DO jj = 2, jpjm1 
     1497       DO ji = 2, jpim1 
     1498          IF ( lconv(ji,jj) ) THEN    ! Convective 
     1499             ! Alan is this right?  Yes, it's a bit different from the previous relationship 
     1500             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1501             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1502             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1503             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1504             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1505             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
     1506             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1507                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1508 
     1509             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1510                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1511                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1512                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1513             ! 
     1514             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 
     1515 
     1516             IF ( ln_osm_mle ) THEN 
     1517                  !  zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 
     1518                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL 
     1519                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1520                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     1521                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     1522                ELSE 
     1523                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     1524                ENDIF 
     1525                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1526                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     1527                   ! OSBL is deepening, entrainment > restratification 
     1528                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1529                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1530                   ELSE 
     1531                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     1532                   ENDIF 
     1533                ELSE 
     1534                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     1535                   zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     1536 
     1537 
     1538                ENDIF 
     1539 
     1540             ELSE 
     1541                ! Fox-Kemper not used. 
     1542 
     1543                  zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     1544                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1545                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1546                ! added ajgn 23 July as temporay fix 
     1547 
     1548             ENDIF 
     1549 
     1550             zdhdt_2(ji,jj) = 0._wp 
     1551 
     1552                ! commented out ajgn 23 July as temporay fix 
     1553!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1554! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 
     1555!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1556!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 
     1557!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1558!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 
     1559!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 
     1560!                     ! Alan no idea what this should be? 
     1561!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 
     1562!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 
     1563!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 
     1564!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 
     1565!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 
     1566!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 
     1567!                 ENDIF 
     1568            ELSE                        ! Stable 
     1569                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     1570                zdhdt_2(ji,jj) = 0._wp 
     1571                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     1572                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     1573                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     1574                ELSE 
     1575                    zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     1576                ENDIF 
     1577                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 
     1578            ENDIF 
     1579         END DO 
     1580      END DO 
     1581    END SUBROUTINE zdf_osm_calculate_dhdt 
     1582 
     1583    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     1584     !!--------------------------------------------------------------------- 
     1585     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
     1586     !! 
     1587     !! ** Purpose : Increments hbl. 
     1588     !! 
     1589     !! ** Method  : If thechange in hbl exceeds one model level the change is 
     1590     !!              is calculated by moving down the grid, changing the buoyancy 
     1591     !!              jump. This is to ensure that the change in hbl does not 
     1592     !!              overshoot a stable layer. 
     1593     !! 
     1594     !!---------------------------------------------------------------------- 
     1595 
     1596 
     1597    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl. 
     1598 
     1599    INTEGER :: jk, jj, ji, jm 
     1600    REAL(wp) :: zhbl_s, zvel_max, zdb 
     1601    REAL(wp) :: zthermal, zbeta 
     1602 
     1603     DO jj = 2, jpjm1 
     1604         DO ji = 2, jpim1 
     1605           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     1606! 
     1607! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     1608! 
     1609              zhbl_s = hbl(ji,jj) 
     1610              jm = imld(ji,jj) 
     1611              zthermal = rab_n(ji,jj,1,jp_tem) 
     1612              zbeta = rab_n(ji,jj,1,jp_sal) 
     1613 
     1614 
     1615              IF ( lconv(ji,jj) ) THEN 
     1616!unstable 
     1617 
     1618                 IF( ln_osm_mle ) THEN 
     1619                    zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1620                 ELSE 
     1621 
     1622                    zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     1623                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1624 
     1625                 ENDIF 
     1626 
     1627                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1628                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
     1629                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 
     1630                         &  0.0 ) + zvel_max 
     1631 
     1632 
     1633                    IF ( ln_osm_mle ) THEN 
     1634                       zhbl_s = zhbl_s + MIN( & 
     1635                        & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1636                        & e3w_n(ji,jj,jm) ) 
     1637                    ELSE 
     1638                      zhbl_s = zhbl_s + MIN( & 
     1639                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1640                        & e3w_n(ji,jj,jm) ) 
     1641                    ENDIF 
     1642 
     1643                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1644 
     1645                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
     1646                 END DO 
     1647                 hbl(ji,jj) = zhbl_s 
     1648                 ibld(ji,jj) = jm 
     1649             ELSE 
     1650! stable 
     1651                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1652                    zdb = MAX( & 
     1653                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 
     1654                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 
     1655                         & 0.0 ) + & 
     1656             &       2.0 * zvstr(ji,jj)**2 / zhbl_s 
     1657 
     1658                    ! Alan is thuis right? I have simply changed hbli to hbl 
     1659                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 
     1660                    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) ) ) * & 
     1661               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 
     1662                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 
     1663                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 
     1664 
     1665                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1666                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
     1667                 END DO 
     1668             ENDIF   ! IF ( lconv ) 
     1669             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 
     1670             ibld(ji,jj) = MAX(jm, 4 ) 
     1671           ELSE 
     1672! change zero or one model level. 
     1673             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 
     1674           ENDIF 
     1675           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     1676         END DO 
     1677      END DO 
     1678 
     1679    END SUBROUTINE zdf_osm_timestep_hbl 
     1680 
     1681    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     1682      !!--------------------------------------------------------------------- 
     1683      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     1684      !! 
     1685      !! ** Purpose : Calculates thickness of the pycnocline 
     1686      !! 
     1687      !! ** Method  : The thickness is calculated from a prognostic equation 
     1688      !!              that relaxes the pycnocine thickness to a diagnostic 
     1689      !!              value. The time change is calculated assuming the 
     1690      !!              thickness relaxes exponentially. This is done to deal 
     1691      !!              with large timesteps. 
     1692      !! 
     1693      !!---------------------------------------------------------------------- 
     1694 
     1695      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
     1696      ! 
     1697      INTEGER :: jj, ji 
     1698      INTEGER :: inhml 
     1699      REAL(wp) :: zari, ztau, zddhdt 
     1700 
     1701 
     1702      DO jj = 2, jpjm1 
     1703         DO ji = 2, jpim1 
     1704            IF ( lconv(ji,jj) ) THEN 
     1705 
     1706               IF( ln_osm_mle ) THEN 
     1707                  IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
     1708                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
     1709                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     1710                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1711                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1712                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1713                        ELSE                                                     ! unstable 
     1714                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1715                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1716                        ENDIF 
     1717                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1718                        zddhdt = zari * hbl(ji,jj) 
     1719                     ELSE 
     1720                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1721                        zddhdt = 0.2 * hbl(ji,jj) 
     1722                     ENDIF 
     1723                  ELSE 
     1724                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1725                     zddhdt = 0.2 * hbl(ji,jj) 
     1726                  ENDIF 
     1727               ELSE ! ln_osm_mle 
     1728                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     1729                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1730                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1731                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1732                     ELSE                                                     ! unstable 
     1733                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1734                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1735                     ENDIF 
     1736                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1737                     zddhdt = zari * hbl(ji,jj) 
     1738                  ELSE 
     1739                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1740                     zddhdt = 0.2 * hbl(ji,jj) 
     1741                  ENDIF 
     1742 
     1743               END IF  ! ln_osm_mle 
     1744 
     1745               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1746               ! Alan: this hml is never defined or used 
     1747            ELSE   ! IF (lconv) 
     1748               ztau = hbl(ji,jj) / zvstr(ji,jj) 
     1749               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     1750                  ! boundary layer deepening 
     1751                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1752                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     1753                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     1754                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     1755                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
     1756                  ELSE 
     1757                     zddhdt = 0.2 * hbl(ji,jj) 
     1758                  ENDIF 
     1759               ELSE     ! IF(dhdt < 0) 
     1760                  zddhdt = 0.2 * hbl(ji,jj) 
     1761               ENDIF    ! IF (dhdt >= 0) 
     1762               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1763               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse 
     1764               ! Alan: this hml is never defined or used -- do we need it? 
     1765            ENDIF       ! IF (lconv) 
     1766 
     1767            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     1768            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     1769            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     1770            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     1771            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1772         END DO 
     1773      END DO 
     1774 
     1775    END SUBROUTINE zdf_osm_pycnocline_thickness 
     1776 
     1777 
     1778   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     1779      !!---------------------------------------------------------------------- 
     1780      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  *** 
     1781      !! 
     1782      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 
     1783      !! 
     1784      !! ** Method  : 
     1785      !! 
     1786      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1787      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1788 
     1789 
     1790      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 
     1791      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
     1792      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy 
     1793 
     1794      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1795      INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     1796      REAL(wp)                         :: zc 
     1797      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value 
     1798      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH 
     1799      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 
     1800      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv 
     1801!!---------------------------------------------------------------------- 
     1802      ! 
     1803      !                                      !==  MLD used for MLE  ==! 
     1804 
     1805      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point 
     1806      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
     1807      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     1808      DO jk = nlb10, jpkm1 
     1809         DO jj = 1, jpj                ! Mixed layer level: w-level  
     1810            DO ji = 1, jpi 
     1811               ikt = mbkt(ji,jj) 
     1812               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 
     1813               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     1814            END DO 
     1815         END DO 
     1816      END DO 
     1817      DO jj = 1, jpj 
     1818         DO ji = 1, jpi 
     1819            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 
     1820            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1821         END DO 
     1822      END DO 
     1823      ! ensure mld_prof .ge. ibld 
     1824      ! 
     1825      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     1826      ! 
     1827      ztm(:,:) = 0._wp 
     1828      zsm(:,:) = 0._wp 
     1829      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
     1830         DO jj = 1, jpj 
     1831            DO ji = 1, jpi 
     1832               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     1833               ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 
     1834               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 
     1835            END DO 
     1836         END DO 
     1837      END DO 
     1838      ! average temperature and salinity. 
     1839      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1840      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1841      ! calculate horizontal gradients at u & v points 
     1842 
     1843      DO jj = 2, jpjm1 
     1844         DO ji = 1, jpim1 
     1845            zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1846            zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1847            zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 
     1848            ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 
     1849            ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 
     1850         END DO 
     1851      END DO 
     1852 
     1853      DO jj = 1, jpjm1 
     1854         DO ji = 2, jpim1 
     1855            zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1856            zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1857            zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 
     1858            ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 
     1859            ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 
     1860         END DO 
     1861      END DO 
     1862 
     1863      CALL eos_rab(ztsm_midu, zmld_midu, zabu) 
     1864      CALL eos_rab(ztsm_midv, zmld_midv, zabv) 
     1865 
     1866      DO jj = 2, jpjm1 
     1867         DO ji = 1, jpim1 
     1868            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 
     1869         END DO 
     1870      END DO 
     1871      DO jj = 1, jpjm1 
     1872         DO ji = 2, jpim1 
     1873            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 
     1874         END DO 
     1875      END DO 
     1876 
     1877 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
     1878  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     1879      !!---------------------------------------------------------------------- 
     1880      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
     1881      !! 
     1882      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 
     1883      !! 
     1884      !! ** Method  : 
     1885      !! 
     1886      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1887      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1888 
     1889      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle 
     1890      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1891      INTEGER  ::   ii, ij, ik  ! local integers 
     1892      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     1893      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle 
     1894 
     1895      mld_prof(:,:) = 4 
     1896      DO jk = 5, jpkm1 
     1897        DO jj = 2, jpjm1 
     1898          DO ji = 2, jpim1 
     1899            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1900          END DO 
     1901        END DO 
     1902      END DO 
     1903      ! DO jj = 2, jpjm1 
     1904      !    DO ji = 1, jpim1 
     1905      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1906      !    END DO 
     1907      !  END DO 
     1908   ! Timestep mixed layer eddy depth. 
     1909      DO jj = 2, jpjm1 
     1910        DO ji = 2, jpim1 
     1911          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 
     1912          IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 
     1913             hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh )  ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 
     1914          ELSE 
     1915            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 
     1916             ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 
     1917             !   hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     1918             ! ELSE 
     1919             !   hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
     1920             ! ENDIF 
     1921             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1922               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     1923             ELSE 
     1924               hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
     1925             ENDIF 
     1926          ENDIF 
     1927          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
     1928        END DO 
     1929      END DO 
     1930 
     1931      mld_prof = 4 
     1932      DO jk = 5, jpkm1 
     1933        DO jj = 2, jpjm1 
     1934          DO ji = 2, jpim1 
     1935            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1936          END DO 
     1937        END DO 
     1938      END DO 
     1939      DO jj = 2, jpjm1 
     1940         DO ji = 2, jpim1 
     1941            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 
     1942         END DO 
     1943       END DO 
     1944   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
     1945 
     1946      DO jj = 2, jpjm1 
     1947        DO ji = 2, jpim1 
     1948          IF ( lconv(ji,jj) ) THEN 
     1949             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     1950             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 
     1951             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 
     1952             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 
     1953             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 
     1954             zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
     1955                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
     1956             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 
     1957      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
     1958             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1)  
     1959             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf 
     1960          ENDIF 
     1961        END DO 
     1962      END DO 
     1963END SUBROUTINE zdf_osm_mle_parameters 
     1964 
     1965END SUBROUTINE zdf_osm 
    13641966 
    13651967 
     
    13781980     INTEGER  ::   ios            ! local integer 
    13791981     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     1982     REAL z1_t2 
    13801983     !! 
    13811984     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    13821985          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
    1383           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 
     1986          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
     1987          & ,ln_osm_mle 
     1988! Namelist for Fox-Kemper parametrization. 
     1989      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
     1990           & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 
     1991 
    13841992     !!---------------------------------------------------------------------- 
    13851993     ! 
     
    13972005        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 
    13982006        WRITE(numout,*) '~~~~~~~~~~~~' 
    1399         WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters' 
    1400         WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la 
     2007        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters' 
     2008        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la 
     2009        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    14012010        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
    14022011        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    1403         WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes 
     2012        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
    14042013        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave 
    14052014        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave 
     
    14232032     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    14242033 
    1425      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     2034 
     2035     IF( ln_osm_mle ) THEN 
     2036! Initialise Fox-Kemper parametrization 
     2037         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme 
     2038         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 
     2039903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 
     2040 
     2041         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 
     2042         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 
     2043904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 
     2044         IF(lwm) WRITE ( numond, namosm_mle ) 
     2045 
     2046         IF(lwp) THEN                     ! Namelist print 
     2047            WRITE(numout,*) 
     2048            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 
     2049            WRITE(numout,*) '~~~~~~~~~~~~~' 
     2050            WRITE(numout,*) '   Namelist namosm_mle : ' 
     2051            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle 
     2052            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce 
     2053            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm' 
     2054            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's' 
     2055            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg' 
     2056            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c 
     2057            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
     2058            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
     2059         ENDIF         ! 
     2060     ENDIF 
     2061      ! 
     2062      IF(lwp) THEN 
     2063         WRITE(numout,*) 
     2064         IF( ln_osm_mle ) THEN 
     2065            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 
     2066            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation' 
     2067            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation' 
     2068         ELSE 
     2069            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 
     2070         ENDIF 
     2071      ENDIF 
     2072      ! 
     2073      IF( ln_osm_mle ) THEN                ! MLE initialisation 
     2074         ! 
     2075         rb_c = grav * rn_osm_mle_rho_c /rau0        ! Mixed Layer buoyancy criteria 
     2076         IF(lwp) WRITE(numout,*) 
     2077         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
     2078         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 
     2079         ! 
     2080         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            ! 
     2081! 
     2082         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation 
     2083            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  ) 
     2084            ! 
     2085         ENDIF 
     2086         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
     2087         z1_t2 = 2.e-5 
     2088         do jj=1,jpj 
     2089            do ji = 1,jpi 
     2090               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 
     2091            end do 
     2092         end do 
     2093         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 
     2094         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  ) 
     2095         ! 
     2096      ENDIF 
     2097 
     2098     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 
     2099 
    14262100 
    14272101     IF( ln_zdfddm) THEN 
     
    15362210     REAL(wp) ::   zN2_c           ! local scalar 
    15372211     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    1538      INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top) 
     2212     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 
    15392213     !!---------------------------------------------------------------------- 
    15402214     ! 
     
    15512225           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    15522226        END IF 
     2227 
    15532228        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    1554         id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. ) 
     2229        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. ) 
    15552230        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return 
    15562231           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 
    1557            CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios  ) 
    1558            WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file' 
     2232           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  ) 
     2233           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file' 
    15592234           RETURN 
    1560         ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate 
     2235        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate 
    15612236           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 
    15622237        END IF 
     
    15702245         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
    15712246         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    1572          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
     2247         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios ) 
    15732248        RETURN 
    15742249     END IF 
     
    15782253     !!----------------------------------------------------------------------------- 
    15792254     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 
    1580      ALLOCATE( imld_rst(jpi,jpj) ) 
    15812255     ! w-level of the mixing and mixed layers 
    15822256     CALL eos_rab( tsn, rab_n ) 
     
    15992273     DO jj = 1, jpj 
    16002274        DO ji = 1, jpi 
    1601            iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
     2275           iiki = MAX(4,imld_rst(ji,jj)) 
     2276           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth 
     2277           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth 
    16032278        END DO 
    16042279     END DO 
    1605      hbl = MAX(hbl,epsln) 
    1606      hbli(:,:) = hbl(:,:) 
    1607      DEALLOCATE( imld_rst ) 
     2280 
     2281     IF( ln_osm_mle ) hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2282 
    16082283     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
     2284     wn(:,:,:) = 0._wp 
     2285     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    16092286   END SUBROUTINE osm_rst 
    16102287 
     
    16342311      ENDIF 
    16352312 
    1636       ! add non-local temperature and salinity flux 
    16372313      DO jk = 1, jpkm1 
    16382314         DO jj = 2, jpjm1 
     
    16482324      END DO 
    16492325 
    1650  
    1651       ! save the non-local tracer flux trends for diagnostic 
     2326      ! save the non-local tracer flux trends for diagnostics 
    16522327      IF( l_trdtra )   THEN 
    16532328         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    16542329         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    1655 !!bug gm jpttdzdf ==> jpttosm 
    1656          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    1657          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     2330 
     2331         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 
     2332         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 
    16582333         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    16592334      ENDIF 
     
    17232398 
    17242399   !!====================================================================== 
     2400 
    17252401END MODULE zdfosm 
Note: See TracChangeset for help on using the changeset viewer.