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 14798 for NEMO/branches/2021 – NEMO

Changeset 14798 for NEMO/branches/2021


Ignore:
Timestamp:
2021-05-06T13:09:44+02:00 (3 years ago)
Author:
smueller
Message:

Upgrade of a range of local arrays to module arrays and various adjustments to improve compliance with coding conventions (ticket #2353)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90

    r14785 r14798  
    5656   !!      zdf_osm_mle_parameters               : timestep MLE depth and calculate MLE fluxes 
    5757   !!   zdf_osm_init   : initialization, namelist read, and parameters control 
     58   !!      zdf_osm_alloc                        : memory allocation 
    5859   !!   osm_rst        : read (or initialize) and write osmosis restart fields 
    5960   !!   tra_osm        : compute and add to the T & S trend the non-local flux 
    … …  
    6162   !!   dyn_osm        : compute and add to u & v trensd the non-local flux 
    6263   !!---------------------------------------------------------------------- 
    63    USE oce            ! ocean dynamics and active tracers 
    64                       ! uses ww from previous time step (which is now wb) to calculate hbl 
    65    USE dom_oce        ! ocean space and time domain 
    66    USE zdf_oce        ! ocean vertical physics 
    67    USE sbc_oce        ! surface boundary condition: ocean 
    68    USE sbcwave        ! surface wave parameters 
    69    USE phycst         ! physical constants 
    70    USE eosbn2         ! equation of state 
    71    USE traqsr         ! details of solar radiation absorption 
    72    USE zdfdrg, ONLY : rCdU_bot   ! bottom friction velocity 
    73    USE zdfddm         ! double diffusion mixing (avs array) 
    74    USE iom            ! I/O library 
    75    USE lib_mpp        ! MPP library 
    76    USE trd_oce        ! ocean trends definition 
    77    USE trdtra         ! tracers trends 
     64   USE oce                                        ! Ocean dynamics and active tracers 
     65                                                  ! Uses ww from previous time step (which is now wb) to calculate hbl 
     66   USE dom_oce                                    ! Ocean space and time domain 
     67   USE zdf_oce                                    ! Ocean vertical physics 
     68   USE sbc_oce                                    ! Surface boundary condition: ocean 
     69   USE sbcwave                                    ! Surface wave parameters 
     70   USE phycst                                     ! Physical constants 
     71   USE eosbn2                                     ! Equation of state 
     72   USE traqsr                                     ! Details of solar radiation absorption 
     73   USE zdfdrg, ONLY : rCdU_bot                    ! Bottom friction velocity 
     74   USE zdfddm                                     ! Double diffusion mixing (avs array) 
     75   USE iom                                        ! I/O library 
     76   USE lib_mpp                                    ! MPP library 
     77   USE trd_oce                                    ! Ocean trends definition 
     78   USE trdtra                                     ! Tracers trends 
     79   USE in_out_manager                             ! I/O manager 
     80   USE lbclnk                                     ! Ocean lateral boundary conditions (or mpp link) 
     81   USE prtctl                                     ! Print control 
     82   USE lib_fortran                                ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     83   USE timing, ONLY : timing_start, timing_stop   ! Timer 
    7884   ! 
    79    USE in_out_manager ! I/O manager 
    80    USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    81    USE prtctl         ! Print control 
    82    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    83    USE timing, ONLY : timing_start, timing_stop   ! Timer 
    84  
    8585   IMPLICIT NONE 
    8686   PRIVATE 
    87  
    88    PUBLIC   zdf_osm       ! routine called by step.F90 
    89    PUBLIC   zdf_osm_init  ! routine called by nemogcm.F90 
    90    PUBLIC   osm_rst       ! routine called by step.F90 
    91    PUBLIC   tra_osm       ! routine called by step.F90 
    92    PUBLIC   trc_osm       ! routine called by trcstp.F90 
    93    PUBLIC   dyn_osm       ! routine called by step.F90 
    94  
    95    PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90 
    96  
    97    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux 
    98    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux 
    99    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt    !: non-local temperature flux (gamma/<ws>o) 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams    !: non-local salinity flux (gamma/<ws>o) 
    101    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt 
    102    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth 
    103    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline 
    104    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth 
    105  
    106    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 
    107    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML 
    108    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML 
    109    INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer. 
    110  
     87   ! 
     88   ! Public subroutines 
     89   PUBLIC zdf_osm        ! Routine called by step.F90 
     90   PUBLIC zdf_osm_init   ! Routine called by nemogcm.F90 
     91   PUBLIC osm_rst        ! Routine called by step.F90 
     92   PUBLIC tra_osm        ! Routine called by step.F90 
     93   PUBLIC trc_osm        ! Routine called by trcstp.F90 
     94   PUBLIC dyn_osm        ! Routine called by step.F90 
     95   ! 
     96   ! Public variables 
     97   LOGICAL,  PUBLIC                                      ::   ln_osm_mle   !: Flag to activate the Mixed Layer Eddy (MLE) 
     98   !                                                                       !     parameterisation, needed by tra_mle_init in 
     99   !                                                                       !     tramle.F90 
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu        !: Non-local u-momentum flux 
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv        !: Non-local v-momentum flux 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt        !: Non-local temperature flux (gamma/<ws>o) 
     103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams        !: Non-local salinity flux (gamma/<ws>o) 
     104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl          !: Boundary layer depth 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml          !: ML depth 
     106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle         !: Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 
     107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle     !: Zonal buoyancy gradient in ML 
     108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle     !: Meridional buoyancy gradient in ML 
     109   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof     !: Level of base of MLE layer 
     110   ! 
    111111   INTERFACE zdf_osm_velocity_rotation 
    112112      !!--------------------------------------------------------------------- 
    … …  
    116116      MODULE PROCEDURE zdf_osm_velocity_rotation_3d 
    117117   END INTERFACE 
    118  
    119    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: r1_ft       ! inverse of the modified Coriolis parameter at t-pts 
    120  
     118   ! 
     119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   ! Averaging operator for avt 
     120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! Depth of pycnocline 
     121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_ft    ! Inverse of the modified Coriolis parameter at t-pts 
     122   ! 
     123   ! Layer indices 
     124   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nbld   ! Level of boundary layer base 
     125   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmld   ! Level of mixed-layer depth (pycnocline top) 
     126   ! 
     127   ! Layer type 
     128   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   n_ddh   ! Type of shear layer 
     129   !                                                        !    n_ddh=0: active shear layer 
     130   !                                                        !    n_ddh=1: shear layer not active 
     131   !                                                        !    n_ddh=2: shear production low 
     132   ! 
     133   ! Layer flags 
     134   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_conv    ! Unstable/stable bl 
     135   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_shear   ! Shear layers 
     136   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_coup    ! Coupling to bottom 
     137   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_pyc     ! OSBL pycnocline present 
     138   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_flux    ! Surface flux extends below OSBL into MLE layer 
     139   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_mle     ! MLE layer increases in hickness. 
     140   ! 
    121141   ! Scales 
    122    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swth0       ! Surface heat flux (Kinematic) 
    123    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sws0        ! Surface freshwater flux 
    124    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swb0        ! Surface buoyancy flux 
    125    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: suw0        ! Surface u-momentum flux 
    126    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustar      ! Friction velocity 
    127    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: scos_wind   ! Cos angle of surface stress 
    128    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: ssin_wind   ! Sin angle of surface stress 
    129    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swthav      ! Heat flux - bl average 
    130    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swsav       ! Freshwater flux - bl average 
    131    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swbav       ! Buoyancy flux - bl average 
    132    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustke      ! Surface Stokes drift 
    133    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: dstokes     ! Penetration depth of the Stokes drift 
    134    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrl      ! Langmuir velocity scale 
    135    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrc      ! Convective velocity scale 
    136    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sla         ! Trubulent Langmuir number 
    137    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: svstr       ! Velocity scale that tends to sustar for large Langmuir number 
    138    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: shol        ! Stability parameter for boundary layer 
    139  
    140    !                      !!** Namelist  namzdf_osm  ** 
    141    LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la 
    142  
    143    LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 
    144  
    145    REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    146    REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
    147    REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by 
    148    REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl 
    149    LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering 
    150    REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs 
    151    INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt 
    152    INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave 
    153    INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value 
    154    LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la 
    155    LOGICAL  ::   ln_dia_pyc_scl = .FALSE.   ! Output of pycnocline scalar-gradient profiles 
    156    LOGICAL  ::   ln_dia_pyc_shr = .FALSE.   ! Output of pycnocline velocity-shear  profiles 
    157  
    158  
    159    LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing 
    160    REAL(wp) ::   rn_riinfty   = 0.7     ! local Richardson Number limit for shear instability 
    161    REAL(wp) ::   rn_difri    =  0.005   ! maximum shear mixing at Rig = 0    (m2/s) 
    162    LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing 
    163    REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s) 
    164  
     142   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swth0       ! Surface heat flux (Kinematic) 
     143   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sws0        ! Surface freshwater flux 
     144   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swb0        ! Surface buoyancy flux 
     145   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   suw0        ! Surface u-momentum flux 
     146   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sustar      ! Friction velocity 
     147   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scos_wind   ! Cos angle of surface stress 
     148   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssin_wind   ! Sin angle of surface stress 
     149   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swthav      ! Heat flux - bl average 
     150   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swsav       ! Freshwater flux - bl average 
     151   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swbav       ! Buoyancy flux - bl average 
     152   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sustke      ! Surface Stokes drift 
     153   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dstokes     ! Penetration depth of the Stokes drift 
     154   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swstrl      ! Langmuir velocity scale 
     155   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swstrc      ! Convective velocity scale 
     156   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sla         ! Trubulent Langmuir number 
     157   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   svstr       ! Velocity scale that tends to sustar for large Langmuir number 
     158   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shol        ! Stability parameter for boundary layer 
     159   ! 
     160   !            ** Namelist  namzdf_osm  ** 
     161   LOGICAL  ::   ln_use_osm_la                      ! Use namelist rn_osm_la 
     162   REAL(wp) ::   rn_osm_la                          ! Turbulent Langmuir number 
     163   REAL(wp) ::   rn_osm_dstokes                     ! Depth scale of Stokes drift 
     164   REAL(wp) ::   rn_zdfosm_adjust_sd   = 1.0_wp     ! Factor to reduce Stokes drift by 
     165   REAL(wp) ::   rn_osm_hblfrac        = 0.1_wp     ! For nn_osm_wave = 3/4 specify fraction in top of hbl 
     166   LOGICAL  ::   ln_zdfosm_ice_shelter              ! Flag to activate ice sheltering 
     167   REAL(wp) ::   rn_osm_hbl0           = 10.0_wp    ! Initial value of hbl for 1D runs 
     168   INTEGER  ::   nn_ave                             ! = 0/1 flag for horizontal average on avt 
     169   INTEGER  ::   nn_osm_wave = 0                    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into 
     170   !                                                !    sbcwave 
     171   INTEGER  ::   nn_osm_SD_reduce                   ! = 0/1/2 flag for getting effective stokes drift from surface value 
     172   LOGICAL  ::   ln_dia_osm                         ! Use namelist  rn_osm_la 
     173   LOGICAL  ::   ln_dia_pyc_scl        = .FALSE.    ! Output of pycnocline scalar-gradient profiles 
     174   LOGICAL  ::   ln_dia_pyc_shr        = .FALSE.    ! Output of pycnocline velocity-shear  profiles 
     175   LOGICAL  ::   ln_kpprimix           = .TRUE.     ! Shear instability mixing 
     176   REAL(wp) ::   rn_riinfty            = 0.7_wp     ! Local Richardson Number limit for shear instability 
     177   REAL(wp) ::   rn_difri              = 0.005_wp   ! Maximum shear mixing at Rig = 0    (m2/s) 
     178   LOGICAL  ::   ln_convmix            = .TRUE.     ! Convective instability mixing 
     179   REAL(wp) ::   rn_difconv            = 1.0_wp     ! Diffusivity when unstable below BL  (m2/s) 
     180   ! 
    165181#ifdef key_osm_debug 
    166182   INTEGER :: nn_idb = 297, nn_jdb = 193, nn_kdb = 35, nn_narea_db = 109 
    167183   INTEGER :: iloc_db, jloc_db 
    168184#endif 
    169  
     185   ! 
    170186   ! OSMOSIS mixed layer eddy parametrization constants 
    171    INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt 
    172    REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient 
    173    !                                        ! parameters used in nn_osm_mle = 0 case 
    174    REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front 
    175    REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer 
    176    !                                        ! parameters used in nn_osm_mle = 1 case 
    177    REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
    178    LOGICAL  ::   ln_osm_hmle_limit           ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
    179    REAL(wp) ::   rn_osm_hmle_limit           ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
    180    REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
    181    REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
    182    REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 
    183    REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
    184    REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
    185    REAL(wp) ::   rn_osm_bl_thresh          ! Threshold buoyancy for deepening of OSBL base. 
    186    REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
    187  
    188  
    189    !                                    !!! ** General constants  ** 
    190    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
    191    REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    192    REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    193    REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
    194  
    195    INTEGER :: idebug = 236 
    196    INTEGER :: jdebug = 228 
    197  
     187   INTEGER  ::   nn_osm_mle          ! = 0/1 flag for horizontal average on avt 
     188   REAL(wp) ::   rn_osm_mle_ce       ! MLE coefficient 
     189   !   ! Parameters used in nn_osm_mle = 0 case 
     190   REAL(wp) ::   rn_osm_mle_lf       ! Typical scale of mixed layer front 
     191   REAL(wp) ::   rn_osm_mle_time     ! Time scale for mixing momentum across the mixed layer 
     192   !   ! Parameters used in nn_osm_mle = 1 case 
     193   REAL(wp) ::   rn_osm_mle_lat      ! Reference latitude for a 5 km scale of ML front 
     194   LOGICAL  ::   ln_osm_hmle_limit   ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
     195   REAL(wp) ::   rn_osm_hmle_limit   ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 
     196   REAL(wp) ::   rn_osm_mle_rho_c    ! Density criterion for definition of MLD used by FK 
     197   REAL(wp) ::   rb_c                ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 
     198   REAL(wp) ::   rc_f                ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
     199   REAL(wp) ::   rn_osm_mle_thresh   ! Threshold buoyancy for deepening of MLE layer below OSBL base 
     200   REAL(wp) ::   rn_osm_bl_thresh    ! Threshold buoyancy for deepening of OSBL base 
     201   REAL(wp) ::   rn_osm_mle_tau      ! Adjustment timescale for MLE 
     202   ! 
     203   !             ** General constants  ** 
     204   REAL(wp) ::   epsln     = 1.0e-20_wp      ! A small positive number to ensure no div by zero 
     205   REAL(wp) ::   depth_tol = 1.0e-6_wp       ! A small-ish positive number to give a hbl slightly shallower than gdepw 
     206   REAL(wp) ::   pthird    = 1.0_wp/3.0_wp   ! 1/3 
     207   REAL(wp) ::   p2third   = 2.0_wp/3.0_wp   ! 2/3 
     208   ! 
    198209   !! * Substitutions 
    199210#  include "do_loop_substitute.h90" 
    … …  
    205216   !!---------------------------------------------------------------------- 
    206217CONTAINS 
    207  
     218   ! 
    208219   INTEGER FUNCTION zdf_osm_alloc() 
    209220      !!---------------------------------------------------------------------- 
    210221      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    211222      !!---------------------------------------------------------------------- 
    212       ! 
    213223      INTEGER ::   ierr 
    214224      ! 
    215       ALLOCATE( swth0(jpi,jpj),  sws0(jpi,jpj),      swb0(jpi,jpj),      suw0(jpi,jpj),     & 
    216          &      sustar(jpi,jpj), scos_wind(jpi,jpj), ssin_wind(jpi,jpj), swthav(jpi,jpj),   & 
    217          &      swsav(jpi,jpj),  swbav(jpi,jpj),     sustke(jpi,jpj),    swstrl(jpi,jpj),   & 
    218          &      swstrc(jpi,jpj), sla(jpi,jpj),       svstr(jpi,jpj),     shol(jpi,jpj),     & 
    219          &      STAT=ierr ) 
    220       zdf_osm_alloc = ierr 
    221       ! 
    222       ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
    223          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
    224          &       etmean(jpi,jpj,jpk), STAT=ierr ) 
     225      zdf_osm_alloc = 0 
     226      ! 
     227      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), hbl(jpi,jpj), hml(jpi,jpj),   & 
     228         &      hmle(jpi,jpj),      dbdx_mle(jpi,jpj),  dbdy_mle(jpi,jpj),  mld_prof(jpi,jpj),  STAT=ierr ) 
    225229      zdf_osm_alloc = zdf_osm_alloc + ierr 
    226230      ! 
    227       ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 
    228          &       mld_prof(jpi,jpj), STAT=ierr ) 
     231      ALLOCATE( etmean(jpi,jpj,jpk), dh(jpi,jpj), r1_ft(jpi,jpj), STAT=ierr ) 
    229232      zdf_osm_alloc = zdf_osm_alloc + ierr 
    230233      ! 
     234      ALLOCATE( nbld(jpi,jpj), nmld(jpi,jpj), STAT=ierr ) 
     235      zdf_osm_alloc = zdf_osm_alloc + ierr 
     236      ! 
     237      ALLOCATE( n_ddh(jpi,jpj), STAT=ierr ) 
     238      zdf_osm_alloc = zdf_osm_alloc + ierr 
     239      ! 
     240      ALLOCATE( l_conv(jpi,jpj), l_shear(jpi,jpj), l_coup(jpi,jpj), l_pyc(jpi,jpj), l_flux(jpi,jpj), l_mle(jpi,jpj), STAT=ierr ) 
     241      zdf_osm_alloc = zdf_osm_alloc + ierr 
     242      ! 
     243      ALLOCATE( swth0(jpi,jpj),     sws0(jpi,jpj),   swb0(jpi,jpj),  suw0(jpi,jpj),  sustar(jpi,jpj), scos_wind(jpi,jpj),   & 
     244         &      ssin_wind(jpi,jpj), swthav(jpi,jpj), swsav(jpi,jpj), swbav(jpi,jpj), sustke(jpi,jpj), dstokes(jpi,jpj),     & 
     245         &      swstrl(jpi,jpj),    swstrc(jpi,jpj), sla(jpi,jpj),   svstr(jpi,jpj), shol(jpi,jpj),   STAT=ierr ) 
     246      zdf_osm_alloc = zdf_osm_alloc + ierr 
     247      ! 
    231248      CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
    232       IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
     249      IF( zdf_osm_alloc /= 0 ) CALL ctl_warn( 'zdf_osm_alloc: failed to allocate zdf_osm arrays' ) 
    233250      ! 
    234251   END FUNCTION zdf_osm_alloc 
    235  
    236  
    237    SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) 
     252   ! 
     253   SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm,   & 
     254      &                p_avt ) 
    238255      !!---------------------------------------------------------------------- 
    239256      !!                   ***  ROUTINE zdf_osm  *** 
    … …  
    305322      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK 
    306323 
    307       LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl 
    308       LOGICAL, DIMENSION(jpi,jpj)  :: lshear    ! Shear layers 
    309       LOGICAL, DIMENSION(jpi,jpj)  :: lcoup     ! Coupling to bottom 
    310       LOGICAL, DIMENSION(jpi,jpj)  :: lpyc      ! OSBL pycnocline present 
    311       LOGICAL, DIMENSION(jpi,jpj)  :: lflux     ! surface flux extends below OSBL into MLE layer. 
    312       LOGICAL, DIMENSION(jpi,jpj)  :: lmle      ! MLE layer increases in hickness. 
    313  
    314324      ! mixed-layer variables 
    315325 
    316       INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 
    317       INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 
    318326      INTEGER, DIMENSION(jpi,jpj) :: jp_ext ! offset for external level 
    319       INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 
    320327 
    321328      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    … …  
    373380      ! 
    374381      IF( ln_timing ) CALL timing_start('zdf_osm') 
    375       ibld(:,:)   = 0     ; imld(:,:)  = 0 
     382      nbld(:,:)   = 0      ; nmld(:,:)  = 0 
    376383      sustar(:,:) = zlarge 
    377384      swstrl(:,:) = zlarge ; svstr(:,:)      = zlarge ; swstrc(:,:)     = zlarge ; suw0(:,:)      = zlarge 
    … …  
    380387      sustke(:,:) = zlarge ; sla(:,:)        = zlarge ; scos_wind(:,:)  = zlarge ; ssin_wind(:,:) = zlarge 
    381388      shol(:,:)   = zlarge ; zalpha_pyc(:,:) = zlarge 
    382       lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE. 
     389      l_pyc(:,:) = .FALSE. ; l_flux(:,:) = .FALSE. ;  l_mle(:,:) = .FALSE. 
    383390      ! mixed layer 
    384391      ! no initialization of zhbl or zhml (or zdh?) 
    … …  
    607614         ! Note sustke and swstrl are not amended. 
    608615         ! 
    609          ! get convective velocity (swstrc), stabilty scale (shol) and logical conection flag lconv 
     616         ! get convective velocity (swstrc), stabilty scale (shol) and logical conection flag l_conv 
    610617         IF ( swbav(ji,jj) > 0.0) THEN 
    611618            swstrc(ji,jj) = ( 2.0 * swbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 
    … …  
    631638      ! BL must be always 4 levels deep. 
    632639      ! For calculation of lateral buoyancy gradients for FK in 
    633       ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 
     640      ! zdf_osm_zmld_horizontal_gradients need halo values for nbld, so must 
    634641      ! previously exist for hbl also. 
    635642 
    … …  
    637644      ! ########################################################################## 
    638645      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 
    639       ibld(:,:) = 4 
     646      nbld(:,:) = 4 
    640647      DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 
    641648         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    642             ibld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 
     649            nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 
    643650         ENDIF 
    644651      END_3D 
    … …  
    646653 
    647654      DO_2D( 0, 0, 0, 0 ) 
    648          zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    649          imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj) - 1, Kmm )) , 1 )) 
    650          zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
     655         zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) 
     656         nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) ) 
     657         zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm) 
    651658         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    652659      END_2D 
    … …  
    657664            & 'Before updating hbl: hbl=', hbl(ji,jj), ' dh=', dh(ji,jj), & 
    658665            &' zhbl =',zhbl(ji,jj) , ' zhml=', zhml(ji,jj), ' zdh=', zdh(ji,jj),& 
    659             &' imld=', imld(ji,jj), ' ibld=', ibld(ji,jj) 
     666            &' imld=', nmld(ji,jj), ' ibld=', nbld(ji,jj) 
    660667 
    661668         WRITE(narea+100,'(a,g11.3,a,2g11.3)') 'Physics: ssh ',ssh(ji,jj,Kmm),' T S surface=',ts(ji,jj,1,jp_tem,Kmm),ts(ji,jj,1,jp_sal,Kmm) 
    662          jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     669         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    663670         WRITE(narea+100,'(a,*(g11.3))') ' T[imld-1..ibld+2] =', ( ts(ji,jj,jk,jp_tem,Kmm), jk=jl,jm ) 
    664671         WRITE(narea+100,'(a,*(g11.3))') ' S[imld-1..ibld+2] =', ( ts(ji,jj,jk,jp_sal,Kmm), jk=jl,jm ) 
    … …  
    675682      ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 
    676683      jp_ext(:,:) = 1   ! ag 19/03 
    677       CALL zdf_osm_vertical_average( Kbb, Kmm,                                          & 
    678          &                           ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl,           & 
    679          &                           jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    680       jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    681       CALL zdf_osm_vertical_average( Kbb, Kmm,                                            & 
    682          &                           imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, jp_ext,   & 
    683          &                           zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
     684      CALL zdf_osm_vertical_average( Kbb,   Kmm,   nbld,  zt_bl,  zs_bl,   & 
     685         &                           zb_bl, zu_bl, zv_bl, jp_ext, zdt_bl,  & 
     686         &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     687      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     688      CALL zdf_osm_vertical_average( Kbb,    Kmm,    nmld - 1, zt_ml,  zs_ml,    & 
     689         &                           zb_ml,  zu_ml,  zv_ml,    jp_ext, zdt_ml,   & 
     690         &                           zds_ml, zdb_ml, zdu_ml,  zdv_ml ) 
    684691#ifdef key_osm_debug 
    685692      IF(narea==nn_narea_db) THEN 
    … …  
    713720 
    714721      ! Determine the state of the OSBL, stable/unstable, shear/no shear 
    715       CALL zdf_osm_osbl_state( Kmm, lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear,   & 
    716          &                     zhbl, zhml, zdh, zdb_bl, zdu_bl, zdv_bl, zdb_ml, zdu_ml, zdv_ml ) 
     722      CALL zdf_osm_osbl_state( Kmm,    zwb_ent, zwb_min, zshear, zhbl,     & 
     723         &                     zhml,   zdh,     zdb_bl,  zdu_bl, zdv_bl,   & 
     724         &                     zdb_ml, zdu_ml,  zdv_ml ) 
    717725 
    718726#ifdef key_osm_debug 
    … …  
    720728         ji=iloc_db; jj=jloc_db 
    721729         WRITE(narea+100,'(2(a,l7),a, i7,/,3(a,g11.3),/)') & 
    722             & 'After zdf_osm_osbl_state: lconv=', lconv(ji,jj), ' lshear=', lshear(ji,jj),  ' j_ddh=', j_ddh(ji,jj),& 
     730            & 'After zdf_osm_osbl_state: lconv=', l_conv(ji,jj), ' lshear=', l_shear(ji,jj),  ' j_ddh=', n_ddh(ji,jj),& 
    723731            & 'zwb_ent=', zwb_ent(ji,jj), ' zwb_min=', zwb_min(ji,jj),  ' zshear=', zshear(ji,jj) 
    724732         FLUSH(narea+100) 
    … …  
    731739            IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
    732740         END_3D 
    733          CALL zdf_osm_vertical_average( Kbb, Kmm,                                            & 
    734             &                           mld_prof, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle ) 
     741         CALL zdf_osm_vertical_average( Kbb,    Kmm,    mld_prof, zt_mle, zs_mle,   & 
     742            &                           zb_mle, zu_mle, zv_mle ) 
    735743 
    736744         DO_2D( 0, 0, 0, 0 ) 
    … …  
    749757 
    750758         !! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients 
    751          CALL zdf_osm_zmld_horizontal_gradients( Kmm, ibld, zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
     759         CALL zdf_osm_zmld_horizontal_gradients( Kmm,   zmld,     zdtdx,    zdtdy, zdsdx,   & 
     760            &                                    zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
    752761         !! Calculate max vertical FK flux zwb_fk & set logical descriptors 
    753          CALL zdf_osm_osbl_state_fk( Kmm, ibld, lconv, lpyc, lflux, lmle, zwb_fk, zt_mle, zs_mle, zb_mle, zhbl,   & 
    754             &                        zhmle, zwb_ent, zt_bl, zs_bl, zb_bl, zdb_bl, zdbds_mle ) 
     762         CALL zdf_osm_osbl_state_fk( Kmm,   zwb_fk, zt_mle,  zs_mle, zb_mle,   & 
     763            &                        zhbl,  zhmle,  zwb_ent, zt_bl,  zs_bl,    & 
     764            &                        zb_bl, zdb_bl, zdbds_mle ) 
    755765         !! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 
    756          CALL zdf_osm_mle_parameters( Kmm, lconv, lmle, mld_prof, zmld, zhmle,   & 
    757             &                         zvel_mle, zdiff_mle, zdbds_mle, zhbl, zb_bl, zwb0tot ) 
     766         CALL zdf_osm_mle_parameters( Kmm,       mld_prof,  zmld, zhmle, zvel_mle,   & 
     767            &                         zdiff_mle, zdbds_mle, zhbl, zb_bl, zwb0tot ) 
    758768#ifdef key_osm_debug 
    759769         IF(narea==nn_narea_db) THEN 
    … …  
    772782      ELSE    ! ln_osm_mle 
    773783         ! FK not selected, Boundary Layer only. 
    774          lpyc(:,:) = .TRUE. 
    775          lflux(:,:) = .FALSE. 
    776          lmle(:,:) = .FALSE. 
     784         l_pyc(:,:) = .TRUE. 
     785         l_flux(:,:) = .FALSE. 
     786         l_mle(:,:) = .FALSE. 
    777787         DO_2D( 0, 0, 0, 0 ) 
    778             IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 
     788            IF ( l_conv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
    779789         END_2D 
    780790      ENDIF   ! ln_osm_mle 
    781791 
    782792      !! External gradient below BL needed both with and w/o FK 
    783       CALL zdf_osm_external_gradients( Kmm, ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
     793      CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
    784794 
    785795      ! Test if pycnocline well resolved 
    786796!      DO_2D( 0, 0, 0, 0 )                                         Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 
    787 !         IF (lconv(ji,jj) ) THEN                                  should account for this. 
    788 !            ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) 
     797!         IF (l_conv(ji,jj) ) THEN                                  should account for this. 
     798!            ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm) 
    789799!            IF ( ztmp > 6 ) THEN 
    790800!               ! pycnocline well resolved 
    … …  
    803813         ji=iloc_db; jj=jloc_db 
    804814         WRITE(narea+100,'(4(a,l7),a,i7,/, 3(a,g11.3),/)') & 
    805             & 'BL logical descriptors: lconv =',lconv(ji,jj),' lpyc=', lpyc(ji,jj),' lflux=', lflux(ji,jj),' lmle=', lmle(ji,jj),& 
     815            & 'BL logical descriptors: lconv =',l_conv(ji,jj),' lpyc=', l_pyc(ji,jj),' lflux=', l_flux(ji,jj),' lmle=', l_mle(ji,jj),& 
    806816            & ' jp_ext=', jp_ext(ji,jj), & 
    807817            & 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj) 
    … …  
    812822      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 
    813823      jp_ext(:,:) = 1   ! ag 19/03 
    814       CALL zdf_osm_vertical_average( Kbb, Kmm,                                          & 
    815          &                           ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl,           & 
    816          &                           jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    817       jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    818       CALL zdf_osm_vertical_average( Kbb, Kmm,                                               & 
    819          &                           imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml,              & 
    820          &                           jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )   ! ag 19/03 
     824      CALL zdf_osm_vertical_average( Kbb,    Kmm,    nbld,   zt_bl,  zs_bl,    & 
     825         &                           zb_bl,  zu_bl,  zv_bl,  jp_ext, zdt_bl,   & 
     826         &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     827      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     828      CALL zdf_osm_vertical_average( Kbb,    Kmm,    nmld - 1, zt_ml,  zs_ml,    & 
     829         &                           zb_ml,  zu_ml,  zv_ml,    jp_ext, zdt_ml,   & 
     830         &                           zds_ml, zdb_ml, zdu_ml,  zdv_ml )   ! ag 19/03 
    821831#ifdef key_osm_debug 
    822832      IF(narea==nn_narea_db) THEN 
    … …  
    836846 
    837847! Rate of change of hbl 
    838       CALL zdf_osm_calculate_dhdt( lconv, lshear, j_ddh, zdhdt, zhbl, zdh, zwb_ent, zwb_min,   & 
    839          &                         zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk, zvel_mle ) 
     848      CALL zdf_osm_calculate_dhdt( zdhdt,  zhbl,         zdh,    zwb_ent,  zwb_min,   & 
     849         &                         zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk,    & 
     850         &                         zvel_mle ) 
    840851      ! Test if surface boundary layer coupled to bottom 
    841       lcoup(:,:) = .FALSE.   ! ag 19/03 
     852      l_coup(:,:) = .FALSE.   ! ag 19/03 
    842853      DO_2D( 0, 0, 0, 0 ) 
    843          zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it 
     854         zhbl_t(ji,jj) = hbl(ji,jj) + ( zdhdt(ji,jj) - ww(ji,jj,nbld(ji,jj)) ) * rn_Dt   ! Certainly need ww here, so subtract it 
    844855         ! adjustment to represent limiting by ocean bottom 
    845856         IF ( mbkt(ji,jj) > 2 ) THEN   ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access 
    846857            IF ( zhbl_t(ji,jj) > gdepw(ji, jj,mbkt(ji,jj)-2,Kmm) ) THEN 
    847858               zhbl_t(ji,jj) = MIN( zhbl_t(ji,jj), gdepw(ji,jj,mbkt(ji,jj)-2,Kmm) )   ! ht(:,:)) 
    848                lpyc(ji,jj) = .FALSE. 
    849                lcoup(ji,jj) = .TRUE.   ! ag 19/03 
     859               l_pyc(ji,jj) = .FALSE. 
     860               l_coup(ji,jj) = .TRUE.   ! ag 19/03 
    850861            END IF 
    851862         END IF 
    … …  
    853864         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    854865            WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3)),2(a,l7)')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),& 
    855                & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_Dt,' delta hbl from w ', ww(ji,jj,ibld(ji,jj))*rn_Dt,   & 
    856                & ' lcoup= ', lcoup(ji,jj), ' lpyc= ', lpyc(ji,jj) 
     866               & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_Dt,' delta hbl from w ', ww(ji,jj,nbld(ji,jj))*rn_Dt,   & 
     867               & ' lcoup= ', l_coup(ji,jj), ' lpyc= ', l_pyc(ji,jj) 
    857868            FLUSH(narea+100) 
    858869         END IF 
    … …  
    860871      END_2D 
    861872 
    862       imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index 
    863       ibld(:,:) = 4 
     873      nmld(:,:) = nbld(:,:)           ! use nmld to hold previous blayer index 
     874      nbld(:,:) = 4 
    864875 
    865876      DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 
    866877         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    867             ibld(ji,jj) = jk 
     878            nbld(ji,jj) = jk 
    868879         ENDIF 
    869880      END_3D 
    … …  
    872883      ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    873884      ! 
    874       CALL zdf_osm_timestep_hbl( Kmm, ibld, imld, lconv, lpyc, zdhdt, zhbl, zhbl_t, zt_bl, zs_bl, zwb_ent, zwb_fk_b ) 
     885      CALL zdf_osm_timestep_hbl( Kmm,   zdhdt,   zhbl, zhbl_t, zt_bl,   & 
     886         &                       zs_bl, zwb_ent, zwb_fk_b ) 
    875887      ! is external level in bounds? 
    876888 
    877889      !   Recalculate BL averages and differences using new BL depth 
    878890      jp_ext(:,:) = 1   ! ag 19/03 
    879       CALL zdf_osm_vertical_average( Kbb, Kmm,                                          & 
    880          &                           ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl,           & 
    881          &                           jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    882  
    883       CALL zdf_osm_pycnocline_thickness( Kmm, ibld, imld, lshear, lconv, j_ddh, zdh, zhml, zdhdt, zdb_bl,   & 
     891      CALL zdf_osm_vertical_average( Kbb,    Kmm,    nbld,   zt_bl,  zs_bl,    & 
     892         &                           zb_bl,  zu_bl,  zv_bl,  jp_ext, zdt_bl,   & 
     893         &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     894 
     895      CALL zdf_osm_pycnocline_thickness( Kmm,  zdh,     zhml,        zdhdt, zdb_bl,   & 
    884896         &                               zhbl, zwb_ent, zdbdz_bl_ext, zwb_fk_b ) 
    885897 
    … …  
    887899 
    888900      DO_2D( 0, 0, 0, 0 ) 
    889          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) >= mbkt(ji,jj) - 2 .or.   & 
    890             & ibld(ji,jj)-imld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0.0_wp ) THEN              ! ag 19/03 
    891             lpyc(ji,jj) = .FALSE.                           ! ag 19/03 
    892             IF ( ibld(ji,jj) >= mbkt(ji,jj) -2 ) THEN 
    893                imld(ji,jj) = ibld(ji,jj) - 1                ! ag 19/03 
    894                zdh(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) - gdepw(ji,jj,imld(ji,jj),Kmm)   ! ag 19/03 
    895                zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)   ! ag 19/03 
     901         IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. nbld(ji,jj) >= mbkt(ji,jj) - 2 .or.   & 
     902            & nbld(ji,jj) - nmld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0.0_wp ) THEN              ! ag 19/03 
     903            l_pyc(ji,jj) = .FALSE.                           ! ag 19/03 
     904            IF ( nbld(ji,jj) >= mbkt(ji,jj) -2 ) THEN 
     905               nmld(ji,jj) = nbld(ji,jj) - 1                ! ag 19/03 
     906               zdh(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) - gdepw(ji,jj,nmld(ji,jj),Kmm)   ! ag 19/03 
     907               zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm)   ! ag 19/03 
    896908               dh(ji,jj) = zdh(ji,jj)                       ! ag 19/03   
    897909               hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)          ! ag 19/03 
    … …  
    899911               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    900912                  WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= F5: ibld(ji,jj) >= mbkt(ji,jj) -2' 
    901                   WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj) 
     913                  WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',nbld(ji,jj),' imld=',nmld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj) 
    902914                  WRITE(narea+100,'(2(a,g11.3))')'dh=',dh(ji,jj),' hml=',hml(ji,jj) 
    903915                  FLUSH(narea+100) 
    … …  
    911923      ! 
    912924      ! Average over the depth of the mixed layer in the convective boundary layer 
    913       !      jp_ext = ibld - imld +1 
     925      !      jp_ext = nbld - nmld +1 
    914926      !     Recalculate ML averages and differences using new ML depth 
    915       jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    916       CALL zdf_osm_vertical_average( Kbb, Kmm,                                               & 
    917          &                           imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml,              & 
    918          &                           jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    919  
    920       CALL zdf_osm_external_gradients( Kmm, ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     927      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     928      CALL zdf_osm_vertical_average( Kbb,    Kmm,    nmld - 1, zt_ml,  zs_ml,    & 
     929         &                           zb_ml,  zu_ml,  zv_ml,    jp_ext, zdt_ml,   & 
     930         &                           zds_ml, zdb_ml, zdu_ml,  zdv_ml ) 
     931 
     932      CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    921933#ifdef key_osm_debug 
    922934      IF(narea==nn_narea_db) THEN 
    … …  
    953965 
    954966      jp_ext(:,:) = 1   ! ag 19/03 
    955       CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, ibld, jp_ext, lconv, lpyc, zdbdz_pyc, zalpha_pyc, zdh,   & 
    956          &                                       zdb_bl, zhbl, zdbdz_bl_ext, zhml, zdb_ml, zdhdt ) 
     967      CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm,    jp_ext, zdbdz_pyc,    zalpha_pyc, zdh,      & 
     968         &                                       zdb_bl, zhbl,   zdbdz_bl_ext, zhml,       zdb_ml,   & 
     969         &                                       zdhdt ) 
    957970#ifdef key_osm_debug 
    958971      IF(narea==nn_narea_db) THEN 
    959972         ji=iloc_db; jj=jloc_db 
    960          jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     973         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    961974         WRITE(narea+100,'(a,l7,/,3(a,g11.3),/)') & 
    962             & 'After pycnocline profiles BL  lpyc=', lpyc(ji,jj),& 
     975            & 'After pycnocline profiles BL  lpyc=', l_pyc(ji,jj),& 
    963976            & 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj), & 
    964977            & 'Pycnocline: zalpha_pyc=', zalpha_pyc(ji,jj) 
    … …  
    975988      ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    976989      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    977       CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, ibld, imld, lconv, lshear, lpyc, lcoup, j_ddh, zdiffut, zviscos,           & 
    978          &                                zhbl, zhml, zdh, zdhdt, zshear, zb_bl, zu_bl, zv_bl, zb_ml, zu_ml, zv_ml, zwb_ent,   & 
    979          &                                zwb_min ) 
     990      CALL zdf_osm_diffusivity_viscosity( Kbb,     Kmm,   zdiffut, zviscos, zhbl,    & 
     991         &                                zhml,    zdh,   zdhdt,   zshear,  zb_bl,   & 
     992         &                                zu_bl,   zv_bl, zb_ml,   zu_ml,   zv_ml,   & 
     993         &                                zwb_ent, zwb_min ) 
    980994#ifdef key_osm_debug 
    981995      IF(narea==nn_narea_db) THEN 
    982996         ji=iloc_db; jj=jloc_db 
    983          jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     997         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    984998         WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) 
    985999         WRITE(narea+100,'(a,*(g11.3))') ' zviscos[imld-1..ibld+2] =', ( zviscos(ji,jj,jk), jk=jl,jm ) 
    … …  
    9921006      ! Calculate non-gradient components of the flux-gradient relationships 
    9931007      ! -------------------------------------------------------------------- 
    994       CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zshear,          & 
    995          &                    zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml,       & 
    996          &                    zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 
     1008      CALL zdf_osm_fgr_terms( Kmm,        jp_ext,  zhbl,         zhml,         zdh,         & 
     1009         &                    zdhdt,      zshear,  zdt_bl,       zds_bl,       zdb_bl,      & 
     1010         &                    zdu_bl,     zdv_bl,  zdt_ml,       zds_ml,       zdb_ml,      & 
     1011         &                    zdu_ml,     zdv_ml,  zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc,   & 
     1012         &                    zalpha_pyc, zdiffut, zviscos ) 
    9971013 
    9981014      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    10041020 
    10051021      ! Rotate non-gradient velocity terms back to model reference frame 
    1006       CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE., 2, ibld ) 
     1022      CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE., 2, nbld ) 
    10071023 
    10081024      ! KPP-style Ri# mixing 
    … …  
    10101026         jkflt = jpk 
    10111027         DO_2D( 0, 0, 0, 0 ) 
    1012             IF ( ibld(ji,jj) < jkflt ) jkflt = ibld(ji,jj) 
     1028            IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj) 
    10131029         END_2D 
    10141030         DO jk = jkflt+1, jpkm1 
    … …  
    10231039            END_2D 
    10241040            DO_2D( 0, 0, 0, 0 ) 
    1025                IF ( jk > ibld(ji,jj) ) THEN 
     1041               IF ( jk > nbld(ji,jj) ) THEN 
    10261042                  ! Shear prod. at w-point weightened by mask 
    10271043                  zesh2  =  ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    … …  
    10421058      IF( ln_convmix) THEN 
    10431059         DO_2D( 0, 0, 0, 0 ) 
    1044             DO jk = ibld(ji,jj) + 1, jpkm1 
     1060            DO jk = nbld(ji,jj) + 1, jpkm1 
    10451061               IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) 
    10461062            END DO 
    … …  
    10501066      IF(narea==nn_narea_db) THEN 
    10511067         ji=iloc_db; jj=jloc_db 
    1052          jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     1068         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    10531069         WRITE(narea+100,'(a)') ' After including KPP Ri# diffusivity & viscosity' 
    10541070         WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) 
    … …  
    10631079      IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
    10641080         DO_2D( 0, 0, 0, 0 ) 
    1065             IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 
     1081            IF ( l_flux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 
    10661082               ! Calculate MLE flux contribution from surface fluxes 
    1067                DO jk = 1, ibld(ji,jj) 
     1083               DO jk = 1, nbld(ji,jj) 
    10681084                  znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) 
    10691085                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd ) 
    … …  
    10921108         IF(narea==nn_narea_db) THEN 
    10931109            ji=iloc_db; jj=jloc_db 
    1094             jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     1110            jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    10951111            WRITE(narea+100,'(a)') ' After including FK diffusivity & non-local terms' 
    10961112            WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) 
    … …  
    11341150      IF(narea==nn_narea_db) THEN 
    11351151         ji=iloc_db; jj=jloc_db 
    1136          jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     1152         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    11371153         WRITE(narea+100,'(a)') ' Final diffusivity & viscosity, & non-local terms' 
    11381154         WRITE(narea+100,'(a,*(g11.3))') ' p_avt[imld-1..ibld+2] =', ( p_avt(ji,jj,jk), jk=jl,jm ) 
    … …  
    11751191         IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*swth0 )     ! upward BL-avged turb buoyancy flux 
    11761192         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth 
    1177          IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1193         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )               ! boundary-layer max k 
    11781194         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
    11791195         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
    … …  
    11971213         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    11981214         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    1199          IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1215         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*nmld )               ! index for ML depth internal to zdf_osm routine 
    12001216         IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext )         ! =1 if pycnocline resolved internal to zdf_osm routine 
    1201          IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*j_ddh )            !    index forpyc thicknessh internal to zdf_osm routine 
     1217         IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*n_ddh )            !    index forpyc thicknessh internal to zdf_osm routine 
    12021218         IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear )         !    shear production of TKE internal to zdf_osm routine 
    12031219         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    … …  
    12251241   END SUBROUTINE zdf_osm 
    12261242 
    1227    SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm,                           & 
    1228       &                                 knlev, pt, ps, pb, pu, pv,          & 
    1229       &                                 kp_ext, pdt, pds, pdb, pdu, pdv ) 
     1243   SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt,     ps,    & 
     1244      &                                 pb,  pu,  pv,    kp_ext, pdt,   & 
     1245      &                                 pds, pdb, pdu,  pdv ) 
    12301246      !!--------------------------------------------------------------------- 
    12311247      !!                ***  ROUTINE zdf_vertical_average  *** 
    … …  
    12981314         pu(ji,jj) = pu(ji,jj) / zthick(ji,jj) 
    12991315         pv(ji,jj) = pv(ji,jj) / zthick(ji,jj) 
    1300          zthermal  = rab_n(ji,jj,1,jp_tem)   ! ideally use ibld not 1?? 
     1316         zthermal  = rab_n(ji,jj,1,jp_tem)   ! ideally use nbld not 1?? 
    13011317         zbeta     = rab_n(ji,jj,1,jp_sal) 
    13021318         pb(ji,jj) = grav * zthermal * pt(ji,jj) - grav * zbeta * ps(ji,jj) 
    … …  
    13151331               pdv(ji,jj) = pv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) /              & 
    13161332                  &                        MAX(1.0_wp , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 
    1317                zthermal   = rab_n(ji,jj,1,jp_tem)   ! ideally use ibld not 1?? 
     1333               zthermal   = rab_n(ji,jj,1,jp_tem)   ! ideally use nbld not 1?? 
    13181334               zbeta      = rab_n(ji,jj,1,jp_sal) 
    13191335               pdb(ji,jj) = grav * zthermal * pdt(ji,jj) - grav * zbeta * pds(ji,jj) 
    … …  
    14151431   END SUBROUTINE zdf_osm_velocity_rotation_3d 
    14161432 
    1417    SUBROUTINE zdf_osm_osbl_state( Kmm, ldconv, ldshear, k_ddh, pwb_ent, pwb_min, pshear,   & 
    1418       &                           phbl, phml, pdh, pdb_bl, pdu_bl, pdv_bl, pdb_ml, pdu_ml, pdv_ml ) 
     1433   SUBROUTINE zdf_osm_osbl_state( Kmm,    pwb_ent, pwb_min, pshear, phbl,     & 
     1434      &                           phml,   pdh,     pdb_bl,  pdu_bl, pdv_bl,   & 
     1435      &                           pdb_ml, pdu_ml,  pdv_ml ) 
    14191436      !!--------------------------------------------------------------------- 
    14201437      !!                 ***  ROUTINE zdf_osm_osbl_state  *** 
    … …  
    14291446      !!---------------------------------------------------------------------- 
    14301447      INTEGER,                  INTENT(in   ) ::   Kmm       ! Ocean time-level index 
    1431       LOGICAL,  DIMENSION(:,:), INTENT(  out) ::   ldconv    ! BL stability flags 
    1432       LOGICAL,  DIMENSION(:,:), INTENT(  out) ::   ldshear   ! Shear layers 
    1433       INTEGER,  DIMENSION(:,:), INTENT(  out) ::   k_ddh     ! k_ddh=0: active shear layer 
    1434       !                                                      ! k_ddh=1: shear layer not active 
    1435       !                                                      ! k_ddh=2: shear production low 
    14361448      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base 
    14371449      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_min   !    of well-mixed layer 
    … …  
    14671479      ! 
    14681480      ! Initialise INTENT(  out) arrays 
    1469       ldconv(:,:)  = .FALSE. 
    1470       ldshear(:,:) = .FALSE. 
    1471       k_ddh(:,:)   = 1 
     1481      l_conv(:,:)  = .FALSE. 
     1482      l_shear(:,:) = .FALSE. 
     1483      n_ddh(:,:)   = 1 
    14721484      pwb_ent(:,:) = zlarge 
    14731485      pwb_min(:,:) = zlarge 
    14741486      pshear(:,:)  = zlarge 
    14751487      ! 
    1476       ! Determins stability and set flag ldconv 
     1488      ! Determins stability and set flag l_conv 
    14771489      DO_2D( 0, 0, 0, 0 ) 
    14781490         IF ( shol(ji,jj) < 0._wp ) THEN 
    1479             ldconv(ji,jj) = .TRUE. 
     1491            l_conv(ji,jj) = .TRUE. 
    14801492         ELSE 
    1481             ldconv(ji,jj) = .FALSE. 
     1493            l_conv(ji,jj) = .FALSE. 
    14821494         ENDIF 
    14831495      END_2D 
    … …  
    14961508      ! 
    14971509      DO_2D( 0, 0, 0, 0 ) 
    1498          IF ( ldconv(ji,jj) ) THEN 
     1510         IF ( l_conv(ji,jj) ) THEN 
    14991511            IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
    15001512               zri_p(ji,jj) = MAX (  SQRT( pdb_bl(ji,jj) * pdh(ji,jj) / MAX( pdu_bl(ji,jj)**2 + pdv_bl(ji,jj)**2, 1e-8_wp ) ) *   & 
    … …  
    15281540#endif 
    15291541               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    1530                ! Test ensures k_ddh=0 is not selected. Change to zri_p<27 when  ! 
     1542               ! Test ensures n_ddh=0 is not selected. Change to zri_p<27 when  ! 
    15311543               ! full code available                                          ! 
    15321544               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    … …  
    15341546                  IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN 
    15351547                     ! Growing shear layer 
    1536                      k_ddh(ji,jj) = 0 
    1537                      ldshear(ji,jj) = .TRUE. 
     1548                     n_ddh(ji,jj) = 0 
     1549                     l_shear(ji,jj) = .TRUE. 
    15381550                  ELSE 
    1539                      k_ddh(ji,jj) = 1 
     1551                     n_ddh(ji,jj) = 1 
    15401552                     !             IF ( zri_b <= 1.5 .and. pshear(ji,jj) > 0._wp ) THEN 
    15411553                     ! Shear production large enough to determine layer charcteristics, but can't maintain a shear layer 
    1542                      ldshear(ji,jj) = .TRUE. 
     1554                     l_shear(ji,jj) = .TRUE. 
    15431555                     !             ELSE 
    15441556                  END IF 
    15451557               ELSE 
    1546                   k_ddh(ji,jj) = 2 
    1547                   ldshear(ji,jj) = .FALSE. 
     1558                  n_ddh(ji,jj) = 2 
     1559                  l_shear(ji,jj) = .FALSE. 
    15481560               END IF 
    15491561               ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline 
    15501562               !               pshear(ji,jj) = 0.5 * pshear(ji,jj) 
    1551                !               ldshear(ji,jj) = .FALSE. 
     1563               !               l_shear(ji,jj) = .FALSE. 
    15521564               !            ENDIF 
    15531565            ELSE   ! pdb_bl test, note pshear set to zero 
    1554                k_ddh(ji,jj) = 2 
    1555                ldshear(ji,jj) = .FALSE. 
     1566               n_ddh(ji,jj) = 2 
     1567               l_shear(ji,jj) = .FALSE. 
    15561568            ENDIF 
    15571569         ENDIF 
    … …  
    15601572      ! Calculate entrainment buoyancy flux due to surface fluxes. 
    15611573      DO_2D( 0, 0, 0, 0 ) 
    1562          IF ( ldconv(ji,jj) ) THEN 
     1574         IF ( l_conv(ji,jj) ) THEN 
    15631575            zwcor        = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln 
    15641576            zrf_conv     = TANH( ( swstrc(ji,jj) / zwcor )**0.69_wp ) 
    … …  
    15871599      ! 
    15881600      DO_2D( 0, 0, 0, 0 ) 
    1589          IF ( ldshear(ji,jj) ) THEN 
    1590             IF ( ldconv(ji,jj) ) THEN 
     1601         IF ( l_shear(ji,jj) ) THEN 
     1602            IF ( l_conv(ji,jj) ) THEN 
    15911603               ! Unstable OSBL 
    15921604               zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * pshear(ji,jj) 
    … …  
    15971609               END IF 
    15981610#endif 
    1599                IF ( k_ddh(ji,jj) == 0 ) THEN 
     1611               IF ( n_ddh(ji,jj) == 0 ) THEN 
    16001612                  ! Developing shear layer, additional shear production possible. 
    16011613 
    … …  
    16081620#ifdef key_osm_debug 
    16091621                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1610                      WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state k_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, & 
     1622                     WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, & 
    16111623                        & '  zshear=',pshear(ji,jj),'  zshear_u=', pshear_u 
    16121624                     FLUSH(narea+100) 
    … …  
    16161628               pwb_ent(ji,jj) = pwb_ent(ji,jj) + zwb_shr 
    16171629               !           pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * zwb0(ji,jj) 
    1618             ELSE   ! IF ( ldconv ) THEN - ENDIF 
     1630            ELSE   ! IF ( l_conv ) THEN - ENDIF 
    16191631               ! Stable OSBL  - shear production not coded for first attempt. 
    1620             ENDIF   ! ldconv 
    1621          END IF   ! ldshear 
    1622          IF ( ldconv(ji,jj) ) THEN 
     1632            ENDIF   ! l_conv 
     1633         END IF   ! l_shear 
     1634         IF ( l_conv(ji,jj) ) THEN 
    16231635            ! Unstable OSBL 
    16241636            pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * 2.0_wp * swbav(ji,jj) 
    1625          END IF  ! ldconv 
     1637         END IF  ! l_conv 
    16261638#ifdef key_osm_debug 
    16271639         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    … …  
    16431655      !! ** Purpose : Calculates the gradients below the OSBL 
    16441656      !! 
    1645       !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1657      !! ** Method  : Uses nbld and ibld_ext to determine levels to calculate the gradient. 
    16461658      !! 
    16471659      !!----------------------------------------------------------------------    
    … …  
    16641676      DO_2D( 0, 0, 0, 0 ) 
    16651677         IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
    1666             zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use ibld not 1?? 
     1678            zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use nbld not 1?? 
    16671679            zbeta    = rab_n(ji,jj,1,jp_sal) 
    16681680            jkb = kbase(ji,jj) 
    … …  
    16821694   END SUBROUTINE zdf_osm_external_gradients 
    16831695 
    1684    SUBROUTINE zdf_osm_calculate_dhdt( ldconv, ldshear, k_ddh, pdhdt, phbl, pdh, pwb_ent, pwb_min,   & 
    1685       &                               pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk, pvel_mle ) 
     1696   SUBROUTINE zdf_osm_calculate_dhdt( pdhdt,  phbl,         pdh,    pwb_ent,  pwb_min,   & 
     1697      &                               pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk,    & 
     1698      &                               pvel_mle ) 
    16861699      !!--------------------------------------------------------------------- 
    16871700      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
    … …  
    16921705      !! 
    16931706      !!---------------------------------------------------------------------- 
    1694       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
    1695       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldshear        ! Shear layers 
    1696       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! k_ddh=0: active shear layer 
    1697       !                                                              ! k_ddh=1: shear layer not active 
    1698       !                                                              ! k_ddh=2: shear production low 
    16991707      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl 
    17001708      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    … …  
    17281736      DO_2D( 0, 0, 0, 0 ) 
    17291737         ! 
    1730          IF ( ldshear(ji,jj) ) THEN 
     1738         IF ( l_shear(ji,jj) ) THEN 
    17311739            ! 
    1732             IF ( ldconv(ji,jj) ) THEN   ! Convective 
     1740            IF ( l_conv(ji,jj) ) THEN   ! Convective 
    17331741               ! 
    17341742               IF ( ln_osm_mle ) THEN 
    … …  
    17591767                        END IF 
    17601768#endif 
    1761                         IF ( k_ddh(ji,jj) == 1 ) THEN 
     1769                        IF ( n_ddh(ji,jj) == 1 ) THEN 
    17621770                           IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN 
    17631771                              zari = MIN( 1.5_wp * pdb_bl(ji,jj) /                                                   & 
    … …  
    17801788                           END IF 
    17811789#endif 
    1782                         ELSE IF ( k_ddh(ji,jj) == 0 ) THEN   ! Growing shear layer 
     1790                        ELSE IF ( n_ddh(ji,jj) == 0 ) THEN   ! Growing shear layer 
    17831791                           zddhdt = -1.0_wp * a_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
    17841792                              &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
    … …  
    17861794                        ELSE 
    17871795                           zddhdt = 0.0_wp 
    1788                         ENDIF   ! k_ddh 
     1796                        ENDIF   ! n_ddh 
    17891797                        pdhdt(ji,jj) = pdhdt(ji,jj) + zalpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
    17901798                           &                            pdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) /   & 
    … …  
    18051813               ENDIF   ! ln_osm_mle 
    18061814               ! 
    1807             ELSE   ! ldconv - Stable 
     1815            ELSE   ! l_conv - Stable 
    18081816               ! 
    18091817               pdhdt(ji,jj) = ( 0.06_wp + 0.52_wp * shol(ji,jj) / 2.0_wp ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 
    … …  
    18161824               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
    18171825               ! 
    1818             ENDIF   ! ldconv 
     1826            ENDIF   ! l_conv 
    18191827            ! 
    1820          ELSE   ! ldshear 
     1828         ELSE   ! l_shear 
    18211829            ! 
    1822             IF ( ldconv(ji,jj) ) THEN   ! Convective 
     1830            IF ( l_conv(ji,jj) ) THEN   ! Convective 
    18231831               ! 
    18241832               IF ( ln_osm_mle ) THEN 
    … …  
    18601868               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
    18611869               ! 
    1862             ENDIF  ! ldconv 
     1870            ENDIF  ! l_conv 
    18631871            ! 
    1864          ENDIF ! ldshear 
     1872         ENDIF ! l_shear 
    18651873#ifdef key_osm_debug 
    18661874         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    … …  
    18811889   END SUBROUTINE zdf_osm_calculate_dhdt 
    18821890 
    1883    SUBROUTINE zdf_osm_timestep_hbl( Kmm, kbld, kmld, ldconv, ldpyc, pdhdt, phbl, phbl_t, pt_bl, ps_bl, pwb_ent, pwb_fk_b ) 
     1891   SUBROUTINE zdf_osm_timestep_hbl( Kmm,   pdhdt,   phbl, phbl_t, pt_bl,   & 
     1892      &                             ps_bl, pwb_ent, pwb_fk_b ) 
    18841893      !!--------------------------------------------------------------------- 
    18851894      !!                ***  ROUTINE zdf_osm_timestep_hbl  *** 
    … …  
    18941903      !!---------------------------------------------------------------------- 
    18951904      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index 
    1896       INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kbld       ! BL base layer 
    1897       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld       ! ML base layer 
    1898       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv     ! BL stability flags 
    1899       LOGICAL,  DIMENSION(:,:),    INTENT(inout) ::   ldpyc      ! Pycnocline flags 
    19001905      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl 
    19011906      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phbl       ! BL depth 
    … …  
    19161921#ifdef key_osm_debug 
    19171922         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1918             WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',kmld(ji,jj),' trial ibld=', kbld(ji,jj) 
     1923            WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',nmld(ji,jj),' trial ibld=', nbld(ji,jj) 
    19191924            FLUSH(narea+100) 
    19201925         END IF 
    19211926#endif 
    1922          IF ( kbld(ji,jj) - kmld(ji,jj) > 1 ) THEN 
     1927         IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN 
    19231928            ! 
    19241929            ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    19251930            ! 
    19261931            zhbl_s   = hbl(ji,jj) 
    1927             jm       = kmld(ji,jj) 
     1932            jm       = nmld(ji,jj) 
    19281933            zthermal = rab_n(ji,jj,1,jp_tem) 
    19291934            zbeta    = rab_n(ji,jj,1,jp_sal) 
    19301935            ! 
    1931             IF ( ldconv(ji,jj) ) THEN   ! Unstable 
     1936            IF ( l_conv(ji,jj) ) THEN   ! Unstable 
    19321937               ! 
    19331938               IF( ln_osm_mle ) THEN 
    … …  
    19441949               END IF 
    19451950#endif 
    1946                DO jk = kmld(ji,jj), kbld(ji,jj) 
     1951               DO jk = nmld(ji,jj), nbld(ji,jj) 
    19471952                  zdb = MAX( grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -   & 
    19481953                     &                zbeta    * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + zvel_max 
    … …  
    19501955                  IF ( ln_osm_mle ) THEN 
    19511956                     zhbl_s = zhbl_s + MIN( rn_Dt * ( ( -1.0_wp * pwb_ent(ji,jj) - 2.0_wp * pwb_fk_b(ji,jj) ) / zdb ) /   & 
    1952                         &                   REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
     1957                        &                   REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
    19531958                  ELSE 
    19541959                     zhbl_s = zhbl_s + MIN( rn_Dt * ( -1.0_wp * pwb_ent(ji,jj) / zdb ) /   & 
    1955                         &                   REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
     1960                        &                   REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
    19561961                  ENDIF 
    19571962                  !                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
    19581963                  IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN 
    19591964                     zhbl_s = MIN( zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol ) 
    1960                      ldpyc(ji,jj) = .FALSE. 
     1965                     l_pyc(ji,jj) = .FALSE. 
    19611966                  ENDIF 
    19621967                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
    … …  
    19641969                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    19651970                     WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
    1966                      WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',ldpyc(ji,jj) 
     1971                     WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',l_pyc(ji,jj) 
    19671972                     FLUSH(narea+100) 
    19681973                  END IF 
    … …  
    19701975               END DO 
    19711976               hbl(ji,jj)  = zhbl_s 
    1972                kbld(ji,jj) = jm 
     1977               nbld(ji,jj) = jm 
    19731978            ELSE   ! Stable 
    19741979#ifdef key_osm_debug 
    … …  
    19781983               END IF 
    19791984#endif 
    1980                DO jk = kmld(ji,jj), kbld(ji,jj) 
     1985               DO jk = nmld(ji,jj), nbld(ji,jj) 
    19811986                  zdb = MAX(  grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -               & 
    19821987                     &                 zbeta    * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) +   & 
    … …  
    19901995                     &           ( 0.725_wp + 0.225_wp * EXP( -7.5_wp * shol(ji,jj) ) ) 
    19911996                  pdhdt(ji,jj) = pdhdt(ji,jj) + swbav(ji,jj) 
    1992                   zhbl_s = zhbl_s + MIN( pdhdt(ji,jj) / zdb * rn_Dt / REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ),   & 
     1997                  zhbl_s = zhbl_s + MIN( pdhdt(ji,jj) / zdb * rn_Dt / REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ),   & 
    19931998                     &                   e3w(ji,jj,jm,Kmm) ) 
    19941999                   
    … …  
    19962001                  IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 
    19972002                     zhbl_s      = MIN( zhbl_s,  gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - depth_tol ) 
    1998                      ldpyc(ji,jj) = .FALSE. 
     2003                     l_pyc(ji,jj) = .FALSE. 
    19992004                  ENDIF 
    20002005                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
    … …  
    20022007                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    20032008                     WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
    2004                      WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' shol',shol(ji,jj),' zdhdt',pdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',ldpyc(ji,jj) 
     2009                     WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' shol',shol(ji,jj),' zdhdt',pdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',l_pyc(ji,jj) 
    20052010                     FLUSH(narea+100) 
    20062011                  END IF 
    20072012#endif 
    20082013               END DO 
    2009             ENDIF   ! IF ( ldconv ) 
     2014            ENDIF   ! IF ( l_conv ) 
    20102015            hbl(ji,jj)  = MAX( zhbl_s, gdepw(ji,jj,4,Kmm) ) 
    2011             kbld(ji,jj) = MAX( jm, 4 ) 
     2016            nbld(ji,jj) = MAX( jm, 4 ) 
    20122017         ELSE 
    20132018            ! change zero or one model level. 
    20142019            hbl(ji,jj) = MAX( phbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) 
    20152020         ENDIF 
    2016          phbl(ji,jj) = gdepw(ji,jj,kbld(ji,jj),Kmm) 
     2021         phbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) 
    20172022#ifdef key_osm_debug 
    20182023         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    2019             WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', phbl(ji,jj),' ibld=', kbld(ji,jj) 
     2024            WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', phbl(ji,jj),' ibld=', nbld(ji,jj) 
    20202025            FLUSH(narea+100) 
    20212026         END IF 
    … …  
    20272032   END SUBROUTINE zdf_osm_timestep_hbl 
    20282033 
    2029    SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, kbld, kmld, ldshear, ldconv, k_ddh, pdh, phml, pdhdt, pdb_bl, phbl, pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 
     2034   SUBROUTINE zdf_osm_pycnocline_thickness( Kmm,  pdh,     phml,         pdhdt, pdb_bl,   & 
     2035      &                                     phbl, pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 
    20302036      !!--------------------------------------------------------------------- 
    20312037      !!            ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
    … …  
    20412047      !!---------------------------------------------------------------------- 
    20422048      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    2043       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
    2044       INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kmld           ! ML base layer 
    2045       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldshear        ! Shear layers 
    2046       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
    2047       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! k_ddh=0: active shear layer 
    2048       !                                                              ! k_ddh=1: shear layer not active 
    2049       !                                                              ! k_ddh=2: shear production low 
    20502049      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   pdh            ! Pycnocline thickness 
    20512050      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phml           ! ML depth 
    … …  
    20702069      DO_2D( 0, 0, 0, 0 ) 
    20712070         ! 
    2072          IF ( ldshear(ji,jj) ) THEN 
     2071         IF ( l_shear(ji,jj) ) THEN 
    20732072            ! 
    2074             IF ( ldconv(ji,jj) ) THEN 
     2073            IF ( l_conv(ji,jj) ) THEN 
    20752074               ! 
    20762075               IF ( pdb_bl(ji,jj) > 1e-15_wp ) THEN 
    2077                   IF ( k_ddh(ji,jj) == 0 ) THEN 
     2076                  IF ( n_ddh(ji,jj) == 0 ) THEN 
    20782077                     zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
    20792078                     ! ddhdt for pycnocline determined in osm_calculate_dhdt 
    … …  
    21052104               END IF 
    21062105               ! 
    2107             ELSE   ! ldconv 
     2106            ELSE   ! l_conv 
    21082107               ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 
    21092108               ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln) 
    … …  
    21252124            ENDIF 
    21262125            ! 
    2127          ELSE   ! ldshear = .FALSE., calculate ddhdt here 
     2126         ELSE   ! l_shear = .FALSE., calculate ddhdt here 
    21282127            ! 
    2129             IF ( ldconv(ji,jj) ) THEN 
     2128            IF ( l_conv(ji,jj) ) THEN 
    21302129               ! 
    21312130               IF( ln_osm_mle ) THEN 
    … …  
    21692168               IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 
    21702169               ! Alan: this hml is never defined or used 
    2171             ELSE   ! IF (ldconv) 
     2170            ELSE   ! IF (l_conv) 
    21722171               ! 
    21732172               ztau = hbl(ji,jj) / MAX( svstr(ji,jj), epsln ) 
    … …  
    21872186               IF ( pdhdt(ji,jj) < 0.0_wp .AND. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref   ! Can be a problem with dh>hbl for 
    21882187               !                                                                                !    rapid collapse 
    2189             END IF   ! IF (ldconv) 
     2188            END IF   ! IF (l_conv) 
    21902189            ! 
    2191          END IF   ! ldshear 
     2190         END IF   ! l_shear 
    21922191         ! 
    21932192         hml(ji,jj)  = hbl(ji,jj) - dh(ji,jj) 
    2194          inhml       = MAX( INT( dh(ji,jj) / MAX( e3t(ji,jj,kbld(ji,jj)-1,Kmm), 1e-3_wp ) ), 1 ) 
    2195          kmld(ji,jj) = MAX( kbld(ji,jj) - inhml, 3 ) 
    2196          phml(ji,jj) = gdepw(ji,jj,kmld(ji,jj),Kmm) 
     2193         inhml       = MAX( INT( dh(ji,jj) / MAX( e3t(ji,jj,nbld(ji,jj)-1,Kmm), 1e-3_wp ) ), 1 ) 
     2194         nmld(ji,jj) = MAX( nbld(ji,jj) - inhml, 3 ) 
     2195         phml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm) 
    21972196         pdh(ji,jj)  = phbl(ji,jj) - phml(ji,jj) 
    21982197#ifdef key_osm_debug 
    21992198         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    22002199            WRITE(narea+100,'(4(a,g11.3),2(a,i7),/,5(a,g11.3),/)') 'end of zdf_osm_pycnocline_thickness:hml=',hml(ji,jj), & 
    2201                & '  zhml=',phml(ji,jj),' zdh=', pdh(ji,jj), '  dh=', dh(ji,jj), ' imld=', kmld(ji,jj), ' inhml=', inhml, & 
     2200               & '  zhml=',phml(ji,jj),' zdh=', pdh(ji,jj), '  dh=', dh(ji,jj), ' imld=', nmld(ji,jj), ' inhml=', inhml, & 
    22022201               & 'zvel_max=', zvel_max, ' ztau=', ztau,' zdh_ref=', zdh_ref, ' zar=', zari, ' zddhdt=', zddhdt 
    22032202            FLUSH(narea+100) 
    … …  
    22112210   END SUBROUTINE zdf_osm_pycnocline_thickness 
    22122211 
    2213    SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kbld, kp_ext, ldconv, ldpyc, pdbdz, palpha, pdh, pdb_bl, phbl, pdbdz_bl_ext, phml, pdb_ml, pdhdt ) 
     2212   SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm,    kp_ext, pdbdz,        palpha, pdh,      & 
     2213      &                                             pdb_bl, phbl,   pdbdz_bl_ext, phml,   pdb_ml,   & 
     2214      &                                             pdhdt ) 
    22142215      !!--------------------------------------------------------------------- 
    22152216      !!       ***  ROUTINE zdf_osm_pycnocline_buoyancy_profiles  *** 
    … …  
    22212222      !!---------------------------------------------------------------------- 
    22222223      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    2223       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
    22242224      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kp_ext         ! External-level offsets 
    2225       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
    2226       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
    22272225      REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pdbdz          ! Gradients in the pycnocline 
    22282226      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   palpha 
    … …  
    22482246      DO_2D( 0, 0, 0, 0 ) 
    22492247         ! 
    2250          IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     2248         IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    22512249            ! 
    2252             IF ( ldconv(ji,jj) ) THEN   ! Convective conditions 
     2250            IF ( l_conv(ji,jj) ) THEN   ! Convective conditions 
    22532251               ! 
    2254                IF ( ldpyc(ji,jj) ) THEN 
     2252               IF ( l_pyc(ji,jj) ) THEN 
    22552253                  ! 
    22562254                  zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 
    … …  
    22682266                  zbgrad = palpha(ji,jj) * pdb_ml(ji,jj) * ztmp + pdbdz_bl_ext(ji,jj) 
    22692267                  zgamma_b_nd = pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / MAX( pdb_ml(ji,jj), epsln ) 
    2270                   DO jk = 2, kbld(ji,jj) 
     2268                  DO jk = 2, nbld(ji,jj) 
    22712269                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) * ztmp 
    22722270                     IF ( znd <= zzeta_m ) THEN 
    … …  
    23002298                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 
    23012299                        zbgrad = pdb_bl(ji,jj) * ztmp 
    2302                         DO jk = 2, kbld(ji,jj) 
     2300                        DO jk = 2, nbld(ji,jj) 
    23032301                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
    23042302                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
    … …  
    23072305                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 
    23082306                        zbgrad = pdb_bl(ji,jj) * ztmp 
    2309                         DO jk = 2, kbld(ji,jj) 
     2307                        DO jk = 2, nbld(ji,jj) 
    23102308                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
    23112309                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
    … …  
    23232321               END IF         ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    23242322               ! 
    2325             END IF            ! IF (ldconv) 
     2323            END IF            ! IF (l_conv) 
    23262324            ! 
    2327          END IF   ! IF ( kbld(ji,jj) < mbkt(ji,jj) ) 
     2325         END IF   ! IF ( nbld(ji,jj) < mbkt(ji,jj) ) 
    23282326         ! 
    23292327      END_2D 
    … …  
    23372335   END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 
    23382336 
    2339    SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, kbld, kmld, ldconv, ldshear, ldpyc, ldcoup, k_ddh, pdiffut, pviscos,       & 
    2340       &                                      phbl, phml, pdh, pdhdt, pshear, pb_bl, pu_bl, pv_bl, pb_ml, pu_ml, pv_ml, pwb_ent,   & 
    2341       &                                      pwb_min ) 
     2337   SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb,     Kmm,   pdiffut, pviscos, phbl,    & 
     2338      &                                      phml,    pdh,   pdhdt,   pshear,  pb_bl,   & 
     2339      &                                      pu_bl,   pv_bl, pb_ml,   pu_ml,   pv_ml,   & 
     2340      &                                      pwb_ent, pwb_min ) 
    23422341      !!--------------------------------------------------------------------- 
    23432342      !!           ***  ROUTINE zdf_osm_diffusivity_viscosity  *** 
    … …  
    23502349      !!---------------------------------------------------------------------- 
    23512350      INTEGER,                     INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
    2352       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
    2353       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld           ! ML base layer 
    2354       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
    2355       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldshear        ! Shear layers 
    2356       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
    2357       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldcoup         ! Coupling to bottom 
    2358       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! Type of shear layer 
    23592351      REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
    23602352      REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pviscos        ! Viscosity 
    … …  
    23982390      ! 
    23992391      DO_2D( 0, 0, 0, 0 ) 
    2400          IF ( ldconv(ji,jj) ) THEN 
     2392         IF ( l_conv(ji,jj) ) THEN 
    24012393            ! 
    24022394            zvel_sc_pyc = ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25_wp * pshear(ji,jj) * phbl(ji,jj) )**pthird 
    … …  
    24152407#endif 
    24162408            ! 
    2417             IF ( ldpyc(ji,jj) ) THEN 
     2409            IF ( l_pyc(ji,jj) ) THEN 
    24182410               zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj) 
    24192411               zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 *                           & 
    … …  
    24282420#endif 
    24292421               ! 
    2430                IF ( ldshear(ji,jj) .AND. k_ddh(ji,jj) /= 2 ) THEN 
     2422               IF ( l_shear(ji,jj) .AND. n_ddh(ji,jj) /= 2 ) THEN 
    24312423                  ztmp = rn_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj) 
    24322424                  zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + ztmp 
    … …  
    24712463               zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj)   ! ag 19/03 
    24722464               zvispyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03 
    2473                IF(ldcoup(ji,jj) ) THEN   ! ag 19/03 
     2465               IF(l_coup(ji,jj) ) THEN   ! ag 19/03 
    24742466                  ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub| 
    24752467                  !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90 
    … …  
    25352527      ! 
    25362528      DO_2D( 0, 0, 0, 0 ) 
    2537       IF ( ldconv(ji,jj) ) THEN 
    2538          DO jk = 2, kmld(ji,jj)   ! Mixed layer diffusivity 
     2529      IF ( l_conv(ji,jj) ) THEN 
     2530         DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity 
    25392531            zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    25402532            pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
    … …  
    25452537         ! Coupling to bottom 
    25462538         ! 
    2547          IF ( ldcoup(ji,jj) ) THEN                                                         ! ag 19/03 
    2548             DO jk = mbkt(ji,jj), kmld(ji,jj), -1                                           ! ag 19/03 
     2539         IF ( l_coup(ji,jj) ) THEN                                                         ! ag 19/03 
     2540            DO jk = mbkt(ji,jj), nmld(ji,jj), -1                                           ! ag 19/03 
    25492541               zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )   ! ag 19/03 
    25502542               pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03 
    … …  
    25532545         ENDIF                                                                             ! ag 19/03 
    25542546         ! Pycnocline 
    2555          IF ( ldpyc(ji,jj) ) THEN  
     2547         IF ( l_pyc(ji,jj) ) THEN  
    25562548            ! Diffusivity and viscosity profiles in the pycnocline given by 
    2557             ! cubic polynomial. Note, if ldpyc TRUE can't be coupled to seabed. 
     2549            ! cubic polynomial. Note, if l_pyc TRUE can't be coupled to seabed. 
    25582550            za_cubic = 0.5_wp 
    25592551            zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
    … …  
    25682560            zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic ) 
    25692561            zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic 
    2570             DO jk = kmld(ji,jj) , kbld(ji,jj) 
     2562            DO jk = nmld(ji,jj) , nbld(ji,jj) 
    25712563               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp ) 
    25722564               ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ) 
    … …  
    25812573            END DO 
    25822574!                  IF ( pdhdt(ji,jj) > 0._wp ) THEN 
    2583 !                     zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 
    2584 !                     zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 
     2575!                     zdiffut(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
     2576!                     zviscos(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
    25852577!                  ELSE 
    2586 !                     zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
    2587 !                     zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     2578!                     zdiffut(ji,jj,nbld(ji,jj)) = 0._wp 
     2579!                     zviscos(ji,jj,nbld(ji,jj)) = 0._wp 
    25882580!                  ENDIF 
    25892581         ENDIF 
    25902582      ELSE 
    25912583         ! Stable conditions 
    2592          DO jk = 2, kbld(ji,jj) 
     2584         DO jk = 2, nbld(ji,jj) 
    25932585            zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
    25942586            pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp 
    … …  
    25972589         ! 
    25982590         IF ( pdhdt(ji,jj) > 0.0_wp ) THEN 
    2599             pdiffut(ji,jj,kbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, kbld(ji,jj), Kmm) 
    2600             pviscos(ji,jj,kbld(ji,jj)) = pdiffut(ji,jj,kbld(ji,jj)) 
     2591            pdiffut(ji,jj,nbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, nbld(ji,jj), Kmm) 
     2592            pviscos(ji,jj,nbld(ji,jj)) = pdiffut(ji,jj,nbld(ji,jj)) 
    26012593         ENDIF 
    2602       ENDIF   ! End if ( ldconv ) 
     2594      ENDIF   ! End if ( l_conv ) 
    26032595      ! 
    26042596      END_2D 
    … …  
    26082600   END SUBROUTINE zdf_osm_diffusivity_viscosity 
    26092601 
    2610    SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, pshear,             & 
    2611       &                          pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml,                  & 
    2612       &                          pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) 
     2602   SUBROUTINE zdf_osm_fgr_terms( Kmm,        kp_ext,  phbl,         phml,         pdh,         & 
     2603      &                          pdhdt,      pshear,  pdt_bl,       pds_bl,       pdb_bl,      & 
     2604      &                          pdu_bl,     pdv_bl,  pdt_ml,       pds_ml,       pdb_ml,      & 
     2605      &                          pdu_ml,     pdv_ml,  pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc,   & 
     2606      &                          palpha_pyc, pdiffut, pviscos ) 
    26132607      !!--------------------------------------------------------------------- 
    26142608      !!                 ***  ROUTINE zdf_osm_fgr_terms *** 
    … …  
    26202614      !!---------------------------------------------------------------------- 
    26212615      INTEGER,                     INTENT(in   ) ::   Kmm            ! Time-level index 
    2622       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
    2623       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld           ! ML base layer 
    26242616      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kp_ext         ! Offset for external level 
    2625       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
    2626       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
    2627       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! Type of shear layer 
    26282617      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    26292618      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
    … …  
    26872676      jkm_mld = 0 
    26882677      DO_2D( 0, 0, 0, 0 ) 
    2689          IF ( kbld(ji,jj) > jkm_bld ) jkm_bld = kbld(ji,jj) 
    2690          IF ( kmld(ji,jj) < jkf_mld ) jkf_mld = kmld(ji,jj) 
    2691          IF ( kmld(ji,jj) > jkm_mld ) jkm_mld = kmld(ji,jj) 
     2678         IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj) 
     2679         IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj) 
     2680         IF ( nmld(ji,jj) > jkm_mld ) jkm_mld = nmld(ji,jj) 
    26922681      END_2D 
    26932682      ! 
    26942683      ! Stokes term in scalar flux, flux-gradient relationship 
    26952684      ! ------------------------------------------------------ 
    2696       WHERE ( ldconv(A2D(0)) ) 
     2685      WHERE ( l_conv(A2D(0)) ) 
    26972686         zsc_wth_1(:,:) = swstrl(A2D(0))**3 * swth0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    26982687         zsc_ws_1(:,:)  = swstrl(A2D(0))**3 * sws0(A2D(0))  / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    … …  
    27022691      ENDWHERE 
    27032692      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
    2704          IF ( ldconv(ji,jj) ) THEN 
    2705             IF ( jk <= kmld(ji,jj) ) THEN 
     2693         IF ( l_conv(ji,jj) ) THEN 
     2694            IF ( jk <= nmld(ji,jj) ) THEN 
    27062695               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    27072696               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   & 
    … …  
    27112700            END IF 
    27122701         ELSE   ! Stable conditions 
    2713             IF ( jk <= kbld(ji,jj) ) THEN 
     2702            IF ( jk <= nbld(ji,jj) ) THEN 
    27142703               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    27152704               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) *   & 
    … …  
    27182707                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 
    27192708            END IF 
    2720          END IF   ! Check on ldconv 
     2709         END IF   ! Check on l_conv 
    27212710      END_3D 
    27222711      ! 
    … …  
    27292718      ! svstr since term needs to go to zero as swstrl goes to zero) 
    27302719      ! --------------------------------------------------------------------- 
    2731       WHERE ( ldconv(A2D(0)) ) 
     2720      WHERE ( l_conv(A2D(0)) ) 
    27322721         zsc_uw_1(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) /   & 
    27332722            &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 
    … …  
    27422731      ENDWHERE 
    27432732      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
    2744          IF ( ldconv(ji,jj) ) THEN 
    2745             IF ( jk <= kmld(ji,jj) ) THEN 
     2733         IF ( l_conv(ji,jj) ) THEN 
     2734            IF ( jk <= nmld(ji,jj) ) THEN 
    27462735               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    27472736               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp   * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) +     & 
    … …  
    27522741            END IF 
    27532742         ELSE   ! Stable conditions 
    2754             IF ( jk <= kbld(ji,jj) ) THEN   ! Corrected to ibld 
     2743            IF ( jk <= nbld(ji,jj) ) THEN   ! Corrected to nbld 
    27552744               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    27562745               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) *             & 
    … …  
    27622751      IF(narea==nn_narea_db) THEN 
    27632752         ji=iloc_db; jj=jloc_db 
    2764          jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) 
     2753         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    27652754         WRITE(narea+100,'(a,g11.3)')'Stokes contrib to ghamt/s:  zsc_wth_1=',zsc_wth_1(ji,jj), '  zsc_ws_1=',zsc_ws_1(ji,jj) 
    27662755         WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
    27672756         WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) 
    2768          IF( ldconv(ji,jj) ) THEN 
     2757         IF( l_conv(ji,jj) ) THEN 
    27692758            WRITE(narea+100,'(3(a,g11.3))')'Stokes contrib to ghamu/v:  zsc_uw_1=',zsc_uw_1(ji,jj), '  zsc_vw_1=',zsc_vw_1(ji,jj), & 
    27702759               &' zsc_uw_2=',zsc_uw_2(ji,jj) 
    … …  
    27822771      ! (X0.3) and pressure (X0.5)] 
    27832772      ! ---------------------------------------------------------------------- 
    2784       WHERE ( ldconv(A2D(0)) ) 
     2773      WHERE ( l_conv(A2D(0)) ) 
    27852774         zsc_wth_1(:,:) = swbav(A2D(0)) * swth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) /   & 
    27862775            &             ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    … …  
    27922781      ENDWHERE 
    27932782      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
    2794          IF ( ldconv(ji,jj) ) THEN 
    2795             IF ( jk <= kmld(ji,jj) ) THEN 
     2783         IF ( l_conv(ji,jj) ) THEN 
     2784            IF ( jk <= nmld(ji,jj) ) THEN 
    27962785               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    27972786               ! Calculate turbulent time scale 
    … …  
    28062795            END IF 
    28072796         ELSE   ! Stable conditions 
    2808             IF ( jk <= kbld(ji,jj) ) THEN 
     2797            IF ( jk <= nbld(ji,jj) ) THEN 
    28092798               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
    28102799               ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj) 
    … …  
    28132802      END_3D 
    28142803      DO_2D( 0, 0, 0, 0 ) 
    2815          IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 
     2804         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 
    28162805            ztau_sc_u(ji,jj)    = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *                             & 
    28172806               &                ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.5_wp ) 
    … …  
    28412830      END_2D 
    28422831      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
    2843          IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2832         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    28442833            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    28452834            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) -                                                                                & 
    … …  
    28522841         END IF 
    28532842      END_3D 
    2854       jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) 
     2843      jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    28552844      IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    28562845         WRITE(narea+100,'(3(a,g11.3))')'lpyc= lconv=T: ztau_sc_u=',ztau_sc_u(ji,jj),' zwth_ent=',zwth_ent(ji,jj),   & 
    … …  
    28622851      END IF 
    28632852      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
    2864          IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2853         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    28652854            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    28662855#endif 
    2867             IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) .AND. kbld(ji,jj) - kmld(ji,jj) > 3 ) THEN 
     2856            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) .AND. nbld(ji,jj) - nmld(ji,jj) > 3 ) THEN 
    28682857               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              & 
    28692858                  &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
    … …  
    28822871      ! 
    28832872      zsc_vw_1(:,:) = 0.0_wp 
    2884       WHERE ( ldconv(A2D(0)) ) 
     2873      WHERE ( l_conv(A2D(0)) ) 
    28852874         zsc_uw_1(:,:) = -1.0_wp * swb0(A2D(0)) * sustar(A2D(0))**2 * phml(A2D(0)) /   & 
    28862875            &            ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    … …  
    28912880      ENDWHERE 
    28922881      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
    2893          IF ( ldconv(ji,jj) ) THEN 
    2894             IF ( jk <= kmld(ji,jj) ) THEN 
     2882         IF ( l_conv(ji,jj) ) THEN 
     2883            IF ( jk <= nmld(ji,jj) ) THEN 
    28952884               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    28962885               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp *   & 
    … …  
    29002889            END IF 
    29012890         ELSE   ! Stable conditions 
    2902             IF ( jk <= kbld(ji,jj) ) THEN 
     2891            IF ( jk <= nbld(ji,jj) ) THEN 
    29032892               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 
    29042893               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
    … …  
    29082897      ! 
    29092898      DO_2D( 0, 0, 0, 0 ) 
    2910          IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 
    2911             IF ( k_ddh(ji,jj) == 0 ) THEN 
     2899         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 
     2900            IF ( n_ddh(ji,jj) == 0 ) THEN 
    29122901               ! Place holding code. Parametrization needs checking for these conditions. 
    29132902               zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
    … …  
    29272916      END_2D 
    29282917      DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array. 
    2929          IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2918         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    29302919            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    29312920            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) *                 & 
    … …  
    29352924               &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   & 
    29362925               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
    2937          END IF   ! ldconv .AND. ldpyc 
     2926         END IF   ! l_conv .AND. l_pyc 
    29382927      END_3D 
    29392928      ! 
    … …  
    29412930      IF(narea==nn_narea_db) THEN 
    29422931         ji=iloc_db; jj=jloc_db 
    2943          jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) 
     2932         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    29442933         WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc contribs to ghamt/s:  zsc_wth_1=',zsc_wth_1(ji,jj), '  zsc_ws_1=',zsc_ws_1(ji,jj) 
    29452934         WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
    29462935         WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) 
    2947          IF( ldconv(ji,jj) ) THEN 
     2936         IF( l_conv(ji,jj) ) THEN 
    29482937            WRITE(narea+100,'(3(a,g11.3))')'Stokes + buoy + pyc contribs to ghamu/v:  zsc_uw_1=',zsc_uw_1(ji,jj), '  zsc_vw_1=',zsc_vw_1(ji,jj), & 
    29492938               &' zsc_uw_2=',zsc_uw_2(ji,jj) 
    … …  
    29662955      ! (X0.3) ] 
    29672956      ! ----------------------------------------------------------------------- 
    2968       WHERE ( ldconv(A2D(0)) ) 
     2957      WHERE ( l_conv(A2D(0)) ) 
    29692958         zsc_wth_1(:,:) = swth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 
    29702959         zsc_ws_1(:,:)  = sws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 
    2971          WHERE ( ldpyc(A2D(0)) )   ! Pycnocline scales 
     2960         WHERE ( l_pyc(A2D(0)) )   ! Pycnocline scales 
    29722961            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pdt_ml(A2D(0)) 
    29732962            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pds_ml(A2D(0)) 
    … …  
    29782967      END WHERE 
    29792968      DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) 
    2980          IF ( ldconv(ji,jj) ) THEN 
    2981             IF ( ( jk > 1 ) .AND. ( jk <= kmld(ji,jj) ) ) THEN 
     2969         IF ( l_conv(ji,jj) ) THEN 
     2970            IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN 
    29822971               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    29832972               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   & 
    … …  
    29882977            ! 
    29892978            ! may need to comment out lpyc block 
    2990             IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN   ! Pycnocline 
     2979            IF ( l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN   ! Pycnocline 
    29912980               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    29922981               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) 
    … …  
    29952984         ELSE 
    29962985            IF( pdhdt(ji,jj) > 0. ) THEN 
    2997                IF ( ( jk > 1 ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2986               IF ( ( jk > 1 ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    29982987                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    29992988                  znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
    … …  
    30072996      END_3D 
    30082997      ! 
    3009       WHERE ( ldconv(A2D(0)) ) 
     2998      WHERE ( l_conv(A2D(0)) ) 
    30102999         zsc_uw_1(:,:) = sustar(A2D(0))**2 
    30113000         zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phml(A2D(0)) 
    … …  
    30193008      ENDWHERE 
    30203009      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
    3021          IF ( ldconv(ji,jj) ) THEN 
    3022             IF ( jk <= kmld(ji,jj) ) THEN 
     3010         IF ( l_conv(ji,jj) ) THEN 
     3011            IF ( jk <= nmld(ji,jj) ) THEN 
    30233012               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    30243013               zznd_d  = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    … …  
    30313020            END IF 
    30323021         ELSE 
    3033             IF ( jk <= kbld(ji,jj) ) THEN 
     3022            IF ( jk <= nbld(ji,jj) ) THEN 
    30343023               znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
    30353024               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    … …  
    30513040      IF(narea==nn_narea_db) THEN 
    30523041         ji=iloc_db; jj=jloc_db 
    3053          jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) 
     3042         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    30543043         WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/s:  zsc_wth_1=',zsc_wth_1(ji,jj), '  zsc_ws_1=',zsc_ws_1(ji,jj) 
    3055          IF (ldpyc(ji,jj)) WRITE(narea+100,'(2(a,g11.3))') 'zsc_wth_pyc=', zsc_wth_pyc(ji,jj), '  zsc_wth_pyc=',zsc_wth_pyc(ji,jj) 
     3044         IF (l_pyc(ji,jj)) WRITE(narea+100,'(2(a,g11.3))') 'zsc_wth_pyc=', zsc_wth_pyc(ji,jj), '  zsc_wth_pyc=',zsc_wth_pyc(ji,jj) 
    30563045         WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
    30573046         WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) 
    3058          IF( ldconv(ji,jj) ) THEN 
     3047         IF( l_conv(ji,jj) ) THEN 
    30593048            WRITE(narea+100,'(2(a,g11.3))')'Unstable; transport contrib to ghamu/v:  zsc_uw_1=',zsc_uw_1(ji,jj), '  zsc_vw_1=',zsc_vw_1(ji,jj) 
    30603049         ELSE 
    … …  
    30833072      ! of the boundary layer. 
    30843073      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
    3085          IF ( ( .NOT. ldconv(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     3074         IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    30863075            znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about 
    30873076            IF ( znd >= 0.0_wp ) THEN 
    … …  
    31043093      END IF 
    31053094      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
    3106          IF ( ldconv (ji,jj) ) THEN 
     3095         IF ( l_conv (ji,jj) ) THEN 
    31073096            ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
    31083097            !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
    … …  
    31143103            !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
    31153104            !                       &      )/ (zdh(ji,jj)  + epsln ) 
    3116             !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     3105            !                  DO jk = 2, nbld(ji,jj) - 1 + ibld_ext 
    31173106            !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
    31183107            !                     IF ( znd <= 0.0 ) THEN 
    … …  
    31253114            !                  END DO 
    31263115         ELSE   ! Stable conditions 
    3127             IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     3116            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    31283117               ! Pycnocline profile only defined when depth steady of increasing. 
    31293118               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
    … …  
    31343123                        zsgrad = pds_bl(ji,jj) * ztmp 
    31353124                        zbgrad = pdb_bl(ji,jj) * ztmp 
    3136                         IF ( jk <= kbld(ji,jj) ) THEN 
     3125                        IF ( jk <= nbld(ji,jj) ) THEN 
    31373126                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
    31383127                           zdtdz_pyc =  ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
    … …  
    31503139                        zsgrad = pds_bl(ji,jj) * ztmp 
    31513140                        zbgrad = pdb_bl(ji,jj) * ztmp 
    3152                         IF ( jk <= kbld(ji,jj) ) THEN 
     3141                        IF ( jk <= nbld(ji,jj) ) THEN 
    31533142                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
    31543143                           zdtdz_pyc =  ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
    … …  
    31723161      END IF 
    31733162      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
    3174          IF ( .NOT. ldconv (ji,jj) ) THEN 
    3175             IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     3163         IF ( .NOT. l_conv (ji,jj) ) THEN 
     3164            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    31763165               zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj) 
    31773166               zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj) 
    3178                IF ( jk <= kbld(ji,jj) ) THEN 
     3167               IF ( jk <= nbld(ji,jj) ) THEN 
    31793168                  znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
    31803169                  IF ( znd < 1.0 ) THEN 
    … …  
    32073196      ! 
    32083197      DO_2D( 0, 0, 0, 0 ) 
    3209          ghamt(ji,jj,kbld(ji,jj)) = 0.0_wp 
    3210          ghams(ji,jj,kbld(ji,jj)) = 0.0_wp 
    3211          ghamu(ji,jj,kbld(ji,jj)) = 0.0_wp 
    3212          ghamv(ji,jj,kbld(ji,jj)) = 0.0_wp 
     3198         ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp 
     3199         ghams(ji,jj,nbld(ji,jj)) = 0.0_wp 
     3200         ghamu(ji,jj,nbld(ji,jj)) = 0.0_wp 
     3201         ghamv(ji,jj,nbld(ji,jj)) = 0.0_wp 
    32133202      END_2D 
    32143203#ifdef key_osm_debug 
    32153204      IF(narea==nn_narea_db) THEN 
    32163205         ji=iloc_db; jj=jloc_db 
    3217          jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) 
     3206         jl = nmld(ji,jj) - 1; jm = MIN(nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    32183207         WRITE(narea+100,'(a)')'Tweak gham[uv] to go to zero near surface, add pycnocline viscosity/diffusivity  & set=0 at ibld' 
    32193208         WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
    … …  
    32363225   END SUBROUTINE zdf_osm_fgr_terms 
    32373226 
    3238    SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, kbld, pmld, pdtdx, pdtdy, pdsdx, pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle ) 
     3227   SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm,   pmld,      pdtdx,     pdtdy, pdsdx,   & 
     3228      &                                          pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle ) 
    32393229      !!---------------------------------------------------------------------- 
    32403230      !!          ***  ROUTINE zdf_osm_zmld_horizontal_gradients  *** 
    … …  
    32503240      !!---------------------------------------------------------------------- 
    32513241      INTEGER,                  INTENT(in   ) ::   Kmm          ! Time-level index 
    3252       INTEGER,  DIMENSION(:,:), INTENT(in   ) ::   kbld         ! BL base layer 
    32533242      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == ! 
    32543243      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization 
    … …  
    32863275      END_3D 
    32873276      DO_2D( 1, 1, 1, 1 ) 
    3288          mld_prof(ji,jj) = MAX( mld_prof(ji,jj), kbld(ji,jj) )   ! Ensure mld_prof .ge. kbld 
     3277         mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) )   ! Ensure mld_prof .ge. nbld 
    32893278         pmld(ji,jj)     = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
    32903279      END_2D 
    … …  
    33373326   END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
    33383327 
    3339    SUBROUTINE zdf_osm_osbl_state_fk( Kmm, kbld, ldconv, ldpyc, ldflux, ldmle, pwb_fk, pt_mle, ps_mle, pb_mle, phbl,   & 
    3340       &                              phmle, pwb_ent, pt_bl, ps_bl, pb_bl, pdb_bl, pdbds_mle ) 
     3328   SUBROUTINE zdf_osm_osbl_state_fk( Kmm,  pwb_fk,  pt_mle,  ps_mle, pb_mle,   & 
     3329      &                              phbl, phmle,   pwb_ent, pt_bl,  ps_bl,    & 
     3330      &                              pb_bl, pdb_bl, pdbds_mle ) 
    33413331      !!--------------------------------------------------------------------- 
    33423332      !!               ***  ROUTINE zdf_osm_osbl_state_fk  *** 
    33433333      !! 
    33443334      !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is 
    3345       !!              returned in the logicals ldpyc, ldflux and ldmle. Used 
     3335      !!              returned in the logicals l_pyc, l_flux and ldmle. Used 
    33463336      !!              with Fox-Kemper scheme. 
    3347       !!                ldpyc  :: determines whether pycnocline flux-grad 
     3337      !!                l_pyc  :: determines whether pycnocline flux-grad 
    33483338      !!                          relationship needs to be determined 
    3349       !!                ldflux :: determines whether effects of surface flux 
     3339      !!                l_flux :: determines whether effects of surface flux 
    33503340      !!                          extend below the base of the OSBL 
    33513341      !!                ldmle  :: determines whether the layer with MLE is 
    … …  
    33583348      ! Outputs 
    33593349      INTEGER,                  INTENT(in   ) ::   Kmm         ! Time-level index 
    3360       INTEGER,  DIMENSION(:,:), INTENT(in   ) ::   kbld        ! BL base layer 
    3361       LOGICAL,  DIMENSION(:,:), INTENT(in   ) ::   ldconv      ! BL stability flags 
    3362       LOGICAL,  DIMENSION(:,:), INTENT(inout) ::   ldpyc 
    3363       LOGICAL,  DIMENSION(:,:), INTENT(inout) ::   ldflux      ! Surface flux extends below OSBL into MLE layer 
    3364       LOGICAL,  DIMENSION(:,:), INTENT(inout) ::   ldmle       ! MLE layer increases in thickness 
    33653350      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pwb_fk 
    33663351      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_mle      ! Temperature average over the depth of the MLE layer 
    … …  
    33973382      DO_2D( 0, 0, 0, 0 ) 
    33983383         ! 
    3399          IF ( ldconv(ji,jj) ) THEN 
     3384         IF ( l_conv(ji,jj) ) THEN 
    34003385            IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN 
    34013386               pt_mle(ji,jj) = ( pt_mle(ji,jj) * phmle(ji,jj) - pt_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 
    … …  
    34083393               zthermal = rab_n(ji,jj,1,jp_tem) 
    34093394               zbeta    = rab_n(ji,jj,1,jp_sal) 
    3410                DO jk = kbld(ji,jj), mld_prof(ji,jj) 
     3395               DO jk = nbld(ji,jj), mld_prof(ji,jj) 
    34113396                  zbuoy         = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) 
    34123397                  zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    … …  
    34323417      DO_2D( 0, 0, 0, 0 ) 
    34333418         ! 
    3434          IF ( ldconv(ji,jj) ) THEN 
     3419         IF ( l_conv(ji,jj) ) THEN 
    34353420            IF ( -2.0_wp * pwb_fk(ji,jj) / pwb_ent(ji,jj) > 0.5_wp ) THEN 
    34363421               IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN   ! MLE layer growing 
    34373422                  IF ( znd_param (ji,jj) > 100.0_wp ) THEN   ! Thermocline present 
    3438                      ldflux(ji,jj) = .FALSE. 
    3439                      ldmle(ji,jj)  = .FALSE. 
     3423                     l_flux(ji,jj) = .FALSE. 
     3424                     l_mle(ji,jj)  = .FALSE. 
    34403425                  ELSE   ! Thermocline not present 
    3441                      ldflux(ji,jj) = .TRUE. 
    3442                      ldmle(ji,jj)  = .TRUE. 
     3426                     l_flux(ji,jj) = .TRUE. 
     3427                     l_mle(ji,jj)  = .TRUE. 
    34433428                  ENDIF  ! znd_param > 100 
    34443429                  ! 
    34453430                  IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
    3446                      ldpyc(ji,jj) = .FALSE. 
     3431                     l_pyc(ji,jj) = .FALSE. 
    34473432                  ELSE 
    3448                      ldpyc(ji,jj) = .TRUE. 
     3433                     l_pyc(ji,jj) = .TRUE. 
    34493434                  ENDIF 
    34503435               ELSE   ! MLE layer restricted to OSBL or just below 
    34513436                  IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN   ! Weak stratification MLE layer can grow 
    3452                      ldpyc(ji,jj)  = .FALSE. 
    3453                      ldflux(ji,jj) = .TRUE. 
    3454                      ldmle(ji,jj)  = .TRUE. 
     3437                     l_pyc(ji,jj)  = .FALSE. 
     3438                     l_flux(ji,jj) = .TRUE. 
     3439                     l_mle(ji,jj)  = .TRUE. 
    34553440                  ELSE   ! Strong stratification 
    3456                      ldpyc(ji,jj)  = .TRUE. 
    3457                      ldflux(ji,jj) = .FALSE. 
    3458                      ldmle(ji,jj)  = .FALSE. 
     3441                     l_pyc(ji,jj)  = .TRUE. 
     3442                     l_flux(ji,jj) = .FALSE. 
     3443                     l_mle(ji,jj)  = .FALSE. 
    34593444                  END IF   ! pdb_bl < rn_mle_thresh_bl and 
    34603445               END IF   ! phmle > 1.2 phbl 
    34613446            ELSE 
    3462                ldpyc(ji,jj)  = .TRUE. 
    3463                ldflux(ji,jj) = .FALSE. 
    3464                ldmle(ji,jj)  = .FALSE. 
    3465                IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) ldpyc(ji,jj) = .FALSE. 
     3447               l_pyc(ji,jj)  = .TRUE. 
     3448               l_flux(ji,jj) = .FALSE. 
     3449               l_mle(ji,jj)  = .FALSE. 
     3450               IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
    34663451            END IF   !  -2.0 * pwb_fk(ji,jj) / pwb_ent > 0.5 
    34673452         ELSE   ! Stable Boundary Layer 
    3468             ldpyc(ji,jj)  = .FALSE. 
    3469             ldflux(ji,jj) = .FALSE. 
    3470             ldmle(ji,jj)  = .FALSE. 
    3471          END IF   ! ldconv 
     3453            l_pyc(ji,jj)  = .FALSE. 
     3454            l_flux(ji,jj) = .FALSE. 
     3455            l_mle(ji,jj)  = .FALSE. 
     3456         END IF   ! l_conv 
    34723457#ifdef key_osm_debug 
    34733458         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    34743459            WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',pwb_ent(ji,jj), & 
    34753460               & '  zhmle=',phmle(ji,jj), ' zhbl=', phbl(ji,jj), & 
    3476                & ' lpyc= ', ldpyc(ji,jj), ' lflux= ', ldflux(ji,jj),  ' lmle= ', ldmle(ji,jj), ' lconv= ', ldconv(ji,jj) 
     3461               & ' lpyc= ', l_pyc(ji,jj), ' lflux= ', l_flux(ji,jj),  ' lmle= ', l_mle(ji,jj), ' lconv= ', l_conv(ji,jj) 
    34773462            FLUSH(narea+100) 
    34783463         END IF 
    … …  
    34853470   END SUBROUTINE zdf_osm_osbl_state_fk 
    34863471 
    3487    SUBROUTINE zdf_osm_mle_parameters( Kmm, ldconv, ldmle, kmld_prof, pmld, phmle, pvel_mle, pdiff_mle, pdbds_mle, phbl, pb_bl, pwb0tot ) 
     3472   SUBROUTINE zdf_osm_mle_parameters( Kmm,       kmld_prof, pmld, phmle, pvel_mle,   & 
     3473      &                               pdiff_mle, pdbds_mle, phbl, pb_bl, pwb0tot ) 
    34883474      !!---------------------------------------------------------------------- 
    34893475      !!               ***  ROUTINE zdf_osm_mle_parameters  *** 
    … …  
    35003486      !!---------------------------------------------------------------------- 
    35013487      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index 
    3502       LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv      ! BL stability flags 
    3503       LOGICAL,  DIMENSION(:,:),    INTENT(inout) ::   ldmle       ! MLE layer increases in thickness 
    35043488      INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kmld_prof 
    35053489      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == ! 
    … …  
    35273511      ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE 
    35283512      DO_2D( 0, 0, 0, 0 ) 
    3529          IF ( ldconv(ji,jj) ) THEN 
     3513         IF ( l_conv(ji,jj) ) THEN 
    35303514            ztmp =  r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf 
    35313515            ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt 
    … …  
    35363520      ! Timestep mixed layer eddy depth 
    35373521      DO_2D( 0, 0, 0, 0 ) 
    3538          IF ( ldmle(ji,jj) ) THEN   ! MLE layer growing 
     3522         IF ( l_mle(ji,jj) ) THEN   ! MLE layer growing 
    35393523            ! Buoyancy gradient at base of MLE layer 
    35403524            zthermal = rab_n(ji,jj,1,jp_tem) 
Note: See TracChangeset for help on using the changeset viewer.