Ignore:
Timestamp:
2021-02-26T20:30:58+01:00 (5 months ago)
Author:
smueller
Message:

Synchronisation of the OSMOSIS boundary layer scheme with the version developed in branch /NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0: adoption of changesets [14405,14407,14408,14412,14514] (ticket #2353)

File:
1 edited

Legend:

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

    r14551 r14554  
    285285      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    286286      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
    287       REAL(wp), DIMENSION(jpi,jpj) :: zddhdt                                    ! correction to dhdt due to internal structure. 
    288287      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients 
    289288      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext              ! external temperature/salinity and buoyancy gradients 
     
    299298      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
    300299      ! Flux-gradient relationship variables 
    301       REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 
     300      REAL(wp), DIMENSION(jpi, jpj) :: zshear   ! Shear production 
    302301 
    303302      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep 
     
    370369      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    371370 
    372       zddhdt(:,:) = 0._wp 
    373371      ! hbl = MAX(hbl,epsln) 
    374372      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    560558!      jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 
    561559      CALL zdf_osm_vertical_average( Kbb, Kmm,                                               & 
    562          &                           ibld, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, ibld-imld+1,   & 
     560         &                           imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, ibld-imld+1,   & 
    563561         &                           zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    564562! Velocity components in frame aligned with surface stress. 
     
    566564      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    567565! Determine the state of the OSBL, stable/unstable, shear/no shear 
    568       CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
     566      CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 
    569567 
    570568      IF ( ln_osm_mle ) THEN 
     
    586584         CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 
    587585         CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
    588          CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
     586         CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
    589587      ELSE    ! ln_osm_mle 
    590588! FK not selected, Boundary Layer only. 
     
    622620         &                           ibld-imld+1, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    623621! Rate of change of hbl 
    624       CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     622      CALL zdf_osm_calculate_dhdt( zdhdt ) 
    625623      DO_2D( 0, 0, 0, 0 ) 
    626624       zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it 
     
    685683      ! Calculate non-gradient components of the flux-gradient relationships 
    686684      ! -------------------------------------------------------------------- 
    687       CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, ibld_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear,   & 
     685      CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear,             & 
    688686         &                    zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla,     & 
    689687         &                    zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml,                & 
     
    928926              zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 
    929927 
    930               IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     928              IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 
    931929                zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 
    932930              ENDIF 
     
    938936              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) 
    939937              zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 
    940               IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 
     938              IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 
    941939                zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 
    942940              ENDIF 
     
    944942              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) ) ) 
    945943              zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 
    946               zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 
     944              zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) ) 
    947945 
    948946              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 
     
    10221020  END SUBROUTINE zdf_osm_diffusivity_viscosity 
    10231021 
    1024   SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 
     1022  SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 
    10251023 
    10261024!!--------------------------------------------------------------------- 
     
    10391037     REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 
    10401038     REAL(wp), DIMENSION(jpi,jpj) :: zshear  ! production of TKE due to shear across the pycnocline 
    1041      REAL(wp), DIMENSION(jpi,jpj) :: zri_i  ! Interfacial Richardson Number 
    10421039 
    10431040! Local Variables 
     
    10461043 
    10471044     REAL(wp), DIMENSION(jpi,jpj) :: zekman 
    1048      REAL(wp) :: zri_p, zri_b   ! Richardson numbers 
     1045     REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b   ! Richardson numbers 
    10491046     REAL(wp) :: zshear_u, zshear_v, zwb_shr 
    10501047     REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
    10511048 
    1052      REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 
    1053      REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 
    1054      REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 
     1049     REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8 
     1050     REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03 
    10551051     REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 
    10561052     REAL, PARAMETER :: rn_ri_p_thresh = 27.0 
     1053     REAL, PARAMETER :: zri_c = 0.25 
     1054     REAL, PARAMETER :: zek = 4.0 
    10571055     REAL, PARAMETER :: zrot=0._wp  ! dummy rotation rate of surface stress. 
    10581056 
     
    10671065     END_2D 
    10681066 
    1069      zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 
    1070  
    1071      WHERE ( lconv ) 
    1072        zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 
    1073      END WHERE 
     1067     zekman(:,:) = EXP( -1.0_wp * zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 
    10741068 
    10751069     zshear(:,:) = 0._wp 
     
    10791073      IF ( lconv(ji,jj) ) THEN 
    10801074         IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    1081            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 & 
    1082                 & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp ) 
    1083  
    1084            zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 
    1085  
    1086            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 ) ) 
     1075            zri_p(ji,jj) = 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 & 
     1076               & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp ) 
     1077 
     1078            IF ( ff_t(ji,jj) >= 0.0_wp ) THEN 
     1079               ! Northern hemisphere 
     1080               zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( -1.0_wp * zdv_ml(ji,jj), 1e-5_wp)**2 ) 
     1081            ELSE 
     1082               ! Southern hemisphere 
     1083               zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX(           zdv_ml(ji,jj), 1e-5_wp)**2 ) 
     1084            END IF 
     1085            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 ) ) 
     1086            ! Stability dependence 
     1087            zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) ) 
    10871088!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    10881089! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  ! 
    10891090! full code available                                          ! 
    10901091!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    1091            IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 
     1092            IF ( zshear(ji,jj) > 1e-10 ) THEN 
     1093               IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN 
    10921094! Growing shear layer 
    1093              j_ddh(ji,jj) = 0 
    1094              lshear(ji,jj) = .TRUE. 
    1095            ELSE 
    1096              j_ddh(ji,jj) = 1 
    1097              IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 
     1095                  j_ddh(ji,jj) = 0 
     1096                  lshear(ji,jj) = .TRUE. 
     1097               ELSE 
     1098                  j_ddh(ji,jj) = 1 
     1099!                 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 
    10981100! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 
    1099                lshear(ji,jj) = .TRUE. 
    1100              ELSE 
     1101                  lshear(ji,jj) = .TRUE. 
     1102!             ELSE 
     1103               END IF 
     1104            ELSE 
     1105               j_ddh(ji,jj) = 2 
     1106               lshear(ji,jj) = .FALSE. 
     1107            END IF 
    11011108! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 
    1102                zshear(ji,jj) = 0.5 * zshear(ji,jj) 
    1103                lshear(ji,jj) = .FALSE. 
    1104              ENDIF 
    1105            ENDIF 
     1109!               zshear(ji,jj) = 0.5 * zshear(ji,jj) 
     1110!               lshear(ji,jj) = .FALSE. 
     1111!             ENDIF 
    11061112         ELSE                ! zdb_bl test, note zshear set to zero 
    11071113           j_ddh(ji,jj) = 2 
     
    11281134         &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
    11291135        END IF 
    1130         zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
    1131                &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
    1132                &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
    1133                &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1136        zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * zwbav(ji,jj) & 
     1137               &                  - zalpha_s * zrf_shear * zustar(ji,jj)**3 / zhml(ji,jj) & 
     1138               &         + zr_stokes * ( zalpha_s * EXP( -1.5_wp * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1139               &                                         - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
    11341140          ! 
    11351141      ENDIF 
     
    11421148        IF ( lconv(ji,jj) ) THEN 
    11431149! Unstable OSBL 
    1144            zwb_shr = -za_wb_s * zshear(ji,jj) 
     1150           zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * zshear(ji,jj) 
    11451151           IF ( j_ddh(ji,jj) == 0 ) THEN 
    11461152 
    1147 ! Developing shear layer, additional shear production possible. 
    1148  
    1149              zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) /  zhbl(ji,jj), 0._wp ) 
    1150              zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 
    1151              zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 
    1152  
    1153              zwb_shr = -za_wb_s * zshear(ji,jj) 
     1153! ! Developing shear layer, additional shear production possible. 
     1154 
     1155!              zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) /  zhbl(ji,jj), 0._wp ) 
     1156!              zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) 
     1157!              zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 
     1158 
     1159!              zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) 
     1160!              zwb_shr = MAX( zwb_shr, -0.25 * zshear_u ) 
    11541161 
    11551162           ENDIF 
    11561163           zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
    1157            zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1164!           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
    11581165        ELSE    ! IF ( lconv ) THEN - ENDIF 
    11591166! Stable OSBL  - shear production not coded for first attempt. 
    11601167        ENDIF  ! lconv 
    1161       ELSE  ! lshear 
    1162         IF ( lconv(ji,jj) ) THEN 
     1168      END IF  ! lshear 
     1169      IF ( lconv(ji,jj) ) THEN 
    11631170! Unstable OSBL 
    1164            zwb_shr = -za_wb_s * zshear(ji,jj) 
    1165            zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 
    1166            zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
    1167         ENDIF  ! lconv 
    1168       ENDIF    ! lshear 
     1171         zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 
     1172      END IF  ! lconv 
    11691173     END_2D 
    11701174     IF( ln_timing ) CALL timing_stop('zdf_osm_os') 
     
    11961200           ztemp = zdu(ji,jj) 
    11971201           zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
    1198            zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1202           zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
    11991203        END_2D 
    12001204        IF( ln_timing ) CALL timing_stop('zdf_osm_vr') 
     
    12211225      REAL(wp), DIMENSION(jpi,jpj)  :: znd_param 
    12221226      REAL(wp)                      :: zbuoy, ztmp, zpe_mle_layer 
    1223       REAL(wp)                      :: zpe_mle_ref, zwb_ent, zdbdz_mle_int 
     1227      REAL(wp)                      :: zpe_mle_ref, zdbdz_mle_int 
    12241228 
    12251229      IF( ln_timing ) CALL timing_start('zdf_osm_osf') 
     
    12581262        DO_2D( 0, 0, 0, 0 ) 
    12591263          IF ( lconv(ji,jj) ) THEN 
    1260             zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 
    1261                &                  - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 
    1262                &         + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 
    1263                &         - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
    1264             IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 
     1264            IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN 
    12651265              IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 
    12661266! MLE layer growing 
     
    13351335           jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
    13361336           zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & 
    1337                 &   / e3t(ji,jj,ibld(ji,jj),Kmm) 
     1337                &   / e3t(ji,jj,jkb,Kmm) 
    13381338           zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & 
    1339                 &   / e3t(ji,jj,ibld(ji,jj),Kmm) 
     1339                &   / e3t(ji,jj,jkb,Kmm) 
    13401340           zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
    13411341        ELSE 
     
    13771377                  zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 
    13781378                  zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
    1379                   DO jk = 2, ibld(ji,jj)+ibld_ext 
     1379                  DO jk = 2, ibld(ji,jj) 
    13801380                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 
    13811381                     IF ( znd <= zzeta_m ) THEN 
     
    14261426   END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 
    14271427 
    1428    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     1428   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) 
    14291429     !!--------------------------------------------------------------------- 
    14301430     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     
    14361436     !!---------------------------------------------------------------------- 
    14371437 
    1438     REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt        ! Rate of change of hbl 
     1438    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt        ! Rate of change of hbl 
    14391439 
    14401440    INTEGER :: jj, ji 
    14411441    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 
    1442     REAL(wp) :: zvel_max!, zwb_min 
    1443     REAL(wp) :: zzeta_m = 0.3 
    1444     REAL(wp) :: zgamma_c = 2.0 
    1445     REAL(wp) :: zdhoh = 0.1 
    1446     REAL(wp) :: alpha_bc = 0.5 
    1447     REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 
     1442    REAL(wp) :: zvel_max, zddhdt 
     1443    REAL(wp), PARAMETER :: zzeta_m  = 0.3_wp 
     1444    REAL(wp), PARAMETER :: zgamma_c = 2.0_wp 
     1445    REAL(wp), PARAMETER :: zdhoh    = 0.1_wp 
     1446    REAL(wp), PARAMETER :: zalpha_b = 0.3_wp 
     1447    REAL(wp), PARAMETER :: a_ddh    = 2.5_wp, a_ddh_2 = 3.5 ! also in pycnocline_depth 
    14481448 
    14491449    IF( ln_timing ) CALL timing_start('zdf_osm_cd') 
     
    14661466                ! OSBL is deepening, entrainment > restratification 
    14671467                IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
    1468 ! *** Used for shear Needs to be changed to work stabily 
    1469 !                zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 
    1470 !                zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 
    1471 !                zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 
    1472 !                za_1 = 1.0 / zgamma_b**2 - 0.017 
    1473 !                za_2 = 1.0 / zgamma_b**3 - 0.0025 
    1474 !                zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 
    1475                    zpsi = 0._wp 
    1476                    zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    1477                    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) 
     1468                   zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1469                   zpsi = -1.0_wp * zalpha_pyc(ji,jj) * ( zwb0(ji,jj) - ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj) 
     1470                   zpsi = zpsi - 4.0_wp * ( 1.0_wp + zdh(ji,jj) / zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) 
     1471                   zpsi = zalpha_b * MAX( zpsi, 0.0_wp ) 
     1472                   zdhdt(ji,jj) = -1.0_wp * ( zwb_ent(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15_wp ) )   & 
     1473                      &           + zpsi / MAX( zdb_bl(ji,jj), 1e-15 ) 
    14781474                   IF ( j_ddh(ji,jj) == 1 ) THEN 
    14791475                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 
     
    14831479                     ENDIF 
    14841480! Relaxation to dh_ref = zari * hbl 
    1485                      zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
    1486  
    1487                    ELSE  ! j_ddh == 0 
     1481                     zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     1482 
     1483                  ELSE IF ( j_ddh(ji,jj) == 0 ) THEN 
    14881484! Growing shear layer 
    1489                      zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
    1490                    ENDIF ! j_ddh 
    1491                      zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 
     1485                     zddhdt = -1.0_wp * a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     1486                     zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1e-8_wp ) ) * zddhdt 
     1487                  ELSE 
     1488                     zddhdt = 0.0_wp 
     1489                  ENDIF ! j_ddh 
     1490                  zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0_wp * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) )   & 
     1491                     &                          * zddhdt / zdb_bl(ji,jj) 
    14921492                ELSE    ! zdb_bl >0 
    14931493                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     
    15161516                 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
    15171517             ELSE 
    1518                  zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     1518                 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    15191519             ENDIF 
    15201520             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
     
    15631563                 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 
    15641564             ELSE 
    1565                  zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
     1565                 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    15661566             ENDIF 
    15671567             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
     
    16951695      INTEGER :: jj, ji 
    16961696      INTEGER :: inhml 
    1697       REAL(wp) :: zari, ztau, zdh_ref 
    1698       REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 
     1697      REAL(wp) :: zari, ztau, zdh_ref, zddhdt 
     1698      REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 
    16991699 
    17001700      IF( ln_timing ) CALL timing_start('zdf_osm_pt') 
     
    17031703      IF ( lshear(ji,jj) ) THEN 
    17041704         IF ( lconv(ji,jj) ) THEN 
     1705          IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 
    17051706           IF ( j_ddh(ji,jj) == 0 ) THEN 
    17061707! ddhdt for pycnocline determined in osm_calculate_dhdt 
    1707              dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt 
     1708!             zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 
     1709!             zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 
     1710! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer 
     1711!             dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) ) 
     1712             ztau = MAX( zdb_bl(ji,jj) * ( 0.625_wp * hbl(ji,jj) ) / ( a_ddh * MAX( -1.0_wp * zwb_ent(ji,jj), 1e-12_wp ) ), 2.0_wp * rn_Dt ) 
     1713             dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.625_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 
     1714             IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8_wp * zhbl(ji,jj) 
    17081715           ELSE 
    17091716! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
     
    17171724             IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 
    17181725           ENDIF 
    1719  
     1726          ELSE 
     1727           ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5_wp * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt ) 
     1728           dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.2_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 
     1729           IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj) 
     1730          END IF 
    17201731         ELSE ! lconv 
    17211732! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 
     
    19101921 
    19111922 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
    1912   SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
     1923  SUBROUTINE zdf_osm_mle_parameters( pmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 
    19131924      !!---------------------------------------------------------------------- 
    19141925      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
     
    19211932      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
    19221933 
     1934      REAL(wp), DIMENSION(jpi,jpj)     :: pmld   ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
    19231935      INTEGER, DIMENSION(jpi,jpj)      :: mld_prof 
    19241936      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 
     
    19611973        hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) 
    19621974       IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 
     1975       ! For now try just set hmle to zmld 
     1976       hmle(ji,jj) = pmld(ji,jj) 
    19631977      END_2D 
    19641978 
     
    20812095   END SUBROUTINE zdf_osm_vertical_average 
    20822096 
    2083    SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, kbld_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear,   & 
     2097   SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear,             & 
    20842098      &                          pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla,       & 
    20852099      &                          pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml,                  & 
     
    20972111      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer 
    20982112      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kp_ext         ! Offset for external level 
    2099       INTEGER,                    INTENT(in   ) ::   kbld_ext       ! Offset for external level 
    21002113      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags 
    21012114      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
     
    22002213            IF ( jk <= kbld(ji,jj) ) THEN 
    22012214               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    2202                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   & 
     2215               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) *   & 
    22032216                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 
    2204                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   & 
     2217               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) *   & 
    22052218                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 
    22062219            END IF 
     
    22512264      ! ---------------------------------------------------------------------- 
    22522265      WHERE ( ldconv(A2D(0)) ) 
    2253          zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   & 
     2266         zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) /   & 
    22542267            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
    2255          zsc_ws_1(:,:)  = pwbav(A2D(0)) * pws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   & 
     2268         zsc_ws_1(:,:)  = pwbav(A2D(0)) * pws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) /   & 
    22562269            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
    22572270      ELSEWHERE 
     
    22632276            IF ( jk <= kmld(ji,jj) ) THEN 
    22642277               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    2265                ! Calculate turbulent length scale 
    2266                zl_c   = 0.9_wp * ( 1.0_wp - EXP( -7.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         & 
    2267                   &     ( 1.0_wp - EXP( -15.0_wp * ( 1.1_wp - zznd_ml ) ) ) 
    2268                zl_l   = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         & 
    2269                   &     ( 1.0_wp - EXP( -5.0_wp  * ( 1.0_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 
     2278               ! Calculate turbulent time scale 
     2279               zl_c   = 0.9_wp * ( 1.0_wp - EXP( -5.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) *                         & 
     2280                  &     ( 1.0_wp - EXP( -15.0_wp * ( 1.2_wp - zznd_ml ) ) ) 
     2281               zl_l   = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) *                         & 
     2282                  &     ( 1.0_wp - EXP( -8.0_wp  * ( 1.15_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 
    22702283               zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) 
    22712284               ! Non-gradient buoyancy terms 
    2272                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.5_wp * zsc_wth_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 
    2273                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.5_wp *  zsc_ws_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 
     2285               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_wth_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) 
     2286               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.4_wp *  zsc_ws_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) 
    22742287            END IF 
    22752288         ELSE   ! Stable conditions 
     
    22882301            zws_ent(ji,jj)      = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *   & 
    22892302               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 
    2290             zbuoy_pyc_sc        = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 
    2291             zdelta_pyc          = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird /   & 
    2292                &                  SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 
    2293             zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   & 
    2294                &                  zdelta_pyc**2 / pdh(ji,jj) 
    2295             zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   & 
    2296                &                  zdelta_pyc**2 / pdh(ji,jj) 
    2297             zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 
     2303            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 
     2304               zbuoy_pyc_sc        = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 
     2305               zdelta_pyc          = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird /   & 
     2306               &                       SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 
     2307               zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   & 
     2308               &                     zdelta_pyc**2 / pdh(ji,jj) 
     2309               zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   & 
     2310               &                     zdelta_pyc**2 / pdh(ji,jj) 
     2311               zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 
     2312            END IF 
    22982313         END IF 
    22992314      END_2D 
     
    23072322               &              0.045_wp * ( ( zws_ent(ji,jj)  * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
    23082323               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 
    2309             ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              & 
    2310                &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
    2311                &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
    2312             ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp  * zws_pyc_sc_1(ji,jj) *                              & 
    2313                &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
    2314                &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
     2324            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 
     2325               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              & 
     2326                  &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
     2327                  &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
     2328               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp  * zws_pyc_sc_1(ji,jj) *                              & 
     2329                  &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
     2330                  &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
     2331            END IF 
    23152332         END IF   ! End of pycnocline 
    23162333      END_3D 
     
    23482365      ! 
    23492366      DO_2D( 0, 0, 0, 0 ) 
    2350          IF ( ldpyc(ji,jj) ) THEN 
     2367         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 
    23512368            IF ( k_ddh(ji,jj) == 0 ) THEN 
    23522369               ! Place holding code. Parametrization needs checking for these conditions. 
    2353                zomega         = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
     2370               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
    23542371               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
    23552372               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
    23562373            ELSE 
    2357                zomega         = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
     2374               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
    23582375               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
    23592376               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
     
    23672384      END_2D 
    23682385      DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array. 
    2369          IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2386         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
    23702387            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    23712388            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) *                 & 
     
    23752392               &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   & 
    23762393               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
    2377          END IF   ! ldpyc 
     2394         END IF   ! ldconv .AND. ldpyc 
    23782395      END_3D 
    23792396      ! 
     
    26082625      ! 
    26092626      DO_2D( 0, 0, 0, 0 ) 
    2610          ghamt(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
    2611          ghams(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
    2612          ghamu(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
    2613          ghamv(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
     2627         ghamt(ji,jj,kbld(ji,jj)) = 0.0_wp 
     2628         ghams(ji,jj,kbld(ji,jj)) = 0.0_wp 
     2629         ghamu(ji,jj,kbld(ji,jj)) = 0.0_wp 
     2630         ghamv(ji,jj,kbld(ji,jj)) = 0.0_wp 
    26142631      END_2D 
    26152632      ! 
Note: See TracChangeset for help on using the changeset viewer.