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

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

Ignore:
Timestamp:
2019-06-20T11:35:10+02:00 (5 years ago)
Author:
cguiavarch
Message:

Merge changes from /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser:r11139-11143

Location:
NEMO/branches/UKMO/NEMO_4.0_OSMOSIS
Files:
9 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/cfgs/SHARED/field_def_nemo-oce.xml

    r10823 r11145  
    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/UKMO/NEMO_4.0_OSMOSIS/src/OCE/C1D/step_c1d.F90

    r10888 r11145  
    5858 
    5959                             indic = 0                ! reset to no error condition 
    60       IF( kstp == nit000 )   CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 
     60                             IF( kstp == nit000 )   THEN 
     61                                CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 
     62                                CALL dia_hth_init    ! extra ML depth diagnostics, thermocline depths 
     63                             END IF 
    6164      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    6265                             CALL iom_setkt( kstp - nit000 + 1, "nemo" )   ! say to iom that we are at time step kstp 
    … …  
    8689      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    8790                         CALL dia_wri( kstp )       ! ocean model: outputs 
    88       IF( lk_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C) 
     91      IF( ll_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C) 
    8992 
    9093 

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/DIA/diahth.F90

    r10888 r11145  
    55   !!====================================================================== 
    66   !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code 
    7    !!                 !  1996-11  (E. Guilyardi)  OPA8  
     7   !!                 !  1996-11  (E. Guilyardi)  OPA8 
    88   !!                 !  1997-08  (G. Madec)  optimization 
    9    !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content  
     9   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
    1010   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    1111   !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
    12    !!---------------------------------------------------------------------- 
    13 #if   defined key_diahth 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_diahth' :                              thermocline depth diag. 
    1612   !!---------------------------------------------------------------------- 
    1713   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer 
    … …  
    2420   USE lib_mpp         ! MPP library 
    2521   USE iom             ! I/O library 
    26    USE timing          ! preformance summary 
    2722 
    2823   IMPLICIT NONE 
    … …  
    3025 
    3126   PUBLIC   dia_hth       ! routine called by step.F90 
    32    PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90 
    33  
    34    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag 
    35     
    36    ! note: following variables should move to local variables once iom_put is always used  
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m] 
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m] 
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W] 
     27   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90 
     28 
     29   LOGICAL, PUBLIC ::   ll_diahth   !: Compute further diagnostics of ML and thermocline depth 
    4130 
    4231   !!---------------------------------------------------------------------- 
    … …  
    4736CONTAINS 
    4837 
    49    FUNCTION dia_hth_alloc() 
    50       !!--------------------------------------------------------------------- 
    51       INTEGER :: dia_hth_alloc 
    52       !!--------------------------------------------------------------------- 
     38 
     39   SUBROUTINE dia_hth( kt ) 
     40     !!--------------------------------------------------------------------- 
     41     !!                  ***  ROUTINE dia_hth  *** 
     42     !! 
     43     !! ** Purpose : Computes 
     44     !!      the mixing layer depth (turbocline): avt = 5.e-4 
     45     !!      the depth of strongest vertical temperature gradient 
     46     !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
     47     !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 
     48     !!      the top of the thermochine: tn = tn(10m) - ztem2 
     49     !!      the pycnocline depth with density criteria equivalent to a temperature variation 
     50     !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     51     !!      the barrier layer thickness 
     52     !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
     53     !!      the depth of the 20 degree isotherm (linear interpolation) 
     54     !!      the depth of the 28 degree isotherm (linear interpolation) 
     55     !!      the heat content of first 300 m 
     56     !! 
     57     !! ** Method : 
     58     !!------------------------------------------------------------------- 
     59     INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     60     !! 
     61     INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
     62     INTEGER                          ::   iid, ilevel           ! temporary integers 
     63     INTEGER, DIMENSION(jpi,jpj) ::   ik20, ik28  ! levels 
     64     REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
     65     REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
     66     REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
     67     REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
     68     REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
     69     REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
     70     REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
     71     REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
     72     REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2 
     73     REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3 
     74     REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     75     REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion 
     76     REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion 
     77     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
     78     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
     79     REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz 
     80     REAL(wp), DIMENSION(jpi,jpj) ::   zthick     ! vertical integration thickness 
     81     REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
     82   ! note: following variables should move to local variables once iom_put is always used 
     83     REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m] 
     84     REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m] 
     85     REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m] 
     86     REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W] 
     87 
     88     IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN 
     89        ! ------------------------------------------------------------- ! 
     90        ! thermocline depth: strongest vertical gradient of temperature ! 
     91        ! turbocline depth (mixing layer depth): avt = zavt5            ! 
     92        ! MLD: rho = rho(1) + zrho3                                     ! 
     93        ! MLD: rho = rho(1) + zrho1                                     ! 
     94        ! ------------------------------------------------------------- ! 
     95        zmaxdzT(:,:) = 0._wp 
     96        IF( nla10 > 1 ) THEN 
     97           DO jj = 1, jpj 
     98              DO ji = 1, jpi 
     99                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     100                 zrho0_3(ji,jj) = zztmp 
     101                 zrho0_1(ji,jj) = zztmp 
     102                 zhth(ji,jj) = zztmp 
     103              END DO 
     104           END DO 
     105        ELSE IF (iom_use("mlddzt")) THEN 
     106           DO jj = 1, jpj 
     107              DO ji = 1, jpi 
     108                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     109                 zhth(ji,jj) = zztmp 
     110              END DO 
     111           END DO 
     112        ELSE 
     113           zhth(:,:) = 0._wp 
     114 
     115        ENDIF 
     116 
     117        DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
     118           DO jj = 1, jpj 
     119              DO ji = 1, jpi 
     120                 ! 
     121                 zzdep = gdepw_n(ji,jj,jk) 
     122                 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     123                 zzdep = zzdep * tmask(ji,jj,1) 
     124 
     125                 IF( zztmp > zmaxdzT(ji,jj) ) THEN 
     126                    zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     127                 ENDIF 
     128 
     129                 IF( nla10 > 1 ) THEN 
     130                    zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
     131                    IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
     132                    IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
     133                 ENDIF 
     134 
     135              END DO 
     136           END DO 
     137        END DO 
     138 
     139        IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline 
     140        IF( nla10 > 1 ) THEN 
     141           IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03 
     142           IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01 
     143        ENDIF 
     144     ENDIF 
     145 
     146     IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
     147          & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN 
     148        DO jj = 1, jpj 
     149           DO ji = 1, jpi 
     150              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     151              zabs2   (ji,jj) = zztmp 
     152              ztm2    (ji,jj) = zztmp 
     153              zrho10_3(ji,jj) = zztmp 
     154              zpycn   (ji,jj) = zztmp 
     155           END DO 
     156        END DO 
     157        ztinv  (:,:) = 0._wp 
     158        zdepinv(:,:) = 0._wp 
     159 
     160        IF (iom_use("pycndep")) THEN 
     161           ! Preliminary computation 
     162           ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
     163           DO jj = 1, jpj 
     164              DO ji = 1, jpi 
     165                 IF( tmask(ji,jj,nla10) == 1. ) THEN 
     166                    zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
     167                         &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
     168                         &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
     169                    zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
     170                         &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
     171                    zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
     172                    zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
     173                    zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     174                    zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     175                 ELSE 
     176                    zdelr(ji,jj) = 0._wp 
     177                 ENDIF 
     178              END DO 
     179           END DO 
     180        ELSE 
     181           zdelr(:,:) = 0._wp 
     182        ENDIF 
     183 
     184        ! ------------------------------------------------------------- ! 
     185        ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
     186        ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
     187        ! MLD: rho = rho10m + zrho3                                     ! 
     188        ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
     189        ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
     190        ! depth of temperature inversion                                ! 
     191        ! ------------------------------------------------------------- ! 
     192        DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
     193           DO jj = 1, jpj 
     194              DO ji = 1, jpi 
     195                 ! 
     196                 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
     197                 ! 
     198                 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
     199                 IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     200                 IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     201                 zztmp = -zztmp                                          ! delta T(10m) 
     202                 IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     203                    ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
     204                 ENDIF 
     205 
     206                 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     207                 IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     208                 IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     209                 ! 
     210              END DO 
     211           END DO 
     212        END DO 
     213 
     214        IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2 
     215        IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2 
     216        IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03 
     217        IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2 
     218        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref) 
     219        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref) 
     220     ENDIF 
     221 
     222     IF(iom_use("20d") .OR. iom_use("28d")) THEN 
     223        ! ----------------------------------- ! 
     224        ! search deepest level above 20C/28C  ! 
     225        ! ----------------------------------- ! 
     226        ik20(:,:) = 1 
     227        ik28(:,:) = 1 
     228        DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
     229           DO jj = 1, jpj 
     230              DO ji = 1, jpi 
     231                 zztmp = tsn(ji,jj,jk,jp_tem) 
     232                 IF( zztmp >= 20. )   ik20(ji,jj) = jk 
     233                 IF( zztmp >= 28. )   ik28(ji,jj) = jk 
     234              END DO 
     235           END DO 
     236        END DO 
     237 
     238        ! --------------------------- ! 
     239        !  Depth of 20C/28C isotherm  ! 
     240        ! --------------------------- ! 
     241        DO jj = 1, jpj 
     242           DO ji = 1, jpi 
     243              ! 
     244              zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
     245              ! 
     246              iid = ik20(ji,jj) 
     247              IF( iid /= 1 ) THEN 
     248                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     249                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     250                      &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
     251                      &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     252                 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
     253              ELSE 
     254                 zhd20(ji,jj) = 0._wp 
     255              ENDIF 
     256              ! 
     257              iid = ik28(ji,jj) 
     258              IF( iid /= 1 ) THEN 
     259                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     260                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     261                      &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
     262                      &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     263                 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
     264              ELSE 
     265                 zhd28(ji,jj) = 0._wp 
     266              ENDIF 
     267 
     268           END DO 
     269        END DO 
     270        CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm 
     271        CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm 
     272     ENDIF 
     273 
     274     ! ----------------------------- ! 
     275     !  Heat content of first 300 m  ! 
     276     ! ----------------------------- ! 
     277     IF (iom_use("hc300")) THEN 
     278        ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
     279        ilevel   = 0 
     280        zthick_0 = 0._wp 
     281        DO jk = 1, jpkm1 
     282           zthick_0 = zthick_0 + e3t_1d(jk) 
     283           IF( zthick_0 < 300. )   ilevel = jk 
     284        END DO 
     285        ! surface boundary condition 
     286        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 
     287        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp 
     288        ENDIF 
     289        ! integration down to ilevel 
     290        DO jk = 1, ilevel 
     291           zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
     292           zhtc3  (:,:) = zhtc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
     293        END DO 
     294        ! deepest layer 
     295        zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
     296        DO jj = 1, jpj 
     297           DO ji = 1, jpi 
     298              zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
     299                   &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
     300           END DO 
     301        END DO 
     302        ! from temperature to heat contain 
     303        zcoef = rau0 * rcp 
     304        zhtc3(:,:) = zcoef * zhtc3(:,:) 
     305        CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content 
     306     ENDIF 
     307     ! 
     308   END SUBROUTINE dia_hth 
     309 
     310 
     311   SUBROUTINE dia_hth_init 
     312      !!--------------------------------------------------------------------------- 
     313      !!                  ***  ROUTINE dia_hth_init  *** 
     314      !! 
     315      !! ** Purpose: Initialization for ML and thermocline depths 
     316      !! 
     317      !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth 
     318      !!--------------------------------------------------------------------------- 
     319       ! 
     320      IF(lwp) THEN 
     321         WRITE(numout,*) 
     322         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 
     323         WRITE(numout,*) '~~~~~~~~~~~~ ' 
     324      ENDIF 
     325      ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR.  & 
     326           & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
     327           & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. & 
     328           & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300") 
     329      IF(lwp) THEN 
     330         WRITE(numout,*) '   output upper ocean diagnostics (T) or not (F)       ll_diahth = ', ll_diahth 
     331      ENDIF 
    53332      ! 
    54       ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 
    55       ! 
    56       CALL mpp_sum ( 'diahth', dia_hth_alloc ) 
    57       IF(dia_hth_alloc /= 0)   CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' ) 
    58       ! 
    59    END FUNCTION dia_hth_alloc 
    60  
    61  
    62    SUBROUTINE dia_hth( kt ) 
    63       !!--------------------------------------------------------------------- 
    64       !!                  ***  ROUTINE dia_hth  *** 
    65       !! 
    66       !! ** Purpose : Computes 
    67       !!      the mixing layer depth (turbocline): avt = 5.e-4 
    68       !!      the depth of strongest vertical temperature gradient 
    69       !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
    70       !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2        
    71       !!      the top of the thermochine: tn = tn(10m) - ztem2  
    72       !!      the pycnocline depth with density criteria equivalent to a temperature variation  
    73       !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)  
    74       !!      the barrier layer thickness 
    75       !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
    76       !!      the depth of the 20 degree isotherm (linear interpolation) 
    77       !!      the depth of the 28 degree isotherm (linear interpolation) 
    78       !!      the heat content of first 300 m 
    79       !! 
    80       !! ** Method :  
    81       !!------------------------------------------------------------------- 
    82       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    83       !! 
    84       INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    85       INTEGER                          ::   iid, ilevel           ! temporary integers 
    86       INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels 
    87       REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    88       REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
    89       REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
    90       REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
    91       REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
    92       REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    93       REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
    94       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2  
    95       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2      
    96       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3       
    97       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    98       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion 
    99       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion 
    100       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
    101       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
    102       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz 
    103       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness  
    104       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
    105       !!---------------------------------------------------------------------- 
    106       IF( ln_timing )   CALL timing_start('dia_hth') 
    107  
    108       IF( kt == nit000 ) THEN 
    109          !                                      ! allocate dia_hth array 
    110          IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 
    111  
    112          IF(.NOT. ALLOCATED(ik20) ) THEN 
    113             ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 
    114                &      zabs2(jpi,jpj),   & 
    115                &      ztm2(jpi,jpj),    & 
    116                &      zrho10_3(jpi,jpj),& 
    117                &      zpycn(jpi,jpj),   & 
    118                &      ztinv(jpi,jpj),   & 
    119                &      zdepinv(jpi,jpj), & 
    120                &      zrho0_3(jpi,jpj), & 
    121                &      zrho0_1(jpi,jpj), & 
    122                &      zmaxdzT(jpi,jpj), & 
    123                &      zthick(jpi,jpj),  & 
    124                &      zdelr(jpi,jpj), STAT=ji) 
    125             CALL mpp_sum('diahth', ji) 
    126             IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 
    127          END IF 
    128  
    129          IF(lwp) WRITE(numout,*) 
    130          IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 
    131          IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    132          IF(lwp) WRITE(numout,*) 
    133       ENDIF 
    134  
    135       ! initialization 
    136       ztinv  (:,:) = 0._wp   
    137       zdepinv(:,:) = 0._wp   
    138       zmaxdzT(:,:) = 0._wp   
    139       DO jj = 1, jpj 
    140          DO ji = 1, jpi 
    141             zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    142             hth     (ji,jj) = zztmp 
    143             zabs2   (ji,jj) = zztmp 
    144             ztm2    (ji,jj) = zztmp 
    145             zrho10_3(ji,jj) = zztmp 
    146             zpycn   (ji,jj) = zztmp 
    147         END DO 
    148       END DO 
    149       IF( nla10 > 1 ) THEN  
    150          DO jj = 1, jpj 
    151             DO ji = 1, jpi 
    152                zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    153                zrho0_3(ji,jj) = zztmp 
    154                zrho0_1(ji,jj) = zztmp 
    155             END DO 
    156          END DO 
    157       ENDIF 
    158        
    159       ! Preliminary computation 
    160       ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    161       DO jj = 1, jpj 
    162          DO ji = 1, jpi 
    163             IF( tmask(ji,jj,nla10) == 1. ) THEN 
    164                zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
    165                   &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
    166                   &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
    167                zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
    168                   &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
    169                zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
    170                zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
    171                zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    172                zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    173             ELSE 
    174                zdelr(ji,jj) = 0._wp 
    175             ENDIF 
    176          END DO 
    177       END DO 
    178  
    179       ! ------------------------------------------------------------- ! 
    180       ! thermocline depth: strongest vertical gradient of temperature ! 
    181       ! turbocline depth (mixing layer depth): avt = zavt5            ! 
    182       ! MLD: rho = rho(1) + zrho3                                     ! 
    183       ! MLD: rho = rho(1) + zrho1                                     ! 
    184       ! ------------------------------------------------------------- ! 
    185       DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    186          DO jj = 1, jpj 
    187             DO ji = 1, jpi 
    188                ! 
    189                zzdep = gdepw_n(ji,jj,jk) 
    190                zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    191                zzdep = zzdep * tmask(ji,jj,1) 
    192  
    193                IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
    194                   zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    195                ENDIF 
    196                 
    197                IF( nla10 > 1 ) THEN  
    198                   zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
    199                   IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
    200                   IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
    201                ENDIF 
    202  
    203             END DO 
    204          END DO 
    205       END DO 
    206        
    207       CALL iom_put( "mlddzt", hth )            ! depth of the thermocline 
    208       IF( nla10 > 1 ) THEN  
    209          CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03 
    210          CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01 
    211       ENDIF 
    212  
    213       ! ------------------------------------------------------------- ! 
    214       ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
    215       ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
    216       ! MLD: rho = rho10m + zrho3                                     ! 
    217       ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
    218       ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
    219       ! depth of temperature inversion                                ! 
    220       ! ------------------------------------------------------------- ! 
    221       DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    222          DO jj = 1, jpj 
    223             DO ji = 1, jpi 
    224                ! 
    225                zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
    226                ! 
    227                zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
    228                IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    229                IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    230                zztmp = -zztmp                                          ! delta T(10m) 
    231                IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    232                   ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
    233                ENDIF 
    234  
    235                zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    236                IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    237                IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    238                ! 
    239             END DO 
    240          END DO 
    241       END DO 
    242  
    243       CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2 
    244       CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2 
    245       CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03 
    246       CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2 
    247       CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)  
    248       CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)  
    249  
    250  
    251       ! ----------------------------------- ! 
    252       ! search deepest level above 20C/28C  ! 
    253       ! ----------------------------------- ! 
    254       ik20(:,:) = 1 
    255       ik28(:,:) = 1 
    256       DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    257          DO jj = 1, jpj 
    258             DO ji = 1, jpi 
    259                zztmp = tsn(ji,jj,jk,jp_tem) 
    260                IF( zztmp >= 20. )   ik20(ji,jj) = jk 
    261                IF( zztmp >= 28. )   ik28(ji,jj) = jk 
    262             END DO 
    263          END DO 
    264       END DO 
    265  
    266       ! --------------------------- ! 
    267       !  Depth of 20C/28C isotherm  ! 
    268       ! --------------------------- ! 
    269       DO jj = 1, jpj 
    270          DO ji = 1, jpi 
    271             ! 
    272             zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
    273             ! 
    274             iid = ik20(ji,jj) 
    275             IF( iid /= 1 ) THEN  
    276                zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    277                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    278                   &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
    279                   &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    280                hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    281             ELSE  
    282                hd20(ji,jj) = 0._wp 
    283             ENDIF 
    284             ! 
    285             iid = ik28(ji,jj) 
    286             IF( iid /= 1 ) THEN  
    287                zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    288                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    289                   &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
    290                   &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    291                hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    292             ELSE  
    293                hd28(ji,jj) = 0._wp 
    294             ENDIF 
    295  
    296          END DO 
    297       END DO 
    298       CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm 
    299       CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm 
    300  
    301       ! ----------------------------- ! 
    302       !  Heat content of first 300 m  ! 
    303       ! ----------------------------- ! 
    304  
    305       ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
    306       ilevel   = 0 
    307       zthick_0 = 0._wp 
    308       DO jk = 1, jpkm1                       
    309          zthick_0 = zthick_0 + e3t_1d(jk) 
    310          IF( zthick_0 < 300. )   ilevel = jk 
    311       END DO 
    312       ! surface boundary condition 
    313       IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
    314       ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                    
    315       ENDIF 
    316       ! integration down to ilevel 
    317       DO jk = 1, ilevel 
    318          zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
    319          htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
    320       END DO 
    321       ! deepest layer 
    322       zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
    323       DO jj = 1, jpj 
    324          DO ji = 1, jpi 
    325             htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
    326                &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
    327          END DO 
    328       END DO 
    329       ! from temperature to heat contain 
    330       zcoef = rau0 * rcp 
    331       htc3(:,:) = zcoef * htc3(:,:) 
    332       CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
    333       ! 
    334       IF( ln_timing )   CALL timing_stop('dia_hth') 
    335       ! 
    336    END SUBROUTINE dia_hth 
    337  
    338 #else 
    339    !!---------------------------------------------------------------------- 
    340    !!   Default option :                                       Empty module 
    341    !!---------------------------------------------------------------------- 
    342    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag 
    343 CONTAINS 
    344    SUBROUTINE dia_hth( kt )         ! Empty routine 
    345       IMPLICIT NONE 
    346       INTEGER, INTENT( in ) :: kt 
    347       WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 
    348    END SUBROUTINE dia_hth 
    349 #endif 
    350  
    351    !!====================================================================== 
     333   END SUBROUTINE dia_hth_init 
    352334END MODULE diahth 

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/DIA/diawri.F90

    r10888 r11145  
    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/UKMO/NEMO_4.0_OSMOSIS/src/OCE/TRD/trd_oce.F90

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

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/TRD/trdtra.F90

    r10888 r11145  
    346346         CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl"  , ptrdx )        ! bottom boundary layer 
    347347                                  CALL iom_put( "strd_bbl"  , ptrdy ) 
     348         CASE( jptra_osm  )   ;   CALL iom_put( "ttrd_osm"  , ptrdx )        ! OSMOSIS non-local forcing 
     349                                  CALL iom_put( "strd_osm"  , ptrdy ) 
    348350         CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc"  , ptrdx )        ! static instability mixing 
    349351                                  CALL iom_put( "strd_npc"  , ptrdy ) 

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/ZDF/zdfosm.F90

    r10888 r11145  
    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/UKMO/NEMO_4.0_OSMOSIS/src/OCE/step.F90

    r10888 r11145  
    101101      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 
    102102                             CALL iom_init(      cxios_context          )  ! for model grid (including passible AGRIF zoom) 
    103          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
     103                             IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
     104                             CALL dia_hth_init    ! extra ML depth diagnostics, thermocline depths 
    104105      ENDIF 
    105106      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    … …  
    212213      ! diagnostics and outputs 
    213214      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    214       IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats 
    215       IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
    216       IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
    217       IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
    218                          CALL dia_ar5 ( kstp )        ! ar5 diag 
     215      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats 
     216      IF( ln_diacfl  )   CALL dia_cfl( kstp )         ! Courant number diagnostics 
     217      IF( ll_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth) 
     218      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports 
     219                         CALL dia_ar5( kstp )         ! ar5 diag 
    219220      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    220221                         CALL dia_wri ( kstp )        ! ocean model: outputs 

  • NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/stpctl.F90

    r10888 r11145  
    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.