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 13858 for NEMO – NEMO

Changeset 13858 for NEMO


Ignore:
Timestamp:
2020-11-24T10:45:09+01:00 (3 years ago)
Author:
agn
Message:

first try of Alan's new code

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/cfgs/SHARED/namelist_ref

    r13400 r13858  
    10931093   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor 
    10941094   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4 
     1095   rn_osm_bl_thresh   = 5.e-5      !Threshold buoyancy for deepening of OSBL base 
    10951096   nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv 
    10961097   ln_dia_osm = .true.         ! output OSMOSIS-OBL variables 
    … …  
    11221123   rn_osm_mle_lat      = 20.       ! reference latitude (degrees) of MLE coef.                        (case rn_mle=1) 
    11231124   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
    1124    rn_osm_mle_thresh = 0.0005      ! delta b criterion used for FK MLD criterion 
     1125   rn_osm_mle_thresh  = 0.0005     ! delta b criterion used for FK MLE criterion 
    11251126   rn_osm_mle_tau     = 172800.    ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
    11261127   ln_osm_hmle_limit   = .false.   ! limit hmle to rn_osm_hmle_limit*hbl 

  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90

    r13684 r13858  
    5858   USE traqsr         ! details of solar radiation absorption 
    5959   USE zdfddm         ! double diffusion mixing (avs array) 
     60!   USE zdfmxl         ! mixed layer depth 
    6061   USE iom            ! I/O library 
    6162   USE lib_mpp        ! MPP library 
    … …  
    133134   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
    134135   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
     136   REAL(wp) ::   rn_osm_bl_thresh          ! Threshold buoyancy for deepening of OSBL base. 
    135137   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
    136138 
    … …  
    249251      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    250252      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     253      REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 
    251254 
    252255 
    … …  
    262265      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    263266      LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl 
     267      LOGICAL, DIMENSION(jpi,jpj)  :: lshear    ! Shear layers 
     268      LOGICAL, DIMENSION(jpi,jpj)  :: lpyc      ! OSBL pycnocline present 
     269      LOGICAL, DIMENSION(jpi,jpj)  :: lflux     ! surface flux extends below OSBL into MLE layer. 
     270      LOGICAL, DIMENSION(jpi,jpj)  :: lmle      ! MLE layer increases in hickness. 
    264271 
    265272      ! mixed-layer variables 
    … …  
    267274      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 
    268275      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 
     276      INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level 
     277      INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 
    269278 
    270279      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients 
    … …  
    279288      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    280289      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
    281       REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure. 
    282       REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients 
    283  
     290      REAL(wp), DIMENSION(jpi,jpj) :: zddhdt                                    ! correction to dhdt due to internal structure. 
     291      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients 
     292      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext              ! external temperature/salinity and buoyancy gradients 
    284293      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
    285294 
    286       REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    287       REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
    288       REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 
    289       REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 
    290       REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
    291       REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
     295      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl  ! averages over the depth of the blayer 
     296      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml  ! averages over the depth of the mixed layer 
     297      REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle  ! averages over the depth of the MLE layer 
     298      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 
     299      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 
     300      REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer 
     301!      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
     302      REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
     303      REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
    292304      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline 
    293305      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline 
    … …  
    295307      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline 
    296308      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline 
    297  
     309      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
    298310      ! Flux-gradient relationship variables 
     311      REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 
    299312 
    300313      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale. 
    301314 
    302       REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity 
     315      REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.   
    303316      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. 
     317      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ 
    304318      REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. 
    305319      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep 
    … …  
    316330      REAL(wp) :: zthick, zz0, zz1 ! temporary variables 
    317331      REAL(wp) :: zvel_max, zhbl_s ! temporary variables 
    318       REAL(wp) :: zfac             ! temporary variable 
     332      REAL(wp) :: zfac, ztmp       ! temporary variable 
    319333      REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift 
    320334      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity 
    321335      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    322  
     336      REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 
     337      REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. 
     338      REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc 
     339      REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max 
    323340      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
    324       REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 
     341      REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 
    325342      REAL(wp) :: zzeta_s = 0._wp 
    326343      REAL(wp) :: zzeta_v = 0.46 
    … …  
    340357      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 
    341358      zhol(:,:)   = 0._wp 
    342       lconv(:,:)  = .FALSE. 
     359      lconv(:,:)  = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE. 
    343360      ! mixed layer 
    344361      ! no initialization of zhbl or zhml (or zdh?) 
    345362      zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp 
    346       zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    347       zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
    348  
    349       zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    350       zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     363      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp 
     364      zv_bl(:,:)   = 0._wp ; zb_bl(:,:)  = 0._wp 
     365      zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     366      zt_mle(:,:)   = 0._wp ; zs_mle(:,:)    = 0._wp ; zu_mle(:,:)   = 0._wp 
     367      zb_mle(:,:) = 0._wp 
     368      zv_ml(:,:)   = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
     369      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
    351370      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp 
    352       zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp 
    353       zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp 
     371      zdb_ml(:,:)  = 0._wp 
     372      zdt_mle(:,:)  = 0._wp ; zds_mle(:,:)  = 0._wp ; zdu_mle(:,:)   = 0._wp 
     373      zdv_mle(:,:)  = 0._wp ; zdb_mle(:,:)  = 0._wp 
     374      zwth_ent = 0._wp ; zws_ent = 0._wp 
    354375      ! 
    355376      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 
    356377      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    357378      ! 
    358       zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 
     379      zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 
    359380 
    360381      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
    … …  
    367388 
    368389      ! Flux-Gradient arrays. 
    369       zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
    370       zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp 
    371390      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp 
    372391      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp 
    … …  
    376395      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    377396 
    378       zdhdt_2(:,:) = 0._wp 
     397      zddhdt(:,:) = 0._wp 
    379398      ! hbl = MAX(hbl,epsln) 
    380399      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    548567              zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 
    549568              zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 
    550               lconv(ji,jj) = .TRUE. 
    551            ELSE 
     569            ELSE 
    552570              zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln ) 
    553               lconv(ji,jj) = .FALSE. 
    554571           ENDIF 
    555572        END DO 
    … …  
    584601            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 
    585602            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     603            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    586604         END DO 
    587605      END DO 
    588606      ! Averages over well-mixed and boundary layer 
    589       ! Alan: do we need zb_nl?, zb_ml? 
    590       CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
    591       CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
    592 ! External gradient 
    593       CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
    594  
     607      jp_ext(:,:) = 2 
     608      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     609!      jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 
     610      CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     611! Velocity components in frame aligned with surface stress. 
     612      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     613      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
     614! Determine the state of the OSBL, stable/unstable, shear/no shear 
     615      CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
    595616 
    596617      IF ( ln_osm_mle ) THEN 
    597          CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
    598          CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
    599        ENDIF 
    600  
     618! Fox-Kemper Scheme 
     619         mld_prof = 4 
     620         DO jk = 5, jpkm1 
     621           DO jj = 2, jpjm1 
     622             DO ji = 2, jpim1 
     623               IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     624             END DO 
     625           END DO 
     626         END DO 
     627         jp_ext_mle(:,:) = 2 
     628        CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) 
     629 
     630         DO jj = 2, jpjm1 
     631           DO ji = 2, jpim1 
     632              zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     633           END DO 
     634         END DO 
     635 
     636!! External gradient 
     637         CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     638         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
     639         CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 
     640         CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
     641         CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
     642      ELSE    ! ln_osm_mle 
     643! FK not selected, Boundary Layer only. 
     644         lpyc(:,:) = .TRUE. 
     645         lflux(:,:) = .FALSE. 
     646         lmle(:,:) = .FALSE. 
     647         DO jj = 2, jpjm1 
     648           DO ji = 2, jpim1 
     649             IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 
     650           END DO 
     651         END DO 
     652      ENDIF   ! ln_osm_mle 
     653 
     654! Test if pycnocline well resolved 
     655      DO jj = 2, jpjm1 
     656        DO ji = 2,jpim1 
     657          IF (lconv(ji,jj) ) THEN 
     658             ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 
     659             IF ( ztmp > 6 ) THEN 
     660      ! pycnocline well resolved 
     661               jp_ext(ji,jj) = 1 
     662             ELSE 
     663      ! pycnocline poorly resolved 
     664               jp_ext(ji,jj) = 0 
     665             ENDIF 
     666          ELSE 
     667      ! Stable conditions 
     668            jp_ext(ji,jj) = 0 
     669          ENDIF 
     670        END DO 
     671      END DO 
     672 
     673      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     674!      jp_ext = ibld-imld+1 
     675      CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
    601676! Rate of change of hbl 
    602       CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
    603       ! Calculate averages over depth of boundary layer 
    604          DO jj = 2, jpjm1 
    605             DO ji = 2, jpim1 
    606                zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 
     677      CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     678      DO jj = 2, jpjm1 
     679        DO ji = 2, jpim1 
     680          zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 
    607681               ! adjustment to represent limiting by ocean bottom 
    608                zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
    609             END DO 
     682          IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN 
     683             zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
     684             lpyc(ji,jj) = .FALSE. 
     685          ENDIF 
    610686         END DO 
     687       END DO 
    611688 
    612689      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index 
    … …  
    626703! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    627704! 
    628       CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
    629       ! Alan: do we need zb_ml? 
    630       CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     705      CALL zdf_osm_timestep_hbl( zdhdt ) 
     706! is external level in bounds? 
     707 
     708      CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    631709! 
    632710! 
     711! Check to see if lpyc needs to be changed  
     712 
    633713      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
     714 
     715      DO jj = 2, jpjm1 
     716        DO ji = 2, jpim1 
     717          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE.  
     718        END DO 
     719      END DO         
     720 
    634721      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    635722! 
    636723    ! Average over the depth of the mixed layer in the convective boundary layer 
    637     ! Alan: do we need zb_ml? 
    638     CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
     724!      jp_ext = ibld - imld +1 
     725      CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    639726    ! rotate mean currents and changes onto wind align co-ordinates 
    640727    ! 
    641     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
    642     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    643      zuw_bse = 0._wp 
    644      zvw_bse = 0._wp 
    645      zwth_ent = 0._wp 
    646      zws_ent = 0._wp 
    647      DO jj = 2, jpjm1 
    648         DO ji = 2, jpim1 
    649            IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 
    650               IF ( lconv(ji,jj) ) THEN 
    651              zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 
    652                       &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 
    653                       &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 
    654             zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 
    655                       &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 
    656                  IF ( zdb_ml(ji,jj) > 0._wp ) THEN 
    657                     zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    658                     zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    659                  ENDIF 
    660               ELSE 
    661                  zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    662                  zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    663               ENDIF 
    664            ENDIF 
    665         END DO 
    666      END DO 
    667  
     728     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     729     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    668730      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    669731      !  Pycnocline gradients for scalars and velocity 
    670732      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    671733 
    672       CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
    673       CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 
     734      CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     735      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 
    674736      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    675737       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    676738       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    677739       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    678  
    679        DO jj = 2, jpjm1 
    680           DO ji = 2, jpim1 
    681              IF ( lconv(ji,jj) ) THEN 
    682                zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    683                zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 
    684                zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    685                zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    686                zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    687                zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 
    688              ELSE 
    689                zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
    690                zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
    691              END IF 
    692           END DO 
    693        END DO 
    694 ! 
    695        DO jj = 2, jpjm1 
    696           DO ji = 2, jpim1 
    697              IF ( lconv(ji,jj) ) THEN 
    698                 DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
    699                     zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    700                     ! 
    701                     zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5 
    702                     ! 
    703                     zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) & 
    704                          &            *                                      ( 1.0 -               0.5 * zznd_ml**2 ) 
    705                 END DO 
    706                 ! pycnocline - if present linear profile 
    707                 IF ( zdh(ji,jj) > 0._wp ) THEN 
    708                    zgamma_b = 6.0 
    709                    DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    710                        zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    711                        ! 
    712                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    713                        ! 
    714                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    715                    END DO 
    716                    IF ( ibld_ext == 0 ) THEN 
    717                        zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
    718                        zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
    719                    ELSE 
    720                        zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
    721                        zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
    722                    ENDIF 
    723                 ENDIF 
    724                 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    725                 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 
    726                 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 
    727                 ! could be taken out, take account of entrainment represents as a diffusivity 
    728                 ! should remove w from here, represents entrainment 
    729              ELSE 
    730              ! stable conditions 
    731                 DO jk = 2, ibld(ji,jj) 
    732                    zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    733                    zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
    734                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    735                 END DO 
    736  
    737                 IF ( ibld_ext == 0 ) THEN 
    738                    zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
    739                    zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
    740                 ELSE 
    741                    zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
    742                    zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
    743                 ENDIF 
    744              ENDIF   ! end if ( lconv ) 
    745              ! 
    746           END DO  ! end of ji loop 
    747        END DO     ! end of jj loop 
     740       CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 
    748741 
    749742       ! 
    … …  
    843836                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    844837                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    845                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     838                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 
    846839                   ! non-gradient buoyancy terms 
    847840                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    848841                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    849842                END DO 
    850              ELSE 
     843                 
     844                IF ( lpyc(ji,jj) ) THEN 
     845                  ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     846                  ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) 
     847                  zwth_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)                   
     848                  zws_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 
     849! Cubic profile used for buoyancy term 
     850                  za_cubic = 0.755 * ztau_sc_u(ji,jj) 
     851                  zb_cubic = 0.25 * ztau_sc_u(ji,jj) 
     852                  DO jk = 2, ibld(ji,jj) 
     853                    zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     854                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
     855 
     856                    ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
     857                  END DO 
     858! 
     859                  zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 
     860                  zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) ) 
     861! 
     862                  zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
     863! 
     864                  zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
     865! 
     866                  zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )  
     867                  DO jk = 2, ibld(ji,jj) 
     868                    zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     869                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     870! 
     871                    ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
     872                  END DO 
     873               ENDIF ! End of pycnocline                   
     874             ELSE ! lconv test - stable conditions 
    851875                DO jk = 2, ibld(ji,jj) 
    852876                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
    … …  
    886910       END DO           ! jj loop 
    887911 
     912       DO jj = 2, jpjm1 
     913         DO ji = 2, jpim1 
     914           IF ( lpyc(ji,jj) ) THEN 
     915             IF ( j_ddh(ji,jj) == 0 ) THEN 
     916! Place holding code. Parametrization needs checking for these conditions. 
     917               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     918               zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
     919               zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
     920             ELSE 
     921               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     922               zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
     923               zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
     924             ENDIF 
     925             zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 
     926             zc_cubic = zuw_bse - zd_cubic 
     927! need ztau_sc_u to be available. Change to array.  
     928             DO jk = imld(ji,jj), ibld(ji,jj) 
     929                zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     930                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     931             END DO 
     932             zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 
     933             zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 
     934             zc_cubic = zvw_bse - zd_cubic 
     935             DO jk = imld(ji,jj), ibld(ji,jj) 
     936               zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj) 
     937               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     938             END DO 
     939           ENDIF  ! lpyc 
     940         END DO ! ji loop 
     941       END DO  ! jj loop 
     942 
    888943       IF(ln_dia_osm) THEN 
    889944          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
    … …  
    892947! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    893948 
    894        WHERE ( lconv ) 
    895           zsc_wth_1 = zwth0 
    896           zsc_ws_1 = zws0 
    897        ELSEWHERE 
    898           zsc_wth_1 = 2.0 * zwthav 
    899           zsc_ws_1 = zws0 
    900        ENDWHERE 
     949       DO jj = 1, jpjm1 
     950          DO ji = 1, jpim1 
     951           
     952            IF ( lconv(ji,jj) ) THEN 
     953              zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 
     954              zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 
     955              IF ( lpyc(ji,jj) ) THEN 
     956! Pycnocline scales 
     957                 zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 
     958                 zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 
     959               ENDIF 
     960            ELSE 
     961              zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 
     962              zsc_ws_1(ji,jj) = zws0(ji,jj) 
     963            ENDIF 
     964          END DO 
     965        END DO 
    901966 
    902967       DO jj = 2, jpjm1 
    … …  
    915980                       &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) ) 
    916981               END DO 
     982! 
     983               IF ( lpyc(ji,jj) ) THEN 
     984! pycnocline 
     985                 DO jk = imld(ji,jj), ibld(ji,jj) 
     986                   zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     987                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )  
     988                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )  
     989                 END DO 
     990              ENDIF 
    917991            ELSE 
    918                DO jk = 2, ibld(ji,jj) 
    919                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    920                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    921                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    922                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
    923                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    924                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
    925                END DO 
     992               IF( zdhdt(ji,jj) > 0. ) THEN 
     993                 DO jk = 2, ibld(ji,jj) 
     994                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     995                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     996                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     997                 &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
     998                    ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     999                  &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
     1000                 END DO 
     1001               ENDIF 
    9261002            ENDIF 
    9271003          ENDDO    ! ji loop 
    … …  
    9841060! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
    9851061 
     1062 
     1063 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 
     1064 
    9861065      DO jj = 2, jpjm1 
    9871066         DO ji = 2, jpim1 
    988             IF ( lconv(ji,jj) ) THEN 
     1067            IF ( .not. lconv(ji,jj) ) THEN 
    9891068               DO jk = 2, ibld(ji,jj) 
    990                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    991                   IF ( znd >= 0.0 ) THEN 
    992                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    993                      ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    994                   ELSE 
    995                      ghamu(ji,jj,jk) = 0._wp 
    996                      ghamv(ji,jj,jk) = 0._wp 
    997                   ENDIF 
    998                END DO 
    999             ELSE 
    1000                DO jk = 2, ibld(ji,jj) 
    1001                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1069                  znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 
    10021070                  IF ( znd >= 0.0 ) THEN 
    10031071                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    … …  
    10121080      END DO 
    10131081 
    1014        IF(ln_dia_osm) THEN 
    1015           IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
    1016           IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
    1017        END IF 
    10181082      ! pynocline contributions 
    1019        ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
    1020        zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 
    10211083       DO jj = 2, jpjm1 
    10221084          DO ji = 2, jpim1 
    1023              IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1085            IF ( .not. lconv(ji,jj) ) THEN 
     1086             IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    10241087                DO jk= 2, ibld(ji,jj) 
    10251088                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    … …  
    10271090                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    10281091                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1029                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
    10301092                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    10311093                END DO 
    10321094             END IF 
     1095            END IF 
    10331096          END DO 
    10341097       END DO 
    1035  
    1036 ! Entrainment contribution. 
     1098      IF(ln_dia_osm) THEN 
     1099          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     1100          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     1101       END IF 
    10371102 
    10381103       DO jj=2, jpjm1 
    10391104          DO ji = 2, jpim1 
    1040             IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 
    1041                DO jk = 1, imld(ji,jj) - 1 
    1042                   znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1043                   ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1044                   ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1045                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1046                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
    1047                END DO 
    1048                DO jk = imld(ji,jj), ibld(ji,jj) 
    1049                   znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    1050                   ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1051                   ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1052                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1053                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1054                 END DO 
    1055              ENDIF 
    1056  
    10571105             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    10581106             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    … …  
    10651113          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
    10661114          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
    1067           IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 
    1068           IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 
    10691115          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
    10701116          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
    … …  
    11551201          DO jj = 2 , jpjm1 
    11561202             DO ji = 2, jpim1 
    1157                 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 
    1158             ! Calculate MLE flux profiles 
    1159                    ! DO jk = 1, mld_prof(ji,jj) 
    1160                    !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
    1161                    !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 
    1162                    !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
    1163                    !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 
    1164                    !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
    1165                    ! END DO 
    1166             ! Calculate MLE flux contribution from surface fluxes 
     1203                 IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 
     1204                ! Calculate MLE flux contribution from surface fluxes 
    11671205                   DO jk = 1, ibld(ji,jj) 
    11681206                     znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 
    … …  
    11761214                    END DO 
    11771215            ! Viscosity for MLEs 
    1178                     DO jk = ibld(ji,jj), mld_prof(ji,jj) 
    1179                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 
     1216                    DO jk = 1, mld_prof(ji,jj) 
     1217                      znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1218                      zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 
    11801219                    END DO 
    1181             ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 
    1182                     jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 
    1183                     jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 
    1184                     DO jk = mld_prof(ji,jj),  jl 
    1185                        zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 
    1186                             & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 
    1187                             & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 
     1220                 ELSE 
     1221! Surface transports limited to OSBL.                  
     1222            ! Viscosity for MLEs 
     1223                    DO jk = 1, mld_prof(ji,jj) 
     1224                      znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1225                      zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 
    11881226                    END DO 
    1189                 ENDIF 
     1227                 ENDIF 
    11901228            END DO 
    11911229          END DO 
    … …  
    13071345 
    13081346CONTAINS 
    1309  
    1310  
    1311    ! Alan: do we need zb? 
    1312    SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1347! subroutine code changed, needs syntax checking. 
     1348  SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 
     1349 
     1350!!--------------------------------------------------------------------- 
     1351     !!                   ***  ROUTINE zdf_osm_diffusivity_viscosity  *** 
     1352     !! 
     1353     !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. 
     1354     !! 
     1355     !! ** Method  :  
     1356     !! 
     1357     !! !!---------------------------------------------------------------------- 
     1358     REAL(wp), DIMENSION(:,:,:) :: zdiffut 
     1359     REAL(wp), DIMENSION(:,:,:) :: zviscos 
     1360! local 
     1361 
     1362! Scales used to calculate eddy diffusivity and viscosity profiles 
     1363      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 
     1364      REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 
     1365      REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 
     1366      REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 
     1367! 
     1368      REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 
     1369       
     1370      REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 
     1371      REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 
     1372      REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 
     1373       
     1374      DO jj = 2, jpjm1 
     1375          DO ji = 2, jpim1 
     1376             IF ( lconv(ji,jj) ) THEN 
     1377              
     1378               zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 
     1379               zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1380               zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 
     1381 
     1382               zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 
     1383               zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 
     1384 
     1385               IF ( lpyc(ji,jj) ) THEN 
     1386                 zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 
     1387 
     1388                 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     1389                   zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 
     1390                 ENDIF 
     1391                
     1392                 zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 
     1393                 zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 
     1394                 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 
     1395                  
     1396                 zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 
     1397                 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 
     1398                 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     1399                   zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 
     1400                 ENDIF 
     1401                 
     1402                 zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 
     1403                 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 
     1404                 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 
     1405 
     1406                 zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third 
     1407                 zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 
     1408               ELSE 
     1409                 zbeta_d_sc = 1.0 
     1410                 zbeta_v_sc = 1.0 
     1411               ENDIF 
     1412             ELSE 
     1413               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1414               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1415             END IF 
     1416          END DO 
     1417       END DO 
     1418! 
     1419       DO jj = 2, jpjm1 
     1420          DO ji = 2, jpim1 
     1421             IF ( lconv(ji,jj) ) THEN 
     1422                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
     1423                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
     1424                    ! 
     1425                    zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
     1426                    ! 
     1427                    zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 
     1428      &            *                                      ( 1.0 - 0.5 * zznd_ml**2 ) 
     1429                END DO 
     1430! pycnocline 
     1431                IF ( lpyc(ji,jj) ) THEN 
     1432! Diffusivity profile in the pycnocline given by cubic polynomial. 
     1433                   za_cubic = 0.5 
     1434                   zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
     1435                   zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & 
     1436                        & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 
     1437                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic  - zb_cubic ) 
     1438                   zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1439                   DO jk = imld(ji,jj) , ibld(ji,jj) 
     1440                     zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1441                         ! 
     1442                     zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 +   zd_cubic * zznd_pyc**3 ) 
     1443 
     1444                     zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) 
     1445                   END DO 
     1446 ! viscosity profiles. 
     1447                   za_cubic = 0.5 
     1448                   zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
     1449                   zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 
     1450                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) 
     1451                   zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1452                   DO jk = imld(ji,jj) , ibld(ji,jj) 
     1453                      zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1454                       zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 
     1455                       zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) 
     1456                   END DO 
     1457                   IF ( zdhdt(ji,jj) > 0._wp ) THEN 
     1458                    zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
     1459                    zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
     1460                   ELSE 
     1461                     zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     1462                     zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     1463                   ENDIF 
     1464                ENDIF 
     1465             ELSE 
     1466             ! stable conditions 
     1467                DO jk = 2, ibld(ji,jj) 
     1468                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1469                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
     1470                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
     1471                END DO 
     1472 
     1473                IF ( zdhdt(ji,jj) > 0._wp ) THEN 
     1474                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 
     1475                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 
     1476                ENDIF 
     1477             ENDIF   ! end if ( lconv ) 
     1478             ! 
     1479          END DO  ! end of ji loop 
     1480       END DO     ! end of jj loop 
     1481        
     1482  END SUBROUTINE zdf_osm_diffusivity_viscosity 
     1483   
     1484  SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
     1485 
     1486!!--------------------------------------------------------------------- 
     1487     !!                   ***  ROUTINE zdf_osm_osbl_state  *** 
     1488     !! 
     1489     !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number 
     1490     !! 
     1491     !! ** Method  :  
     1492     !! 
     1493     !! !!---------------------------------------------------------------------- 
     1494 
     1495     INTEGER, DIMENSION(jpi,jpj) :: j_ddh  ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low. 
     1496      
     1497     LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear 
     1498 
     1499     REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 
     1500     REAL(wp), DIMENSION(jpi,jpj) :: zshear  ! production of TKE due to shear across the pycnocline 
     1501     REAL(wp), DIMENSION(jpi,jpj) :: zri_i  ! Interfacial Richardson Number 
     1502 
     1503! Local Variables 
     1504 
     1505     INTEGER :: jj, ji 
     1506      
     1507     REAL(wp), DIMENSION(jpi,jpj) :: zekman 
     1508     REAL(wp) :: zri_p, zri_b   ! Richardson numbers 
     1509     REAL(wp) :: zshear_u, zshear_v, zwb_shr 
     1510     REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1511 
     1512     REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 
     1513     REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 
     1514     REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04      
     1515     REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 
     1516     REAL, PARAMETER :: rn_ri_p_thresh = 27.0 
     1517     REAL, PARAMETER :: zrot=0._wp  ! dummy rotation rate of surface stress. 
     1518      
     1519! Determins stability and set flag lconv 
     1520     DO jj = 2, jpjm1 
     1521       DO ji = 2, jpim1 
     1522         IF ( zhol(ji,jj) < 0._wp ) THEN 
     1523            lconv(ji,jj) = .TRUE. 
     1524          ELSE 
     1525             lconv(ji,jj) = .FALSE. 
     1526          ENDIF 
     1527       END DO 
     1528     END DO        
     1529  
     1530     zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 
     1531      
     1532     WHERE ( lconv ) 
     1533       zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 
     1534     END WHERE 
     1535 
     1536     zshear(:,:) = 0._wp 
     1537     j_ddh(:,:) = 1      
     1538  
     1539     DO jj = 2, jpjm1 
     1540       DO ji = 2, jpim1 
     1541         IF ( lconv(ji,jj) ) THEN 
     1542            IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1543              zri_p = MAX (  SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) )  *  ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 
     1544                   & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp ) 
     1545             
     1546              zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 
     1547                         
     1548              zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 
     1549!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1550! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  ! 
     1551! full code available                                          ! 
     1552!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1553              IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 
     1554! Growing shear layer 
     1555                j_ddh(ji,jj) = 0 
     1556                lshear(ji,jj) = .TRUE. 
     1557              ELSE 
     1558                j_ddh(ji,jj) = 1 
     1559                IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 
     1560! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 
     1561                  lshear(ji,jj) = .TRUE. 
     1562                ELSE 
     1563! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 
     1564                  zshear(ji,jj) = 0.5 * zshear(ji,jj) 
     1565                  lshear(ji,jj) = .FALSE. 
     1566                ENDIF  
     1567              ENDIF                 
     1568            ELSE                ! zdb_bl test, note zshear set to zero 
     1569              j_ddh(ji,jj) = 2 
     1570              lshear(ji,jj) = .FALSE. 
     1571            ENDIF 
     1572          ENDIF 
     1573       END DO 
     1574     END DO 
     1575  
     1576! Calculate entrainment buoyancy flux due to surface fluxes. 
     1577 
     1578     DO jj = 2, jpjm1 
     1579       DO ji = 2, jpim1 
     1580         IF ( lconv(ji,jj) ) THEN 
     1581           zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1582           zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1583           zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1584           zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
     1585           IF (nn_osm_SD_reduce > 0 ) THEN 
     1586           ! Effective Stokes drift already reduced from surface value 
     1587              zr_stokes = 1.0_wp 
     1588           ELSE 
     1589            ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 
     1590             ! requires further reduction where BL is deep 
     1591              zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1592            &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1593           END IF 
     1594           zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1595                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1596                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1597                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1598             ! 
     1599         ENDIF 
     1600       END DO ! ji loop 
     1601     END DO   ! jj loop  
     1602 
     1603     zwb_min(:,:) = 0._wp 
     1604 
     1605     DO jj = 2, jpjm1 
     1606       DO ji = 2, jpjm1 
     1607         IF ( lshear(ji,jj) ) THEN 
     1608           IF ( lconv(ji,jj) ) THEN 
     1609! Unstable OSBL 
     1610              zwb_shr = -za_wb_s * zshear(ji,jj) 
     1611              IF ( j_ddh(ji,jj) == 0 ) THEN 
     1612 
     1613! Developing shear layer, additional shear production possible. 
     1614 
     1615                zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) /  zhbl(ji,jj), 0._wp ) 
     1616                zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 
     1617                zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 
     1618                 
     1619                zwb_shr = -za_wb_s * zshear(ji,jj) 
     1620                 
     1621              ENDIF                 
     1622              zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
     1623              zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1624           ELSE    ! IF ( lconv ) THEN - ENDIF 
     1625! Stable OSBL  - shear production not coded for first attempt.            
     1626           ENDIF  ! lconv 
     1627         ELSE  ! lshear 
     1628           IF ( lconv(ji,jj) ) THEN 
     1629! Unstable OSBL 
     1630              zwb_shr = -za_wb_s * zshear(ji,jj) 
     1631              zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
     1632              zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1633           ENDIF  ! lconv 
     1634         ENDIF    ! lshear 
     1635       END DO   ! ji 
     1636     END DO     ! jj 
     1637   END SUBROUTINE zdf_osm_osbl_state 
     1638      
     1639      
     1640   SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) 
    13131641     !!--------------------------------------------------------------------- 
    13141642     !!                   ***  ROUTINE zdf_vertical_average  *** 
    … …  
    13221650 
    13231651        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1652        INTEGER, DIMENSION(jpi,jpj) :: jp_ext 
    13241653 
    13251654        ! Alan: do we need zb? 
    1326         REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity 
     1655        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb        ! Average temperature and salinity 
    13271656        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
    13281657        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
    13291658        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
    13301659 
    1331         INTEGER :: jk, ji, jj 
     1660        INTEGER :: jk, ji, jj, ibld_ext 
    13321661        REAL(wp) :: zthick, zthermal, zbeta 
    13331662 
    … …  
    13581687            zu(ji,jj) = zu(ji,jj) / zthick 
    13591688            zv(ji,jj) = zv(ji,jj) / zthick 
    1360             ! Alan: do we need zb? 
    1361             zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 
    1362             zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 
    1363             zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 
    1364                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 
    1365             zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 
    1366                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 
    1367             zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1689            zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 
     1690            ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 
     1691            IF ( ibld_ext < mbkt(ji,jj) ) THEN 
     1692              zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem) 
     1693              zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal) 
     1694              zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) & 
     1695                     &    / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) 
     1696              zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) & 
     1697                     &   / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 
     1698              zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1699            ELSE 
     1700              zdt(ji,jj) = 0._wp 
     1701              zds(ji,jj) = 0._wp 
     1702              zdu(ji,jj) = 0._wp 
     1703              zdv(ji,jj) = 0._wp 
     1704              zdb(ji,jj) = 0._wp 
     1705            ENDIF 
    13681706         END DO 
    13691707        END DO 
    … …  
    13991737    END SUBROUTINE zdf_osm_velocity_rotation 
    14001738 
    1401     SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 
     1739    SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
     1740     !!--------------------------------------------------------------------- 
     1741     !!                   ***  ROUTINE zdf_osm_osbl_state_fk  *** 
     1742     !! 
     1743     !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme. 
     1744     !!  lpyc :: determines whether pycnocline flux-grad relationship needs to be determined 
     1745     !!  lflux :: determines whether effects of surface flux extend below the base of the OSBL 
     1746     !!  lmle  :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.  
     1747     !! 
     1748     !! ** Method  :  
     1749     !! 
     1750     !!  
     1751     !!---------------------------------------------------------------------- 
     1752       
     1753! Outputs 
     1754      LOGICAL,  DIMENSION(jpi,jpj)  :: lpyc, lflux, lmle 
     1755      REAL(wp), DIMENSION(jpi,jpj)  :: zwb_fk 
     1756! 
     1757      REAL(wp), DIMENSION(jpi,jpj)  :: znd_param 
     1758      REAL(wp)                      :: zbuoy, ztmp, zpe_mle_layer 
     1759      REAL(wp)                      :: zpe_mle_ref, zwb_ent, zdbdz_mle_int 
     1760       
     1761      znd_param(:,:) = 0._wp 
     1762 
     1763        DO jj = 2, jpjm1 
     1764          DO ji = 2, jpim1           
     1765             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     1766             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 
     1767          END DO 
     1768        END DO         
     1769        DO jj = 2, jpjm1 
     1770          DO ji = 2, jpim1 
     1771                    ! 
     1772            IF ( lconv(ji,jj) ) THEN 
     1773              IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 
     1774                zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1775                zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1776                zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1777                zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 
     1778! Calculate potential energies of actual profile and reference profile. 
     1779                zpe_mle_layer = 0._wp 
     1780                zpe_mle_ref = 0._wp 
     1781                DO jk = ibld(ji,jj), mld_prof(ji,jj) 
     1782                  zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) ) 
     1783                  zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     1784                  zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     1785                END DO 
     1786! Non-dimensional parameter to diagnose the presence of thermocline 
     1787                    
     1788                znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) 
     1789              ENDIF 
     1790            ENDIF 
     1791          END DO 
     1792        END DO 
     1793 
     1794! Diagnosis 
     1795        DO jj = 2, jpjm1 
     1796           DO ji = 2, jpim1 
     1797             IF ( lconv(ji,jj) ) THEN 
     1798               zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1799                  &                  - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1800                  &         + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 
     1801                  &         - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1802               IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 
     1803                 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 
     1804! MLE layer growing 
     1805                   IF ( znd_param (ji,jj) > 100. ) THEN 
     1806! Thermocline present 
     1807                     lflux(ji,jj) = .FALSE. 
     1808                     lmle(ji,jj) =.FALSE. 
     1809                   ELSE 
     1810! Thermocline not present 
     1811                     lflux(ji,jj) = .TRUE. 
     1812                     lmle(ji,jj) = .TRUE. 
     1813                   ENDIF  ! znd_param > 100 
     1814! 
     1815                   IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
     1816                     lpyc(ji,jj) = .FALSE. 
     1817                   ELSE 
     1818                      lpyc = .TRUE. 
     1819                   ENDIF 
     1820                 ELSE 
     1821! MLE layer restricted to OSBL or just below. 
     1822                   IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
     1823! Weak stratification MLE layer can grow. 
     1824                     lpyc(ji,jj) = .FALSE. 
     1825                     lflux(ji,jj) = .TRUE. 
     1826                     lmle(ji,jj) = .TRUE. 
     1827                   ELSE 
     1828! Strong stratification 
     1829                     lpyc(ji,jj) = .TRUE. 
     1830                     lflux(ji,jj) = .FALSE. 
     1831                     lmle(ji,jj) = .FALSE. 
     1832                   ENDIF ! zdb_bl < rn_mle_thresh_bl and  
     1833                 ENDIF  ! zhmle > 1.2 zhbl 
     1834               ELSE 
     1835                 lpyc(ji,jj) = .TRUE. 
     1836                 lflux(ji,jj) = .FALSE. 
     1837                 lmle(ji,jj) = .FALSE. 
     1838                 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 
     1839               ENDIF !  -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5  
     1840             ELSE 
     1841! Stable Boundary Layer 
     1842               lpyc(ji,jj) = .FALSE. 
     1843               lflux(ji,jj) = .FALSE. 
     1844               lmle(ji,jj) = .FALSE. 
     1845             ENDIF  ! lconv 
     1846           END DO 
     1847        END DO 
     1848    END SUBROUTINE zdf_osm_osbl_state_fk 
     1849 
     1850    SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 
    14021851     !!--------------------------------------------------------------------- 
    14031852     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
    … …  
    14091858     !!---------------------------------------------------------------------- 
    14101859 
     1860     INTEGER, DIMENSION(jpi,jpj)  :: jbase 
    14111861     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
    14121862 
    … …  
    14171867     DO jj = 2, jpjm1 
    14181868        DO ji = 2, jpim1 
    1419            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1869           IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
    14201870              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    14211871              zbeta    = rab_n(ji,jj,1,jp_sal) 
    1422               jkb = ibld(ji,jj) + ibld_ext 
     1872              jkb = jbase(ji,jj) 
    14231873              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
    14241874              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 
    … …  
    14271877                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
    14281878              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1879           ELSE 
     1880              zdtdz(ji,jj) = 0._wp 
     1881              zdsdz(ji,jj) = 0._wp 
     1882              zdbdz(ji,jj) = 0._wp 
    14291883           END IF 
    14301884        END DO 
    … …  
    14321886    END SUBROUTINE zdf_osm_external_gradients 
    14331887 
    1434     SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 
     1888    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 
    14351889 
    14361890     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1891     REAL(wp), DIMENSION(jpi,jpj) :: zalpha 
    14371892 
    14381893     INTEGER :: jk, jj, ji 
    14391894     REAL(wp) :: ztgrad, zsgrad, zbgrad 
    1440      REAL(wp) :: zgamma_b_nd, zgamma_c, znd 
    1441      REAL(wp) :: zzeta_s=0.3, ztmp 
     1895     REAL(wp) :: zgamma_b_nd, znd 
     1896     REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 
     1897     REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 
    14421898 
    14431899     DO jj = 2, jpjm1 
    14441900        DO ji = 2, jpim1 
    1445            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1901           IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    14461902              IF ( lconv(ji,jj) ) THEN  ! convective conditions 
    1447                     IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 
    1448                          & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 
    1449                          &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0 
    1450                        ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
    1451                        ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 
    1452                        zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 
    1453                        zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 
    1454                        zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
    1455                        zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 
    1456                             & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 
    1457                        DO jk = 2, ibld(ji,jj)+ibld_ext 
    1458                           znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 
    1459                           IF ( znd <= zzeta_s ) THEN 
    1460                              zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 
    1461                                   & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
    1462                              zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 
    1463                                   & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
    1464                              zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 
    1465                                   & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
    1466                           ELSE 
    1467                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
    1468                              zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
    1469                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
    1470                           ENDIF 
    1471                        END DO 
    1472                     ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 
     1903                IF ( lpyc(ji,jj) ) THEN 
     1904                   zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
     1905                   zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 
     1906                   zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 
     1907 
     1908                   ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
     1909!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1910! Commented lines in this section are not needed in new code, once tested ! 
     1911! can be removed                                                          ! 
     1912!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1913!                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
     1914!                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
     1915                   zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 
     1916                   zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
     1917                   DO jk = 2, ibld(ji,jj)+ibld_ext 
     1918                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 
     1919                     IF ( znd <= zzeta_m ) THEN 
     1920!                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 
     1921!                &        EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1922!                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 
     1923!                                  & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1924                        zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 
     1925                                  & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1926                     ELSE 
     1927!                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1928!                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1929                         zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1930                     ENDIF 
     1931                  END DO 
     1932               ENDIF ! if no pycnocline pycnocline gradients set to zero 
    14731933              ELSE 
    14741934                 ! stable conditions 
    14751935                 ! if pycnocline profile only defined when depth steady of increasing. 
    1476                  IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     1936                 IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady. 
    14771937                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    14781938                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    … …  
    15021962                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    15031963              ENDIF          ! IF (lconv) 
    1504            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1964           ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 
    15051965        END DO 
    15061966     END DO 
    … …  
    15271987         DO ji = 2, jpim1 
    15281988            ! 
    1529             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1989            IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    15301990               IF ( lconv (ji,jj) ) THEN 
    1531                   ! Unstable conditions 
    1532                   zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
    1533                        &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
    1534                        &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1991                  ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
     1992!                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1993!                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1994!                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
    15351995                  !Alan is this right? 
    1536                   zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
    1537                        &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
    1538                        &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
    1539                        &      )/ (zdh(ji,jj)  + epsln ) 
    1540                   DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
    1541                      znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
    1542                      IF ( znd <= 0.0 ) THEN 
    1543                         zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
    1544                         zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
    1545                      ELSE 
    1546                         zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
    1547                         zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
    1548                      ENDIF 
    1549                   END DO 
     1996!                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1997!                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1998!                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1999!                       &      )/ (zdh(ji,jj)  + epsln ) 
     2000!                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     2001!                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     2002!                     IF ( znd <= 0.0 ) THEN 
     2003!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     2004!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     2005!                     ELSE 
     2006!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     2007!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     2008!                     ENDIF 
     2009!                  END DO 
    15502010               ELSE 
    15512011                  ! stable conditions 
    … …  
    15682028    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
    15692029 
    1570     SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
     2030   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
    15712031     !!--------------------------------------------------------------------- 
    15722032     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
    … …  
    15782038     !!---------------------------------------------------------------------- 
    15792039 
    1580     REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl 
     2040    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt        ! Rate of change of hbl 
    15812041 
    15822042    INTEGER :: jj, ji 
    1583     REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 
    1584     REAL(wp) :: zvel_max, zwb_min 
    1585     REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     2043    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 
     2044    REAL(wp) :: zvel_max!, zwb_min 
    15862045    REAL(wp) :: zzeta_m = 0.3 
    15872046    REAL(wp) :: zgamma_c = 2.0 
    15882047    REAL(wp) :: zdhoh = 0.1 
    15892048    REAL(wp) :: alpha_bc = 0.5 
    1590  
    1591     DO jj = 2, jpjm1 
    1592        DO ji = 2, jpim1 
     2049    REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 
     2050  
     2051  DO jj = 2, jpjm1 
     2052     DO ji = 2, jpim1 
     2053        
     2054       IF ( lshear(ji,jj) ) THEN 
    15932055          IF ( lconv(ji,jj) ) THEN    ! Convective 
    1594              ! Alan is this right?  Yes, it's a bit different from the previous relationship 
    1595              ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 
    1596              !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
    1597              zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
    1598              zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
    1599              zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
    1600              zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
    1601              IF (nn_osm_SD_reduce > 0 ) THEN 
    1602                 ! Effective Stokes drift already reduced from surface value 
    1603                 zr_stokes = 1.0_wp 
    1604              ELSE 
    1605                 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 
    1606                 ! requires further reduction where BL is deep 
    1607                 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
    1608                      &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
    1609              END IF 
    1610  
    1611              zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
    1612                   &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
    1613                   &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
    1614                   &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
    1615              ! 
    1616              zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 
    16172056 
    16182057             IF ( ln_osm_mle ) THEN 
    1619                   !  zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 
    1620                 ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL 
     2058 
    16212059                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     2060          ! Fox-Kemper buoyancy flux average over OSBL 
    16222061                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
    16232062                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
    … …  
    16282067                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
    16292068                   ! OSBL is deepening, entrainment > restratification 
    1630                    IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     2069                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
     2070! *** Used for shear Needs to be changed to work stabily 
     2071!                zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 
     2072!                zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 
     2073!                zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 
     2074!                za_1 = 1.0 / zgamma_b**2 - 0.017 
     2075!                za_2 = 1.0 / zgamma_b**3 - 0.0025 
     2076!                zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 
     2077                      zpsi = 0._wp 
     2078                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     2079                      zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) 
     2080                      IF ( j_ddh(ji,jj) == 1 ) THEN 
     2081                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 
     2082                           zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2083                        ELSE 
     2084                           zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2085                        ENDIF 
     2086! Relaxation to dh_ref = zari * hbl 
     2087                        zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     2088                         
     2089                      ELSE  ! j_ddh == 0 
     2090! Growing shear layer 
     2091                        zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     2092                      ENDIF ! j_ddh 
     2093                        zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 
     2094                   ELSE    ! zdb_bl >0 
     2095                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     2096                   ENDIF 
     2097                ELSE   ! zwb_min + 2*zwb_fk_b < 0 
     2098                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     2099                   zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     2100 
     2101 
     2102                ENDIF 
     2103 
     2104             ELSE 
     2105                ! Fox-Kemper not used. 
     2106 
     2107                  zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     2108                  &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 
     2109                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     2110                ! added ajgn 23 July as temporay fix 
     2111 
     2112             ENDIF  ! ln_osm_mle 
     2113 
     2114            ELSE    ! lconv - Stable 
     2115                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     2116                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     2117                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     2118                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     2119                ELSE 
     2120                    zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     2121                ENDIF 
     2122                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
     2123            ENDIF  ! lconv 
     2124       ELSE ! lshear 
     2125         IF ( lconv(ji,jj) ) THEN    ! Convective 
     2126 
     2127             IF ( ln_osm_mle ) THEN 
     2128 
     2129                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     2130          ! Fox-Kemper buoyancy flux average over OSBL 
     2131                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     2132                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     2133                ELSE 
     2134                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     2135                ENDIF 
     2136                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2137                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     2138                   ! OSBL is deepening, entrainment > restratification 
     2139                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
    16312140                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    16322141                   ELSE 
    … …  
    16432152                ! Fox-Kemper not used. 
    16442153 
    1645                   zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     2154                  zvel_max = -zwb_ent(ji,jj) / & 
    16462155                  &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 
    16472156                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    16482157                ! added ajgn 23 July as temporay fix 
    16492158 
    1650              ENDIF 
    1651  
    1652              zdhdt_2(ji,jj) = 0._wp 
    1653  
    1654                 ! commented out ajgn 23 July as temporay fix 
    1655 !                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
    1656 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 
    1657 !                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    1658 !                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 
    1659 !                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
    1660 !                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 
    1661 !                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 
    1662 !                     ! Alan no idea what this should be? 
    1663 !                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 
    1664 !                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 
    1665 !                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 
    1666 !                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 
    1667 !                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 
    1668 !                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 
    1669 !                 ENDIF 
     2159             ENDIF  ! ln_osm_mle 
     2160 
    16702161            ELSE                        ! Stable 
    16712162                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
    1672                 zdhdt_2(ji,jj) = 0._wp 
    16732163                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
    16742164                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    1675                     zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     2165                    zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 
    16762166                ELSE 
    1677                     zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     2167                    zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    16782168                ENDIF 
    16792169                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
    1680             ENDIF 
    1681          END DO 
    1682       END DO 
     2170            ENDIF  ! lconv 
     2171         ENDIF ! lshear 
     2172       END DO 
     2173     END DO 
    16832174    END SUBROUTINE zdf_osm_calculate_dhdt 
    16842175 
    1685     SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     2176    SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 
    16862177     !!--------------------------------------------------------------------- 
    16872178     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
    … …  
    16972188 
    16982189 
    1699     REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl. 
     2190    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt   ! rates of change of hbl. 
    17002191 
    17012192    INTEGER :: jk, jj, ji, jm 
    … …  
    17352226                    IF ( ln_osm_mle ) THEN 
    17362227                       zhbl_s = zhbl_s + MIN( & 
    1737                         & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     2228                        & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
    17382229                        & e3w_n(ji,jj,jm) ) 
    17392230                    ELSE 
    17402231                      zhbl_s = zhbl_s + MIN( & 
    1741                         & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     2232                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
    17422233                        & e3w_n(ji,jj,jm) ) 
    17432234                    ENDIF 
    17442235 
    1745                     zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
    1746  
     2236!                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     2237                    IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN 
     2238                      zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     2239                      lpyc(ji,jj) = .FALSE. 
     2240                    ENDIF 
    17472241                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
    17482242                 END DO 
    … …  
    17652259                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 
    17662260 
    1767                     zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     2261!                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     2262                    IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 
     2263                      zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     2264                      lpyc(ji,jj) = .FALSE. 
     2265                    ENDIF 
    17682266                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
    17692267                 END DO 
    … …  
    17962294 
    17972295      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
    1798       ! 
     2296       ! 
    17992297      INTEGER :: jj, ji 
    18002298      INTEGER :: inhml 
    1801       REAL(wp) :: zari, ztau, zddhdt 
    1802  
    1803  
    1804       DO jj = 2, jpjm1 
    1805          DO ji = 2, jpim1 
     2299      REAL(wp) :: zari, ztau, zdh_ref 
     2300      REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 
     2301 
     2302    DO jj = 2, jpjm1 
     2303       DO ji = 2, jpim1 
     2304 
     2305         IF ( lshear(ji,jj) ) THEN 
    18062306            IF ( lconv(ji,jj) ) THEN 
     2307              IF ( j_ddh(ji,jj) == 0 ) THEN 
     2308! ddhdt for pycnocline determined in osm_calculate_dhdt 
     2309                dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt 
     2310              ELSE 
     2311! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt  
     2312                IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 
     2313                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2314                ELSE 
     2315                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
     2316                ENDIF 
     2317                ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt ) 
     2318                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     2319                IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 
     2320              ENDIF 
     2321                
     2322            ELSE ! lconv 
     2323! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL  
     2324 
     2325               ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 
     2326               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     2327                  ! boundary layer deepening 
     2328                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     2329                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     2330                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     2331                          & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 ) 
     2332                     zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 
     2333                  ELSE 
     2334                     zdh_ref = 0.2 * hbl(ji,jj) 
     2335                  ENDIF 
     2336               ELSE     ! IF(dhdt < 0) 
     2337                  zdh_ref = 0.2 * hbl(ji,jj) 
     2338               ENDIF    ! IF (dhdt >= 0) 
     2339               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     2340               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse 
     2341               ! Alan: this hml is never defined or used -- do we need it? 
     2342            ENDIF 
     2343              
     2344         ELSE   ! lshear   
     2345! for lshear = .FALSE. calculate ddhdt here 
     2346 
     2347             IF ( lconv(ji,jj) ) THEN 
    18072348 
    18082349               IF( ln_osm_mle ) THEN 
    1809                   IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
     2350                  IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
    18102351                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
    1811                      IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     2352                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 
    18122353                        IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
    1813                            zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    1814                                 (1.0 + zdb_bl(ji,jj)**2 / MAX( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
     2354                           zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    18152355                        ELSE                                                     ! unstable 
    1816                            zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
    1817                                 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     2356                           zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    18182357                        ENDIF 
    18192358                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    1820                         zddhdt = zari * hbl(ji,jj) 
     2359                        zdh_ref = zari * hbl(ji,jj) 
    18212360                     ELSE 
    18222361                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    1823                         zddhdt = 0.2 * hbl(ji,jj) 
     2362                        zdh_ref = 0.2 * hbl(ji,jj) 
    18242363                     ENDIF 
    18252364                  ELSE 
    18262365                     ztau = 0.2 * hbl(ji,jj) /  MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    1827                      zddhdt = 0.2 * hbl(ji,jj) 
     2366                     zdh_ref = 0.2 * hbl(ji,jj) 
    18282367                  ENDIF 
    18292368               ELSE ! ln_osm_mle 
    1830                   IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     2369                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 
    18312370                     IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability 
    1832                         zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 
    1833                              (1.0 + zdb_bl(ji,jj)**2 /  MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
     2371                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    18342372                     ELSE                                                     ! unstable 
    1835                         zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 
    1836                              (1.0 + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 
     2373                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    18372374                     ENDIF 
    18382375                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    1839                      zddhdt = zari * hbl(ji,jj) 
     2376                     zdh_ref = zari * hbl(ji,jj) 
    18402377                  ELSE 
    18412378                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 
    1842                      zddhdt = 0.2 * hbl(ji,jj) 
     2379                     zdh_ref = 0.2 * hbl(ji,jj) 
    18432380                  ENDIF 
    18442381 
    18452382               END IF  ! ln_osm_mle 
    18462383 
    1847                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     2384               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     2385!               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 
     2386               IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 
    18482387               ! Alan: this hml is never defined or used 
    18492388            ELSE   ! IF (lconv) 
    … …  
    18552394                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    18562395                          & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 ) 
    1857                      zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
     2396                     zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 
    18582397                  ELSE 
    1859                      zddhdt = 0.2 * hbl(ji,jj) 
     2398                     zdh_ref = 0.2 * hbl(ji,jj) 
    18602399                  ENDIF 
    18612400               ELSE     ! IF(dhdt < 0) 
    1862                   zddhdt = 0.2 * hbl(ji,jj) 
     2401                  zdh_ref = 0.2 * hbl(ji,jj) 
    18632402               ENDIF    ! IF (dhdt >= 0) 
    1864                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
    1865                IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse 
    1866                ! Alan: this hml is never defined or used -- do we need it? 
     2403               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     2404               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse 
    18672405            ENDIF       ! IF (lconv) 
    1868  
    1869             hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
    1870             inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    1871             imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
    1872             zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    1873             zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    1874          END DO 
    1875       END DO 
     2406         ENDIF  ! lshear 
     2407  
     2408         hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     2409         inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)), 1.e-3) ) , 1 ) 
     2410         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     2411         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     2412         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     2413       END DO 
     2414     END DO 
    18762415 
    18772416    END SUBROUTINE zdf_osm_pycnocline_thickness 
    18782417 
    18792418 
    1880    SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     2419   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
    18812420      !!---------------------------------------------------------------------- 
    18822421      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  *** 
    … …  
    18912430 
    18922431      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 
     2432      REAL(wp), DIMENSION(jpi,jpj)     :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 
    18932433      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
    18942434      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy 
    … …  
    19772517      END DO 
    19782518 
     2519      DO jj = 2, jpjm1 
     2520        DO ji = 2, jpim1 
     2521           ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     2522           zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
     2523                & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
     2524        END DO 
     2525      END DO 
     2526       
    19792527 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
    1980   SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     2528  SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
    19812529      !!---------------------------------------------------------------------- 
    19822530      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
    … …  
    19892537      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
    19902538 
    1991       REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle 
     2539      INTEGER, DIMENSION(jpi,jpj)      :: mld_prof 
     2540      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 
    19922541      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    1993       INTEGER  ::   ii, ij, ik  ! local integers 
     2542      INTEGER  ::   ii, ij, ik, jkb, jkb1  ! local integers 
    19942543      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
    1995       REAL(wp) ::   zdb_mle, ztmp, zdbds_mle 
    1996  
    1997       mld_prof(:,:) = 4 
    1998       DO jk = 5, jpkm1 
    1999         DO jj = 2, jpjm1 
    2000           DO ji = 2, jpim1 
    2001             IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
    2002           END DO 
     2544      REAL(wp) ::  ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle 
     2545 
     2546   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
     2547 
     2548      DO jj = 2, jpjm1 
     2549        DO ji = 2, jpim1 
     2550          IF ( lconv(ji,jj) ) THEN 
     2551             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     2552      ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
     2553             zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
     2554             zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 
     2555          ENDIF 
    20032556        END DO 
    20042557      END DO 
    2005       ! DO jj = 2, jpjm1 
    2006       !    DO ji = 1, jpim1 
    2007       !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
    2008       !    END DO 
    2009       !  END DO 
    20102558   ! Timestep mixed layer eddy depth. 
    20112559      DO jj = 2, jpjm1 
    20122560        DO ji = 2, jpim1 
    2013           zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 
    2014           IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 
    2015              hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh )  ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 
    2016           ELSE 
    2017             ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 
    2018              ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 
    2019              !   hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
    2020              ! ELSE 
    2021              !   hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
    2022              ! ENDIF 
    2023              IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
    2024                hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
    2025              ELSE 
    2026                hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 
    2027              ENDIF 
    2028           ENDIF 
    2029           hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
    2030           IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*zmld(ji,jj)) 
     2561           IF ( lmle(ji,jj) ) THEN  ! MLE layer growing. 
     2562! Buoyancy gradient at base of MLE layer. 
     2563              zthermal = rab_n(ji,jj,1,jp_tem) 
     2564              zbeta    = rab_n(ji,jj,1,jp_sal) 
     2565              jkb = mld_prof(ji,jj) 
     2566              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     2567!               
     2568              zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) ) 
     2569              zdb_mle = zb_bl(ji,jj) - zbuoy  
     2570! Timestep hmle.  
     2571              hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / zdb_mle 
     2572           ELSE 
     2573              IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN 
     2574                 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     2575              ELSE 
     2576                 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau 
     2577              ENDIF 
     2578           ENDIF 
     2579           hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
     2580          IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 
    20312581        END DO 
    20322582      END DO 
    … …  
    20452595         END DO 
    20462596       END DO 
    2047    ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
    2048  
    2049       DO jj = 2, jpjm1 
    2050         DO ji = 2, jpim1 
    2051           IF ( lconv(ji,jj) ) THEN 
    2052              ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
    2053              ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 
    2054              !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 
    2055              ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 
    2056              !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 
    2057              zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
    2058                   & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
    2059              zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 
    2060       ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
    2061              zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
    2062              zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf 
    2063           ENDIF 
    2064         END DO 
    2065       END DO 
    20662597END SUBROUTINE zdf_osm_mle_parameters 
    20672598 
    … …  
    20882619          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 
    20892620          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
    2090           & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, ln_zdfosm_ice_shelter 
     2621          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 
    20912622! Namelist for Fox-Kemper parametrization. 
    20922623      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
    … …  
    21382669        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 
    21392670        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm 
     2671        WRITE(numout,*) '         Threshold used to define BL              rn_osm_bl_thresh  = ', rn_osm_bl_thresh, 'm^2/s' 
    21402672        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix 
    21412673        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty 
    … …  
    21442676        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv 
    21452677     ENDIF 
     2678 
     2679 
     2680     !                              ! Check wave coupling settings ! 
     2681     !                         ! Further work needed - see ticket #2447 ! 
     2682     IF( nn_osm_wave == 2 ) THEN 
     2683        IF (.NOT. ( ln_wave .AND. ln_sdw )) & 
     2684           & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 
     2685     END IF 
    21462686 
    21472687     !                              ! allocate zdfosm arrays 
    … …  
    21712711            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg' 
    21722712            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c 
    2173             WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
     2713            WRITE(numout,*) '         Threshold used to define MLE for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
    21742714            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
    21752715            WRITE(numout,*) '         switch to limit hmle                                      ln_osm_hmle_limit  = ', ln_osm_hmle_limit 
Note: See TracChangeset for help on using the changeset viewer.