Changeset 11143


Ignore:
Timestamp:
2019-06-19T18:15:12+02:00 (18 months ago)
Author:
agn
Message:

version w/o FK or local changes for NEMO repo

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/cfgs/SHARED/field_def_nemo-oce.xml

    r10824 r11143  
    200200        <field id="zws0"                long_name="surface non-local salinity flux"          unit="psu m/s" /> 
    201201        <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
     202        <field id="dh"          long_name="OBL pycnocline thickness"                                     unit="m"        /> 
     203        <field id="hml"          long_name="OBL MLD thickness"                                     unit="m"        /> 
    202204        <field id="hbli"                long_name="initial boundary layer depth"             unit="m"       /> 
    203205        <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
     
    205207        <field id="zwstrc"              long_name="convective velocity scale"                unit="m/s"     /> 
    206208        <field id="zwstrl"              long_name="langmuir velocity scale"                  unit="m/s"     /> 
     209        <field id="zvstr"              long_name="mixed La/wind velocity"                        unit="m/s"     /> 
    207210        <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    208211        <field id="zhbl"                long_name="boundary layer depth"                     unit="m"       /> 
    209212        <field id="zhml"                long_name="mixed layer depth"                        unit="m"       /> 
     213        <field id="zla"                long_name="Langmuir #"                        unit="#"       /> 
    210214        <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2"                        unit="mW"      /> 
    211215        <field id="wind_wave_power"     long_name="U_s \dot  tau"                            unit="mW"      /> 
    212216        <field id="wind_power"          long_name="\rho  u*^3"                               unit="mW"      /> 
    213217 
     218        <field id="ttrd_osm"          long_name="T-trend from non-local terms"         unit=" deg/s"    grid_ref="grid_T_3D"     /> 
     219        <field id="strd_osm"          long_name="S-trend from non-local terms"         unit=" psu/s"        grid_ref="grid_T_3D"   /> 
     220 
    214221        <!-- extra OSMOSIS diagnostics --> 
     222        <field id="ibld"              long_name="boundary-layer depth index"                  unit="#" /> 
     223        <field id="imld"              long_name="mixed-layer depth index"                  unit="#" /> 
    215224        <field id="zwthav"              long_name="av turb flux of T in ml"                  unit="deg m/s" /> 
    216225        <field id="zt_ml"               long_name="av T in ml"                               unit="deg"     /> 
     226        <field id="zdt_bl"               long_name="dT in ml"                               unit="deg"     /> 
     227        <field id="zds_bl"               long_name="dS in ml"                               unit="g/kg"     /> 
     228        <field id="zdb_bl"               long_name="db in ml"                               unit="m/s^2"     /> 
     229        <field id="zdu_bl"               long_name="du in ml"                               unit="m/s"     /> 
     230        <field id="zdv_bl"               long_name="dv in ml"                               unit="m/s"     /> 
    217231        <field id="zwth_ent"            long_name="entrainment turb flux of T"               unit="deg m/s" /> 
     232        <field id="zws_ent"            long_name="entrainment turb flux of S"               unit="g/kg m/s" /> 
     233        <field id="zwb_ent"            long_name="entrainment turb flux of buoyancy"     unit="m^2/s^3" /> 
    218234        <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
    219235        <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
     236 
     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" /> 
    220244      </field_group> 
    221245 
     
    224248        <field id="ghams"       long_name="non-local salinity flux"                          unit="psu m/s" /> 
    225249        <field id="zdtdz_pyc"   long_name="Pycnocline temperature gradient"                  unit=" deg/m"  /> 
     250        <field id="zdsdz_pyc"   long_name="Pycnocline salinity gradient"                  unit=" g/kg/m"  /> 
     251        <field id="zdbdz_pyc"   long_name="Pycnocline buoyancy gradient"                  unit=" s^-2"  /> 
     252 
     253        <!-- extra OSMOSIS diagnostics --> 
     254        <field id="ghamu_00"       long_name="non-local u-momentum flux on T before Stokes"                       unit="m^2/s^2" /> 
     255        <field id="ghamv_00"       long_name="non-local v-momentum flux on T before Stokes"                       unit="m^2/s^2" /> 
     256   <field id="ghamu_0"       long_name="non-local u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 
     257        <field id="ghamu_f"       long_name="non-local u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     258        <field id="ghamv_f"       long_name="non-local v-momentum flux on T after Coriolis"                          unit="m^2/s^2" /> 
     259        <field id="ghamu_b"       long_name="non-local u-momentum flux on T before base"                       unit="m^2/s^2" /> 
     260        <field id="ghamv_b"       long_name="non-local v-momentum flux on T before base"                          unit="m^2/s^2" /> 
     261        <field id="ghamu_1"       long_name="non-local u-momentum flux on T"                       unit="m^2/s^2" /> 
     262        <field id="ghamv_1"       long_name="non-local v-momentum flux on T"                          unit="m^2/s^2" /> 
     263        <field id="zviscos"       long_name="OSMOSIS viscosity on T-points"                          unit="m/s^2" /> 
     264        <field id="zdudz_pyc"       long_name="pycnocline u-shear T-points"                          unit="m/s^2" /> 
     265        <field id="zdvdz_pyc"       long_name="pycnocline v-shear T-points"                          unit="m/s^2" /> 
    226266      </field_group> 
    227267 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90

    r10425 r11143  
    4343   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
    4444   USE zdfmxl         ! mixed layer 
     45   USE zdfosm         ! mixed layer 
    4546   ! 
    4647   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    927928         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    928929      ENDIF 
     930 
     931      IF( ln_zdfosm ) THEN 
     932         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)   )    ! now boundary-layer depth 
     933         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)    )    ! now mixed-layer depth 
     934         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask       )    ! w-level diffusion 
     935         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask       )    ! now w-level viscosity 
     936         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask       )    ! non-local t forcing 
     937         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask       )    ! non-local s forcing 
     938         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask       )    ! non-local u forcing 
     939         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask       )    ! non-local v forcing 
     940      ENDIF 
     941     !     ! CALL histwrite( id_i, "zla", kt, zla*tmask(:,:,1)   , jpi*jpj, idex)         ! now Langmuir # 
     942     !     ! CALL histwrite( id_i, "zvstr", kt, zvstr*tmask(:,:,1)   , jpi*jpj, idex)     ! now mixed velocity scale 
     943     !     ! CALL histwrite( id_i, "zustar", kt, zustar*tmask(:,:,1)   , jpi*jpj, idex)   ! now friction velocity scale 
     944     !     ! CALL histwrite( id_i, "zwstrl", kt, zwstrl*tmask(:,:,1)   , jpi*jpj, idex)   ! now Langmuir velocity scale 
     945     !     ! CALL histwrite( id_i, "zwstrc", kt, zwstrc*tmask(:,:,1)   , jpi*jpj, idex)   ! now convective velocity scale 
     946     !     ! CALL histwrite( id_i, "zwb_ent", kt, zwb_ent*tmask(:,:,1)   , jpi*jpj, idex) ! now upward turb buoyancy entrainment flux 
     947     !     ! CALL histwrite( id_i, "zdb_bl", kt, zdb_bl*tmask(:,:,1)   , jpi*jpj, idex)   ! now db at ml base 
    929948  
    930949#if defined key_si3 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRD/trd_oce.F90

    r10068 r11143  
    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/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRD/trdtra.F90

    r10425 r11143  
    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/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90

    r10425 r11143  
    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 
     
    5158   USE traqsr         ! details of solar radiation absorption 
    5259   USE zdfddm         ! double diffusion mixing (avs array) 
     60   USE zdfmxl         ! mixed layer depth 
    5361   USE iom            ! I/O library 
    5462   USE lib_mpp        ! MPP library 
     
    7785   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt 
    7886   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 
     87   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh   ! depth of pycnocline 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml  ! ML depth 
    8089   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift. 
    81  
    8290   !                      !!** Namelist  namzdf_osm  ** 
    8391   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la 
     
    97105 
    98106   !                                    !!! ** General constants  ** 
    99    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number 
     107   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
     108   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    100109   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    101110   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
     
    114123      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    115124      !!---------------------------------------------------------------------- 
    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 ) 
     125     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
     126          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     127          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     128!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ? 
     129!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     130!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     131!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
     132!     IF ( ln_osm_mle ) THEN 
     133!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 
     134!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 
     135!     ENDIF 
     136 
    119137     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
    120138     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     
    166184      REAL(wp) ::   zbeta, zthermal                                  ! 
    167185      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 
    168       REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       ! 
     186      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      ! 
     187 
    169188      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density 
    170189      INTEGER  ::   jm                          ! dummy loop indices 
     
    191210      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    192211      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     212 
    193213      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift 
    194214      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number 
     
    196216      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 
    197217      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    198       LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 
     218      LOGICAL, DIMENSION(jpi,jpj)  :: lconv    ! unstable/stable bl 
    199219 
    200220      ! mixed-layer variables 
     
    202222      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 
    203223      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 
    204  
    205224      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients 
    206225      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear 
     
    210229      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    211230      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
     231      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure. 
     232      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients 
    212233      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    213234      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
     
    238259      ! Temporary variables 
    239260      INTEGER :: inhml 
    240       INTEGER :: i_lconv_alloc 
    241261      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 
    242262      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables 
     
    248268      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    249269 
     270      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
     271      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 
     272      REAL(wp) :: zzeta_s = 0._wp 
     273      REAL(wp) :: zzeta_v = 0.46 
     274      REAL(wp) :: zabsstke 
     275 
    250276      ! For debugging 
    251277      INTEGER :: ikt 
    252278      !!-------------------------------------------------------------------- 
    253279      ! 
    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  
    257280      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    258281      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
     
    268291      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    269292      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     293 
    270294      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    271295      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     
    277301      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    278302      ! 
     303      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 
    279304      ! Flux-Gradient arrays. 
    280305      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
     
    287312      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    288313 
     314      zdhdt_2(:,:) = 0._wp 
    289315      ! hbl = MAX(hbl,epsln) 
    290316      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    356382     CASE(2) 
    357383        zfac =  2.0_wp * rpi / 16.0_wp 
     384 
    358385        DO jj = 2, jpjm1 
    359386           DO ji = 2, jpim1 
     
    362389              ! It could represent the effects of the spread of wave directions 
    363390              ! 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 ! 
     391              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     392              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     393              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 ! 
    367394           END DO 
    368395        END DO 
     
    375402           ! Langmuir velocity scale (zwstrl), at T-point 
    376403           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 
     404           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 
     405           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 
    385406           ! Velocity scale that tends to zustar for large Langmuir numbers 
    386407           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
     
    389410           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    390411           ! Note zustke and zwstrl are not amended. 
    391            IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    392412           ! 
    393413           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
     
    406426     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 
    407427     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408      ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
    410       ibld(:,:) = 3 
    411       DO jk = 4, jpkm1 
     428     ! BL must be always 4 levels deep. 
     429      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
     430      ibld(:,:) = 4 
     431      DO jk = 5, jpkm1 
    412432         DO jj = 2, jpjm1 
    413433            DO ji = 2, jpim1 
     
    419439      END DO 
    420440 
    421       DO jj = 2, jpjm1                                 !  Vertical slab 
     441      DO jj = 2, jpjm1 
    422442         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 
     443            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     444            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 
     445            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     446         END DO 
    483447      END DO 
    484  
     448      ! Averages over well-mixed and boundary layer 
     449      ! Alan: do we need zb_nl?, zb_ml? 
     450      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     451      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     452! External gradient 
     453      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     454 
     455! Rate of change of hbl 
     456      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
    485457      ! Calculate averages over depth of boundary layer 
    486458      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 
     459      ibld(:,:) = 4 
     460 
     461         DO jj = 2, jpjm1 
     462            DO ji = 2, jpim1 
     463               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 
     464               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
     465               !zdhdt(ji,jj) = MIN(zdhdt(ji,jj), (zhbl_t(ji,jj) - hbl(ji,jj))/rn_rdt + wn(ji,jj,ibld(ji,jj))) 
     466            END DO 
     467         END DO 
     468      ! adjustment to represent limiting by ocean bottom 
    492469 
    493470      DO jk = 4, jpkm1 
     
    495472            DO ji = 2, jpim1 
    496473               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    497                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     474                  ibld(ji,jj) = jk 
    498475               ENDIF 
    499476            END DO 
     
    504481! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    505482! 
    506       DO jj = 2, jpjm1 
    507          DO ji = 2, jpim1 
    508             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     483      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     484      ! Alan: do we need zb_ml? 
     485      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    509486! 
    510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    511487! 
    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 
     488      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
    567489      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     ! 
     490! 
     491    ! Average over the depth of the mixed layer in the convective boundary layer 
     492    ! Alan: do we need zb_ml? 
     493    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    724494    ! rotate mean currents and changes onto wind align co-ordinates 
    725495    ! 
    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  
     496    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     497    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    745498     zuw_bse = 0._wp 
    746499     zvw_bse = 0._wp 
     500     zwth_ent = 0._wp 
     501     zws_ent = 0._wp 
    747502     DO jj = 2, jpjm1 
    748503        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) 
     504           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 
     505              IF ( lconv(ji,jj) ) THEN 
     506             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 
     507                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 
     508                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 
     509            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 
     510                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 
     511                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 
     512                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     513                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     514                 ENDIF 
     515              ELSE 
     516                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
     517                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    754518              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) ) 
    758519           ENDIF 
    759520        END DO 
     
    764525      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    765526 
    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 
     527      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     528      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 
     529      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    847530       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    848531       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    849532       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    850533 
    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 
    862534       DO jj = 2, jpjm1 
    863535          DO ji = 2, jpim1 
     
    868540               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    869541               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) 
     542               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 
    871543             ELSE 
    872544               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    873545               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 
     546             END IF 
     547          END DO 
     548       END DO 
    877549! 
    878550       DO jj = 2, jpjm1 
     
    889561                ! pycnocline - if present linear profile 
    890562                IF ( zdh(ji,jj) > 0._wp ) THEN 
     563                   zgamma_b = 6.0 
    891564                   DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892565                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    893566                       ! 
    894                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     567                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    895568                       ! 
    896                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     569                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    897570                   END DO 
     571                   IF ( ibld_ext == 0 ) THEN 
     572                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     573                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     574                   ELSE 
     575                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     576                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     577                   ENDIF 
    898578                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)) 
     579                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
     580                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)) 
     581                zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)) 
    901582                ! could be taken out, take account of entrainment represents as a diffusivity 
    902583                ! should remove w from here, represents entrainment 
     
    908589                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    909590                END DO 
     591 
     592                IF ( ibld_ext == 0 ) THEN 
     593                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     594                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     595                ELSE 
     596                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     597                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     598                ENDIF 
    910599             ENDIF   ! end if ( lconv ) 
    911 ! 
     600             ! 
    912601          END DO  ! end of ji loop 
    913602       END DO     ! end of jj loop 
     
    952641       END DO     ! end of jj loop 
    953642 
    954  
    955643! 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) 
    956644       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 ) 
     645          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 ) 
     646          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
     647          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 ) 
    960648       ELSEWHERE 
    961649          zsc_uw_1 = zustar**2 
    962           zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 
     650          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    963651       ENDWHERE 
    964  
     652       IF(ln_dia_osm) THEN 
     653          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     654          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     655       END IF 
    965656       DO jj = 2, jpjm1 
    966657          DO ji = 2, jpim1 
     
    1007698                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1008699                        &     * ( 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) 
     700                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
    1010701                   ! non-gradient buoyancy terms 
    1011702                   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 ) 
     
    1020711          END DO   ! ji loop 
    1021712       END DO      ! jj loop 
    1022  
    1023713 
    1024714       WHERE ( lconv ) 
     
    1051741       END DO           ! jj loop 
    1052742 
     743       IF(ln_dia_osm) THEN 
     744          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
     745          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     746       END IF 
    1053747! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    1054748 
     
    1089783       END DO      ! jj loop 
    1090784 
    1091  
    1092785       WHERE ( lconv ) 
    1093786          zsc_uw_1 = zustar**2 
     
    1134827          END DO 
    1135828       END DO 
     829 
     830       IF(ln_dia_osm) THEN 
     831          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
     832          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
     833          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     834          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     835          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     836          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     837       END IF 
    1136838! 
    1137839! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
     
    1165867      END DO 
    1166868 
     869       IF(ln_dia_osm) THEN 
     870          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     871          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     872       END IF 
    1167873      ! pynocline contributions 
    1168874       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
     
    1170876       DO jj = 2, jpjm1 
    1171877          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 
     878             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     879                DO jk= 2, ibld(ji,jj) 
     880                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     881                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     882                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     883                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     884                   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) 
     885                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     886                END DO 
     887             END IF 
     888          END DO 
    1181889       END DO 
    1182890 
     
    1185893       DO jj=2, jpjm1 
    1186894          DO ji = 2, jpim1 
    1187              IF ( lconv(ji,jj) ) THEN 
     895            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 
    1188896               DO jk = 1, imld(ji,jj) - 1 
    1189897                  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 
     898                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
     899                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1192900                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1193901                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
     
    1195903               DO jk = imld(ji,jj), ibld(ji,jj) 
    1196904                  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 ) 
     905                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
     906                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1199907                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1200908                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1201909                END DO 
    1202910             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 
     911 
     912             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     913             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     914             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     915             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1207916          END DO       ! ji loop 
    1208917       END DO          ! jj loop 
    1209918 
    1210  
     919       IF(ln_dia_osm) THEN 
     920          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     921          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     922          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 
     923          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 
     924          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
     925          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
     926          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
     927       END IF 
    1211928       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1212929       ! Need to put in code for contributions that are applied explicitly to 
     
    1232949       IF(ln_dia_osm) THEN 
    1233950          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     951          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     952          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
    1234953       END IF 
    1235954 
     
    12871006 
    12881007       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
    1289        CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. ) 
     1008       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 
    12901009 
    12911010       ! GN 25/8: need to change tmask --> wmask 
     
    13191038        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    13201039        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 
     1040         &                  ghamu, 'U', -1. , ghamv, 'V', -1. ) 
     1041 
     1042      IF(ln_dia_osm) THEN 
    13241043         SELECT CASE (nn_osm_wave) 
    13251044         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1). 
     
    13301049         ! Stokes drift read in from sbcwave  (=2). 
    13311050         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 
     1051            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
     1052            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     1053            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period 
     1054            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height 
     1055            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 
     1056            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum 
     1057            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10 
    13341058            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
    13351059                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
     
    13421066         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0> 
    13431067         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 
     1068         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1069         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
     1070         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
     1071         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base 
     1072         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base 
     1073         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base 
     1074         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth 
     1075         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth 
    13451076         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth 
    13461077         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points 
     
    13481079         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    13491080         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
     1081         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale 
     1082         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir # 
    13501083         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    13511084         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    13521085         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    13531086         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 
     1087         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1088         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    13551089         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 
     1090         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
     1091         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
     1092         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
     1093         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
     1094         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
    13591095      END IF 
    1360       ! Lateral boundary conditions on p_avt  (sign unchanged) 
    1361       CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 
     1096 
     1097CONTAINS 
     1098 
     1099 
     1100   ! Alan: do we need zb? 
     1101   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1102     !!--------------------------------------------------------------------- 
     1103     !!                   ***  ROUTINE zdf_vertical_average  *** 
     1104     !! 
     1105     !! ** Purpose : Determines vertical averages from surface to jnlev. 
     1106     !! 
     1107     !! ** Method  : Averages are calculated from the surface to jnlev. 
     1108     !!              The external level used to calculate differences is ibld+ibld_ext 
     1109     !! 
     1110     !!---------------------------------------------------------------------- 
     1111 
     1112        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1113 
     1114        ! Alan: do we need zb? 
     1115        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity 
     1116        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
     1117        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
     1118        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
     1119 
     1120        INTEGER :: jk, ji, jj 
     1121        REAL(wp) :: zthick, zthermal, zbeta 
     1122 
     1123 
     1124        zt   = 0._wp 
     1125        zs   = 0._wp 
     1126        zu   = 0._wp 
     1127        zv   = 0._wp 
     1128        DO jj = 2, jpjm1                                 !  Vertical slab 
     1129         DO ji = 2, jpim1 
     1130            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1131            zbeta    = rab_n(ji,jj,1,jp_sal) 
     1132               ! average over depth of boundary layer 
     1133            zthick = epsln 
     1134            DO jk = 2, jnlev_av(ji,jj) 
     1135               zthick = zthick + e3t_n(ji,jj,jk) 
     1136               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 
     1137               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     1138               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) & 
     1139                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 
     1140                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 
     1141               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) & 
     1142                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 
     1143                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 
     1144            END DO 
     1145            zt(ji,jj) = zt(ji,jj) / zthick 
     1146            zs(ji,jj) = zs(ji,jj) / zthick 
     1147            zu(ji,jj) = zu(ji,jj) / zthick 
     1148            zv(ji,jj) = zv(ji,jj) / zthick 
     1149            ! Alan: do we need zb? 
     1150            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 
     1151            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 
     1152            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 
     1153                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 
     1154            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 
     1155                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 
     1156            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1157         END DO 
     1158        END DO 
     1159   END SUBROUTINE zdf_osm_vertical_average 
     1160 
     1161   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 
     1162     !!--------------------------------------------------------------------- 
     1163     !!                   ***  ROUTINE zdf_velocity_rotation  *** 
     1164     !! 
     1165     !! ** Purpose : Rotates frame of reference of averaged velocity components. 
     1166     !! 
     1167     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w 
     1168     !! 
     1169     !!---------------------------------------------------------------------- 
     1170 
     1171        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle 
     1172        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current 
     1173        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline 
     1174 
     1175        INTEGER :: ji, jj 
     1176        REAL(wp) :: ztemp 
     1177 
     1178        DO jj = 2, jpjm1 
     1179           DO ji = 2, jpim1 
     1180              ztemp = zu(ji,jj) 
     1181              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 
     1182              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1183              ztemp = zdu(ji,jj) 
     1184              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
     1185              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1186           END DO 
     1187        END DO 
     1188    END SUBROUTINE zdf_osm_velocity_rotation 
     1189 
     1190    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 
     1191     !!--------------------------------------------------------------------- 
     1192     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     1193     !! 
     1194     !! ** Purpose : Calculates the gradients below the OSBL 
     1195     !! 
     1196     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1197     !! 
     1198     !!---------------------------------------------------------------------- 
     1199 
     1200     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
     1201 
     1202     INTEGER :: jj, ji, jkb, jkb1 
     1203     REAL(wp) :: zthermal, zbeta 
     1204 
     1205 
     1206     DO jj = 2, jpjm1 
     1207        DO ji = 2, jpim1 
     1208           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1209              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1210              zbeta    = rab_n(ji,jj,1,jp_sal) 
     1211              jkb = ibld(ji,jj) + ibld_ext 
     1212              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     1213              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 
     1214                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1215              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 
     1216                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1217              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1218           END IF 
     1219        END DO 
     1220     END DO 
     1221    END SUBROUTINE zdf_osm_external_gradients 
     1222 
     1223    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 
     1224 
     1225     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1226 
     1227     INTEGER :: jk, jj, ji 
     1228     REAL(wp) :: ztgrad, zsgrad, zbgrad 
     1229     REAL(wp) :: zgamma_b_nd, zgamma_c, znd 
     1230     REAL(wp) :: zzeta_s=0.3 
     1231 
     1232     DO jj = 2, jpjm1 
     1233        DO ji = 2, jpim1 
     1234           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1235              IF ( lconv(ji,jj) ) THEN  ! convective conditions 
     1236                 IF ( zdb_bl(ji,jj) > 0._wp .AND. zdbdz_ext(ji,jj) > 0._wp ) THEN 
     1237                    ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj) 
     1238                    zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj) 
     1239                    zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj) 
     1240                    zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1241                    zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 
     1242                         & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 
     1243                    DO jk = 2, ibld(ji,jj)+ibld_ext 
     1244                       znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     1245                       IF ( znd <= zzeta_s ) THEN 
     1246                          zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * & 
     1247                               & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1248                          zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * & 
     1249                               & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1250                          zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * & 
     1251                               & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1252                       ELSE 
     1253                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1254                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1255                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1256                       ENDIF 
     1257                    END DO 
     1258                 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 
     1259              ELSE 
     1260                 ! stable conditions 
     1261                 ! if pycnocline profile only defined when depth steady of increasing. 
     1262                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     1263                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1264                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     1265                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
     1266                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
     1267                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
     1268                          DO jk = 2, ibld(ji,jj) 
     1269                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1270                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1271                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1272                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1273                          END DO 
     1274                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1275                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
     1276                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
     1277                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
     1278                          DO jk = 2, ibld(ji,jj) 
     1279                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     1280                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1281                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1282                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1283                          END DO 
     1284                       ENDIF ! IF (zhol >=0.5) 
     1285                    ENDIF    ! IF (zdb_bl> 0.) 
     1286                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1287              ENDIF          ! IF (lconv) 
     1288           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1289        END DO 
     1290     END DO 
     1291 
     1292    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
     1293 
     1294    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
     1295      !!--------------------------------------------------------------------- 
     1296      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
     1297      !! 
     1298      !! ** Purpose : Calculates velocity shear in the pycnocline 
     1299      !! 
     1300      !! ** Method  : 
     1301      !! 
     1302      !!---------------------------------------------------------------------- 
     1303 
     1304      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
     1305 
     1306      INTEGER :: jk, jj, ji 
     1307      REAL(wp) :: zugrad, zvgrad, znd 
     1308      REAL(wp) :: zzeta_v = 0.45 
    13621309      ! 
    1363    END SUBROUTINE zdf_osm 
     1310      DO jj = 2, jpjm1 
     1311         DO ji = 2, jpim1 
     1312            ! 
     1313            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1314               IF ( lconv (ji,jj) ) THEN 
     1315                  ! Unstable conditions 
     1316                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1317                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1318                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1319                  !Alan is this right? 
     1320                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1321                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1322                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1323                       &      )/ (zdh(ji,jj)  + epsln ) 
     1324                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     1325                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     1326                     IF ( znd <= 0.0 ) THEN 
     1327                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     1328                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     1329                     ELSE 
     1330                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     1331                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     1332                     ENDIF 
     1333                  END DO 
     1334               ELSE 
     1335                  ! stable conditions 
     1336                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     1337                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     1338                  DO jk = 2, ibld(ji,jj) 
     1339                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1340                     IF ( znd < 1.0 ) THEN 
     1341                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     1342                     ELSE 
     1343                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     1344                     ENDIF 
     1345                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     1346                  END DO 
     1347               ENDIF 
     1348               ! 
     1349            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1350         END DO 
     1351      END DO 
     1352    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1353 
     1354    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
     1355     !!--------------------------------------------------------------------- 
     1356     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     1357     !! 
     1358     !! ** Purpose : Calculates the rate at which hbl changes. 
     1359     !! 
     1360     !! ** Method  : 
     1361     !! 
     1362     !!---------------------------------------------------------------------- 
     1363 
     1364    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl 
     1365 
     1366    INTEGER :: jj, ji 
     1367    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 
     1368    REAL(wp) :: zvel_max, zwb_min 
     1369    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1370    REAL(wp) :: zzeta_m = 0.3 
     1371    REAL(wp) :: zgamma_c = 2.0 
     1372    REAL(wp) :: zdhoh = 0.1 
     1373    REAL(wp) :: alpha_bc = 0.5 
     1374 
     1375    DO jj = 2, jpjm1 
     1376       DO ji = 2, jpim1 
     1377          IF ( lconv(ji,jj) ) THEN    ! Convective 
     1378             ! Alan is this right?  Yes, it's a bit different from the previous relationship 
     1379             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1380             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1381             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1382             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1383             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1384             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
     1385             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1386                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1387 
     1388             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1389                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1390                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1391                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1392             ! 
     1393             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 
     1394                  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) / & 
     1395                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1396                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1397                ! added ajgn 23 July as temporay fix 
     1398             zdhdt_2(ji,jj) = 0._wp 
     1399 
     1400                ! commented out ajgn 23 July as temporay fix 
     1401!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1402! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 
     1403!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1404!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 
     1405!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1406!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 
     1407!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 
     1408!                     ! Alan no idea what this should be? 
     1409!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 
     1410!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 
     1411!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 
     1412!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 
     1413!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 
     1414!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 
     1415!                 ENDIF 
     1416            ELSE                        ! Stable 
     1417                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     1418                zdhdt_2(ji,jj) = 0._wp 
     1419                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     1420                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     1421                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     1422                ELSE 
     1423                    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) ) 
     1424                ENDIF 
     1425                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 
     1426            ENDIF 
     1427         END DO 
     1428      END DO 
     1429    END SUBROUTINE zdf_osm_calculate_dhdt 
     1430 
     1431    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     1432     !!--------------------------------------------------------------------- 
     1433     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
     1434     !! 
     1435     !! ** Purpose : Increments hbl. 
     1436     !! 
     1437     !! ** Method  : If thechange in hbl exceeds one model level the change is 
     1438     !!              is calculated by moving down the grid, changing the buoyancy 
     1439     !!              jump. This is to ensure that the change in hbl does not 
     1440     !!              overshoot a stable layer. 
     1441     !! 
     1442     !!---------------------------------------------------------------------- 
     1443 
     1444 
     1445    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl. 
     1446 
     1447    INTEGER :: jk, jj, ji, jm 
     1448    REAL(wp) :: zhbl_s, zvel_max, zdb 
     1449    REAL(wp) :: zthermal, zbeta 
     1450 
     1451     DO jj = 2, jpjm1 
     1452         DO ji = 2, jpim1 
     1453           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     1454! 
     1455! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     1456! 
     1457              zhbl_s = hbl(ji,jj) 
     1458              jm = imld(ji,jj) 
     1459              zthermal = rab_n(ji,jj,1,jp_tem) 
     1460              zbeta = rab_n(ji,jj,1,jp_sal) 
     1461 
     1462 
     1463              IF ( lconv(ji,jj) ) THEN 
     1464!unstable 
     1465                    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) / & 
     1466                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1467                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1468                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
     1469                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 
     1470                         &  0.0 ) + zvel_max 
     1471 
     1472                      zhbl_s = zhbl_s + MIN( & 
     1473                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1474                        & e3w_n(ji,jj,jm) ) 
     1475                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1476 
     1477                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
     1478                 END DO 
     1479                 hbl(ji,jj) = zhbl_s 
     1480                 ibld(ji,jj) = jm 
     1481             ELSE 
     1482! stable 
     1483                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1484                    zdb = MAX( & 
     1485                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 
     1486                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 
     1487                         & 0.0 ) + & 
     1488             &       2.0 * zvstr(ji,jj)**2 / zhbl_s 
     1489 
     1490                    ! Alan is thuis right? I have simply changed hbli to hbl 
     1491                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 
     1492                    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) ) ) * & 
     1493               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 
     1494                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 
     1495                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 
     1496 
     1497                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1498                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
     1499                 END DO 
     1500             ENDIF   ! IF ( lconv ) 
     1501             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 
     1502             ibld(ji,jj) = MAX(jm, 4 ) 
     1503           ELSE 
     1504! change zero or one model level. 
     1505             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 
     1506           ENDIF 
     1507           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     1508         END DO 
     1509      END DO 
     1510 
     1511    END SUBROUTINE zdf_osm_timestep_hbl 
     1512 
     1513    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     1514      !!--------------------------------------------------------------------- 
     1515      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     1516      !! 
     1517      !! ** Purpose : Calculates thickness of the pycnocline 
     1518      !! 
     1519      !! ** Method  : The thickness is calculated from a prognostic equation 
     1520      !!              that relaxes the pycnocine thickness to a diagnostic 
     1521      !!              value. The time change is calculated assuming the 
     1522      !!              thickness relaxes exponentially. This is done to deal 
     1523      !!              with large timesteps. 
     1524      !! 
     1525      !!---------------------------------------------------------------------- 
     1526 
     1527      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
     1528      ! 
     1529      INTEGER :: jj, ji 
     1530      INTEGER :: inhml 
     1531      REAL(wp) :: zari, ztau, zddhdt 
     1532 
     1533 
     1534      DO jj = 2, jpjm1 
     1535         DO ji = 2, jpim1 
     1536            IF ( lconv(ji,jj) ) THEN 
     1537               IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     1538                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1539                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1540                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1541                  ELSE                                                     ! unstable 
     1542                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1543                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1544                  ENDIF 
     1545                  ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1546                  zddhdt = zari * hbl(ji,jj) 
     1547               ELSE 
     1548                  ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1549                  zddhdt = 0.2 * hbl(ji,jj) 
     1550               ENDIF 
     1551               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1552               ! Alan: this hml is never defined or used 
     1553            ELSE   ! IF (lconv) 
     1554               ztau = hbl(ji,jj) / zvstr(ji,jj) 
     1555               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     1556                  ! boundary layer deepening 
     1557                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1558                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     1559                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     1560                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     1561                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
     1562                  ELSE 
     1563                     zddhdt = 0.2 * hbl(ji,jj) 
     1564                  ENDIF 
     1565               ELSE     ! IF(dhdt < 0) 
     1566                  zddhdt = 0.2 * hbl(ji,jj) 
     1567               ENDIF    ! IF (dhdt >= 0) 
     1568               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1569               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 
     1570               ! Alan: this hml is never defined or used -- do we need it? 
     1571            ENDIF       ! IF (lconv) 
     1572 
     1573            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     1574            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     1575            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     1576            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     1577            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1578         END DO 
     1579      END DO 
     1580 
     1581    END SUBROUTINE zdf_osm_pycnocline_thickness 
     1582 
     1583END SUBROUTINE zdf_osm 
    13641584 
    13651585 
     
    13781598     INTEGER  ::   ios            ! local integer 
    13791599     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     1600     REAL z1_t2 
    13801601     !! 
    13811602     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    13821603          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
    1383           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 
     1604          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave 
    13841605     !!---------------------------------------------------------------------- 
    13851606     ! 
     
    13971618        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 
    13981619        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 
     1620        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters' 
     1621        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la 
    14011622        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
    14021623        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 
     1624        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
    14041625        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave 
    14051626        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave 
     
    14231644     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    14241645 
    1425      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     1646 
     1647     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 
     1648 
    14261649 
    14271650     IF( ln_zdfddm) THEN 
     
    15361759     REAL(wp) ::   zN2_c           ! local scalar 
    15371760     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) 
     1761     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 
    15391762     !!---------------------------------------------------------------------- 
    15401763     ! 
     
    15511774           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    15521775        END IF 
     1776 
    15531777        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    1554         id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. ) 
     1778        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. ) 
    15551779        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return 
    15561780           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' 
     1781           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  ) 
     1782           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file' 
    15591783           RETURN 
    1560         ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate 
     1784        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate 
    15611785           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 
    15621786        END IF 
     
    15701794         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
    15711795         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    1572          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
     1796         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios ) 
    15731797        RETURN 
    15741798     END IF 
     
    15781802     !!----------------------------------------------------------------------------- 
    15791803     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 
    1580      ALLOCATE( imld_rst(jpi,jpj) ) 
    15811804     ! w-level of the mixing and mixed layers 
    15821805     CALL eos_rab( tsn, rab_n ) 
     
    15991822     DO jj = 1, jpj 
    16001823        DO ji = 1, jpi 
    1601            iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
     1824           iiki = MAX(4,imld_rst(ji,jj)) 
     1825           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth 
     1826           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth 
    16031827        END DO 
    16041828     END DO 
    1605      hbl = MAX(hbl,epsln) 
    1606      hbli(:,:) = hbl(:,:) 
    1607      DEALLOCATE( imld_rst ) 
    16081829     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
     1830     wn(:,:,:) = 0._wp 
     1831     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    16091832   END SUBROUTINE osm_rst 
    16101833 
     
    16341857      ENDIF 
    16351858 
    1636       ! add non-local temperature and salinity flux 
    16371859      DO jk = 1, jpkm1 
    16381860         DO jj = 2, jpjm1 
     
    16481870      END DO 
    16491871 
    1650  
    1651       ! save the non-local tracer flux trends for diagnostic 
     1872      ! save the non-local tracer flux trends for diagnostics 
    16521873      IF( l_trdtra )   THEN 
    16531874         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    16541875         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 ) 
     1876 
     1877         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 
     1878         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 
    16581879         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    16591880      ENDIF 
     
    17231944 
    17241945   !!====================================================================== 
     1946 
    17251947END MODULE zdfosm 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/stpctl.F90

    r10570 r11143  
    6666      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
    6767      REAL(wp)               ::   zzz                 ! local real  
    68       REAL(wp), DIMENSION(9) ::   zmax 
    69       LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     68      REAL(wp), DIMENSION(10) ::   zmax 
     69      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_isnan 
    7070      CHARACTER(len=20) :: clname 
    7171      !!---------------------------------------------------------------------- 
     
    109109      ENDIF 
    110110      ! 
     111      ll_isnan = ANY(ISNAN(tsn)) .OR. ANY(ISNAN(un)) 
     112      IF (ll_isnan) nstop = nstop + 1 
    111113      !                                   !==  test of extrema  ==! 
    112114      IF( ll_wd ) THEN 
     
    124126         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
    125127         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max 
     128      ENDIF 
     129      IF (ll_isnan) THEN 
     130         zmax(10) = 1._wp                                           ! stop indicator 
     131      ELSE 
     132         zmax(10) = 0._wp 
    126133      ENDIF 
    127134      ! 
     
    147154      END IF 
    148155      !                                   !==  error handling  ==! 
    149       IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 
     156      IF( ( (ln_ctl .OR. sn_cfctl%l_runstat) .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 
    150157         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
    151158         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     
    153160         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
    154161         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
    155          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
    156          IF( lk_mpp .AND. ln_ctl ) THEN 
     162         &  zmax(10) >   0._wp  ) ) THEN   ! NaN encounter in the tests 
     163         ! &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
     164         IF( lk_mpp .AND. (ln_ctl .OR. sn_cfctl%l_runstat)) THEN 
    157165            CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih  ) 
    158166            CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  ) 
Note: See TracChangeset for help on using the changeset viewer.