Changeset 12218


Ignore:
Timestamp:
2019-12-12T17:01:31+01:00 (8 months ago)
Author:
agn
Message:

FKOSM mods for zdfosm.F90

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90

    r12216 r12218  
    7676   PUBLIC   tra_osm       ! routine called by step.F90 
    7777   PUBLIC   trc_osm       ! routine called by trcstp.F90 
    78    PUBLIC   dyn_osm       ! routine called by 'step.F90' 
     78   PUBLIC   dyn_osm       ! routine called by step.F90 
     79 
     80   PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90 
    7981 
    8082   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux 
     
    8486   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt 
    8587   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth 
    86    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbli     !: intial boundary layer depth for stable blayer 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline 
     89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth 
    8790   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift. 
     91 
     92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts 
     93   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML 
     96   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer. 
    8897 
    8998   !                      !!** Namelist  namzdf_osm  ** 
    9099   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la 
     100 
     101   LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 
     102 
    91103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    92104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
     
    102114   LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing 
    103115   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s) 
     116 
     117! OSMOSIS mixed layer eddy parametrization constants 
     118   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt 
     119   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient 
     120   !                                        ! parameters used in nn_osm_mle = 0 case 
     121   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front 
     122   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer 
     123   !                                        ! parameters used in nn_osm_mle = 1 case 
     124   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front 
     125   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK 
     126   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
     127   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 
     128   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 
     129   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base. 
     130   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE. 
     131 
    104132 
    105133   !                                    !!! ** General constants  ** 
     
    122150      !!                 ***  FUNCTION zdf_osm_alloc  *** 
    123151      !!---------------------------------------------------------------------- 
    124      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 
    125           &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     & 
    126           &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc ) 
     152     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 
     153          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     154          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     155 
     156     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 
     157          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 
     158 
     159!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ? 
     160!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 
     161!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 
     162!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
     163!     IF ( ln_osm_mle ) THEN 
     164!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 
     165!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 
     166!     ENDIF 
     167 
    127168     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 
    128169     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     
    169210      !! 
    170211      INTEGER ::   ji, jj, jk                   ! dummy loop indices 
     212 
     213      INTEGER ::   jl                   ! dummy loop indices 
     214 
    171215      INTEGER ::   ikbot, jkmax, jkm1, jkp2     ! 
    172216 
     
    200244      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average 
    201245      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux 
     246 
     247 
     248      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL 
     249      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux 
     250      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 
     251      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK 
     252 
    202253      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift 
    203254      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number 
     
    217268      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    218269      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid 
     270 
     271      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 
     272      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid 
     273 
    219274      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    220275      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
     276      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure. 
     277      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients 
     278 
     279      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
     280 
    221281      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer 
    222282      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer 
     
    289349      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
    290350      ! 
     351      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 
     352 
     353      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
     354         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 
     355         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 
     356         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 
     357         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp 
     358      ENDIF 
     359      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
     360 
    291361      ! Flux-Gradient arrays. 
    292362      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp 
     
    432502         END DO 
    433503      END DO 
    434  
     504      ! Averages over well-mixed and boundary layer 
     505      ! Alan: do we need zb_nl?, zb_ml? 
     506      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
     507      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     508! External gradient 
     509      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     510 
     511 
     512      IF ( ln_osm_mle ) THEN 
     513         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     514         CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     515       ENDIF 
     516 
     517! Rate of change of hbl 
     518      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
    435519      ! Calculate averages over depth of boundary layer 
     520         DO jj = 2, jpjm1 
     521            DO ji = 2, jpim1 
     522               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 
     523               ! adjustment to represent limiting by ocean bottom 
     524               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
     525            END DO 
     526         END DO 
     527 
    436528      imld = ibld           ! use imld to hold previous blayer index 
    437       ibld(:,:) = 3 
    438  
    439       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
    440       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    441       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     529      ibld(:,:) = 4 
    442530 
    443531      DO jk = 4, jpkm1 
     
    550638                   ENDIF 
    551639                ENDIF 
    552                 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    553                 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 
     640                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
     641                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 
     642                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 
    554643                ! could be taken out, take account of entrainment represents as a diffusivity 
    555644                ! should remove w from here, represents entrainment 
     
    9191008       END DO 
    9201009 
    921        IF(ln_dia_osm) THEN 
    922           IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
    923        END IF 
    924  
    9251010! KPP-style Ri# mixing 
    9261011       IF( ln_kpprimix) THEN 
     
    9741059          END DO 
    9751060       END IF ! ln_convmix = .true. 
     1061 
     1062 
     1063  
     1064       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing 
     1065          DO jj = 2 , jpjm1 
     1066             DO ji = 2, jpim1 
     1067                IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 
     1068            ! Calculate MLE flux profiles 
     1069                   ! DO jk = 1, mld_prof(ji,jj) 
     1070                   !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1071                   !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 
     1072                   !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1073                   !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 
     1074                   !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 
     1075                   ! END DO 
     1076            ! Calculate MLE flux contribution from surface fluxes 
     1077                   DO jk = 1, ibld(ji,jj) 
     1078                     znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 
     1079                     ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 
     1080                     ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 
     1081                    END DO 
     1082                    DO jk = 1, mld_prof(ji,jj) 
     1083                      znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 
     1084                      ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 
     1085                      ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 
     1086                    END DO 
     1087            ! Viscosity for MLEs 
     1088                    DO jk = ibld(ji,jj), mld_prof(ji,jj) 
     1089                      zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 
     1090                    END DO 
     1091            ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 
     1092                    jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 
     1093                    jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 
     1094                    DO jk = mld_prof(ji,jj),  jl 
     1095                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 
     1096                            & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 
     1097                            & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 
     1098                    END DO 
     1099                ENDIF 
     1100            END DO 
     1101          END DO 
     1102       ENDIF 
     1103 
     1104       IF(ln_dia_osm) THEN 
     1105          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
     1106          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
     1107          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
     1108       END IF 
     1109 
    9761110 
    9771111       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
     
    10581192         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    10591193         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    1060          IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine 
    1061          IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine 
    1062          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML 
     1194         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
     1195         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
     1196         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
     1197         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
     1198         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
     1199 
     1200         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )               ! FK layer depth 
     1201         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )               ! FK target layer depth 
     1202         IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux 
     1203         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML 
     1204         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 
     1205         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )            ! FK dtdx at u-pt 
     1206         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )            ! FK dtdy at v-pt 
     1207         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )            ! FK dtdx at u-pt 
     1208         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )            ! FK dsdy at v-pt 
     1209         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )            ! FK dbdx at u-pt 
     1210         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )            ! FK dbdy at v-pt 
     1211         IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1212         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 
     1213 
    10631214      END IF 
    1064       ! Lateral boundary conditions on p_avt  (sign unchanged) 
    1065       CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 
     1215 
     1216CONTAINS 
     1217 
     1218 
     1219   ! Alan: do we need zb? 
     1220   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 
     1221     !!--------------------------------------------------------------------- 
     1222     !!                   ***  ROUTINE zdf_vertical_average  *** 
     1223     !! 
     1224     !! ** Purpose : Determines vertical averages from surface to jnlev. 
     1225     !! 
     1226     !! ** Method  : Averages are calculated from the surface to jnlev. 
     1227     !!              The external level used to calculate differences is ibld+ibld_ext 
     1228     !! 
     1229     !!---------------------------------------------------------------------- 
     1230 
     1231        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over. 
     1232 
     1233        ! Alan: do we need zb? 
     1234        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity 
     1235        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components 
     1236        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 
     1237        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components. 
     1238 
     1239        INTEGER :: jk, ji, jj 
     1240        REAL(wp) :: zthick, zthermal, zbeta 
     1241 
     1242 
     1243        zt   = 0._wp 
     1244        zs   = 0._wp 
     1245        zu   = 0._wp 
     1246        zv   = 0._wp 
     1247        DO jj = 2, jpjm1                                 !  Vertical slab 
     1248         DO ji = 2, jpim1 
     1249            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1250            zbeta    = rab_n(ji,jj,1,jp_sal) 
     1251               ! average over depth of boundary layer 
     1252            zthick = epsln 
     1253            DO jk = 2, jnlev_av(ji,jj) 
     1254               zthick = zthick + e3t_n(ji,jj,jk) 
     1255               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 
     1256               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     1257               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) & 
     1258                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 
     1259                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 
     1260               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) & 
     1261                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 
     1262                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 
     1263            END DO 
     1264            zt(ji,jj) = zt(ji,jj) / zthick 
     1265            zs(ji,jj) = zs(ji,jj) / zthick 
     1266            zu(ji,jj) = zu(ji,jj) / zthick 
     1267            zv(ji,jj) = zv(ji,jj) / zthick 
     1268            ! Alan: do we need zb? 
     1269            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 
     1270            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 
     1271            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 
     1272                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 
     1273            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 
     1274                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 
     1275            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 
     1276         END DO 
     1277        END DO 
     1278   END SUBROUTINE zdf_osm_vertical_average 
     1279 
     1280   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 
     1281     !!--------------------------------------------------------------------- 
     1282     !!                   ***  ROUTINE zdf_velocity_rotation  *** 
     1283     !! 
     1284     !! ** Purpose : Rotates frame of reference of averaged velocity components. 
     1285     !! 
     1286     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w 
     1287     !! 
     1288     !!---------------------------------------------------------------------- 
     1289 
     1290        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle 
     1291        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current 
     1292        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline 
     1293 
     1294        INTEGER :: ji, jj 
     1295        REAL(wp) :: ztemp 
     1296 
     1297        DO jj = 2, jpjm1 
     1298           DO ji = 2, jpim1 
     1299              ztemp = zu(ji,jj) 
     1300              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 
     1301              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1302              ztemp = zdu(ji,jj) 
     1303              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 
     1304              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 
     1305           END DO 
     1306        END DO 
     1307    END SUBROUTINE zdf_osm_velocity_rotation 
     1308 
     1309    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 
     1310     !!--------------------------------------------------------------------- 
     1311     !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     1312     !! 
     1313     !! ** Purpose : Calculates the gradients below the OSBL 
     1314     !! 
     1315     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     1316     !! 
     1317     !!---------------------------------------------------------------------- 
     1318 
     1319     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
     1320 
     1321     INTEGER :: jj, ji, jkb, jkb1 
     1322     REAL(wp) :: zthermal, zbeta 
     1323 
     1324 
     1325     DO jj = 2, jpjm1 
     1326        DO ji = 2, jpim1 
     1327           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1328              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     1329              zbeta    = rab_n(ji,jj,1,jp_sal) 
     1330              jkb = ibld(ji,jj) + ibld_ext 
     1331              jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
     1332              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 
     1333                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1334              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 
     1335                   &   / e3t_n(ji,jj,ibld(ji,jj)) 
     1336              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
     1337           END IF 
     1338        END DO 
     1339     END DO 
     1340    END SUBROUTINE zdf_osm_external_gradients 
     1341 
     1342    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 
     1343 
     1344     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
     1345 
     1346     INTEGER :: jk, jj, ji 
     1347     REAL(wp) :: ztgrad, zsgrad, zbgrad 
     1348     REAL(wp) :: zgamma_b_nd, zgamma_c, znd 
     1349     REAL(wp) :: zzeta_s=0.3 
     1350 
     1351     DO jj = 2, jpjm1 
     1352        DO ji = 2, jpim1 
     1353           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1354              IF ( lconv(ji,jj) ) THEN  ! convective conditions 
     1355                    IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 
     1356                         & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 
     1357                         &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0 
     1358                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj) 
     1359                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj) 
     1360                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj) 
     1361                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1362                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 
     1363                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 
     1364                       DO jk = 2, ibld(ji,jj)+ibld_ext 
     1365                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
     1366                          IF ( znd <= zzeta_s ) THEN 
     1367                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * & 
     1368                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1369                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * & 
     1370                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1371                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * & 
     1372                                  & EXP( -6.0 * ( znd -zzeta_s )**2 ) 
     1373                          ELSE 
     1374                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1375                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1376                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 
     1377                          ENDIF 
     1378                       END DO 
     1379                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 
     1380              ELSE 
     1381                 ! stable conditions 
     1382                 ! if pycnocline profile only defined when depth steady of increasing. 
     1383                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     1384                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1385                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     1386                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
     1387                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
     1388                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
     1389                          DO jk = 2, ibld(ji,jj) 
     1390                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1391                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1392                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1393                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     1394                          END DO 
     1395                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1396                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
     1397                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
     1398                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
     1399                          DO jk = 2, ibld(ji,jj) 
     1400                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     1401                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1402                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1403                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     1404                          END DO 
     1405                       ENDIF ! IF (zhol >=0.5) 
     1406                    ENDIF    ! IF (zdb_bl> 0.) 
     1407                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1408              ENDIF          ! IF (lconv) 
     1409           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1410        END DO 
     1411     END DO 
     1412 
     1413    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
     1414 
     1415    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
     1416      !!--------------------------------------------------------------------- 
     1417      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
     1418      !! 
     1419      !! ** Purpose : Calculates velocity shear in the pycnocline 
     1420      !! 
     1421      !! ** Method  : 
     1422      !! 
     1423      !!---------------------------------------------------------------------- 
     1424 
     1425      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
     1426 
     1427      INTEGER :: jk, jj, ji 
     1428      REAL(wp) :: zugrad, zvgrad, znd 
     1429      REAL(wp) :: zzeta_v = 0.45 
    10661430      ! 
    1067    END SUBROUTINE zdf_osm 
     1431      DO jj = 2, jpjm1 
     1432         DO ji = 2, jpim1 
     1433            ! 
     1434            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     1435               IF ( lconv (ji,jj) ) THEN 
     1436                  ! Unstable conditions 
     1437                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     1438                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     1439                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     1440                  !Alan is this right? 
     1441                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     1442                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     1443                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     1444                       &      )/ (zdh(ji,jj)  + epsln ) 
     1445                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     1446                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     1447                     IF ( znd <= 0.0 ) THEN 
     1448                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     1449                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     1450                     ELSE 
     1451                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     1452                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     1453                     ENDIF 
     1454                  END DO 
     1455               ELSE 
     1456                  ! stable conditions 
     1457                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     1458                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     1459                  DO jk = 2, ibld(ji,jj) 
     1460                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1461                     IF ( znd < 1.0 ) THEN 
     1462                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     1463                     ELSE 
     1464                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     1465                     ENDIF 
     1466                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     1467                  END DO 
     1468               ENDIF 
     1469               ! 
     1470            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1471         END DO 
     1472      END DO 
     1473    END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1474 
     1475    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 
     1476     !!--------------------------------------------------------------------- 
     1477     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     1478     !! 
     1479     !! ** Purpose : Calculates the rate at which hbl changes. 
     1480     !! 
     1481     !! ** Method  : 
     1482     !! 
     1483     !!---------------------------------------------------------------------- 
     1484 
     1485    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl 
     1486 
     1487    INTEGER :: jj, ji 
     1488    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 
     1489    REAL(wp) :: zvel_max, zwb_min 
     1490    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     1491    REAL(wp) :: zzeta_m = 0.3 
     1492    REAL(wp) :: zgamma_c = 2.0 
     1493    REAL(wp) :: zdhoh = 0.1 
     1494    REAL(wp) :: alpha_bc = 0.5 
     1495 
     1496    DO jj = 2, jpjm1 
     1497       DO ji = 2, jpim1 
     1498          IF ( lconv(ji,jj) ) THEN    ! Convective 
     1499             ! Alan is this right?  Yes, it's a bit different from the previous relationship 
     1500             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 
     1501             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1502             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 
     1503             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 
     1504             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
     1505             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
     1506             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1507                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1508 
     1509             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     1510                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 
     1511                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 
     1512                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 
     1513             ! 
     1514             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 
     1515 
     1516             IF ( ln_osm_mle ) THEN 
     1517                  !  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 + & 
     1518                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL 
     1519                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
     1520                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
     1521                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     1522                ELSE 
     1523                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     1524                ENDIF 
     1525                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1526                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
     1527                   ! OSBL is deepening, entrainment > restratification 
     1528                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1529                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1530                   ELSE 
     1531                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
     1532                   ENDIF 
     1533                ELSE 
     1534                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     1535                   zdhdt(ji,jj) = - zvel_mle(ji,jj) 
     1536 
     1537 
     1538                ENDIF 
     1539 
     1540             ELSE 
     1541                ! Fox-Kemper not used. 
     1542 
     1543                  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) / & 
     1544                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1545                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
     1546                ! added ajgn 23 July as temporay fix 
     1547 
     1548             ENDIF 
     1549 
     1550             zdhdt_2(ji,jj) = 0._wp 
     1551 
     1552                ! commented out ajgn 23 July as temporay fix 
     1553!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 
     1554! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 
     1555!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1556!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 
     1557!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 
     1558!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 
     1559!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 
     1560!                     ! Alan no idea what this should be? 
     1561!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 
     1562!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 
     1563!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 
     1564!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 
     1565!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 
     1566!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 
     1567!                 ENDIF 
     1568            ELSE                        ! Stable 
     1569                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 
     1570                zdhdt_2(ji,jj) = 0._wp 
     1571                IF ( zdhdt(ji,jj) < 0._wp ) THEN 
     1572                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     1573                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 
     1574                ELSE 
     1575                    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) ) 
     1576                ENDIF 
     1577                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 
     1578            ENDIF 
     1579         END DO 
     1580      END DO 
     1581    END SUBROUTINE zdf_osm_calculate_dhdt 
     1582 
     1583    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     1584     !!--------------------------------------------------------------------- 
     1585     !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
     1586     !! 
     1587     !! ** Purpose : Increments hbl. 
     1588     !! 
     1589     !! ** Method  : If thechange in hbl exceeds one model level the change is 
     1590     !!              is calculated by moving down the grid, changing the buoyancy 
     1591     !!              jump. This is to ensure that the change in hbl does not 
     1592     !!              overshoot a stable layer. 
     1593     !! 
     1594     !!---------------------------------------------------------------------- 
     1595 
     1596 
     1597    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl. 
     1598 
     1599    INTEGER :: jk, jj, ji, jm 
     1600    REAL(wp) :: zhbl_s, zvel_max, zdb 
     1601    REAL(wp) :: zthermal, zbeta 
     1602 
     1603     DO jj = 2, jpjm1 
     1604         DO ji = 2, jpim1 
     1605           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     1606! 
     1607! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     1608! 
     1609              zhbl_s = hbl(ji,jj) 
     1610              jm = imld(ji,jj) 
     1611              zthermal = rab_n(ji,jj,1,jp_tem) 
     1612              zbeta = rab_n(ji,jj,1,jp_sal) 
     1613 
     1614 
     1615              IF ( lconv(ji,jj) ) THEN 
     1616!unstable 
     1617 
     1618                 IF( ln_osm_mle ) THEN 
     1619                    zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     1620                 ELSE 
     1621 
     1622                    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) / & 
     1623                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1624 
     1625                 ENDIF 
     1626 
     1627                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1628                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
     1629                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 
     1630                         &  0.0 ) + zvel_max 
     1631 
     1632 
     1633                    IF ( ln_osm_mle ) THEN 
     1634                       zhbl_s = zhbl_s + MIN( & 
     1635                        & 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) ), & 
     1636                        & e3w_n(ji,jj,jm) ) 
     1637                    ELSE 
     1638                      zhbl_s = zhbl_s + MIN( & 
     1639                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
     1640                        & e3w_n(ji,jj,jm) ) 
     1641                    ENDIF 
     1642 
     1643                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1644 
     1645                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
     1646                 END DO 
     1647                 hbl(ji,jj) = zhbl_s 
     1648                 ibld(ji,jj) = jm 
     1649             ELSE 
     1650! stable 
     1651                 DO jk = imld(ji,jj), ibld(ji,jj) 
     1652                    zdb = MAX( & 
     1653                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 
     1654                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 
     1655                         & 0.0 ) + & 
     1656             &       2.0 * zvstr(ji,jj)**2 / zhbl_s 
     1657 
     1658                    ! Alan is thuis right? I have simply changed hbli to hbl 
     1659                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 
     1660                    zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & 
     1661               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 
     1662                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 
     1663                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 
     1664 
     1665                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 
     1666                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
     1667                 END DO 
     1668             ENDIF   ! IF ( lconv ) 
     1669             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 
     1670             ibld(ji,jj) = MAX(jm, 4 ) 
     1671           ELSE 
     1672! change zero or one model level. 
     1673             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 
     1674           ENDIF 
     1675           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     1676         END DO 
     1677      END DO 
     1678 
     1679    END SUBROUTINE zdf_osm_timestep_hbl 
     1680 
     1681    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     1682      !!--------------------------------------------------------------------- 
     1683      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     1684      !! 
     1685      !! ** Purpose : Calculates thickness of the pycnocline 
     1686      !! 
     1687      !! ** Method  : The thickness is calculated from a prognostic equation 
     1688      !!              that relaxes the pycnocine thickness to a diagnostic 
     1689      !!              value. The time change is calculated assuming the 
     1690      !!              thickness relaxes exponentially. This is done to deal 
     1691      !!              with large timesteps. 
     1692      !! 
     1693      !!---------------------------------------------------------------------- 
     1694 
     1695      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness. 
     1696      ! 
     1697      INTEGER :: jj, ji 
     1698      INTEGER :: inhml 
     1699      REAL(wp) :: zari, ztau, zddhdt 
     1700 
     1701 
     1702      DO jj = 2, jpjm1 
     1703         DO ji = 2, jpim1 
     1704            IF ( lconv(ji,jj) ) THEN 
     1705 
     1706               IF( ln_osm_mle ) THEN 
     1707                  IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 
     1708                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 
     1709                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     1710                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1711                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1712                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1713                        ELSE                                                     ! unstable 
     1714                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1715                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1716                        ENDIF 
     1717                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1718                        zddhdt = zari * hbl(ji,jj) 
     1719                     ELSE 
     1720                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1721                        zddhdt = 0.2 * hbl(ji,jj) 
     1722                     ENDIF 
     1723                  ELSE 
     1724                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1725                     zddhdt = 0.2 * hbl(ji,jj) 
     1726                  ENDIF 
     1727               ELSE ! ln_osm_mle 
     1728                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 
     1729                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     1730                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1731                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1732                     ELSE                                                     ! unstable 
     1733                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 
     1734                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 
     1735                     ENDIF 
     1736                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 
     1737                     zddhdt = zari * hbl(ji,jj) 
     1738                  ELSE 
     1739                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1740                     zddhdt = 0.2 * hbl(ji,jj) 
     1741                  ENDIF 
     1742 
     1743               END IF  ! ln_osm_mle 
     1744 
     1745               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1746               ! Alan: this hml is never defined or used 
     1747            ELSE   ! IF (lconv) 
     1748               ztau = hbl(ji,jj) / zvstr(ji,jj) 
     1749               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     1750                  ! boundary layer deepening 
     1751                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     1752                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     1753                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     1754                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     1755                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 
     1756                  ELSE 
     1757                     zddhdt = 0.2 * hbl(ji,jj) 
     1758                  ENDIF 
     1759               ELSE     ! IF(dhdt < 0) 
     1760                  zddhdt = 0.2 * hbl(ji,jj) 
     1761               ENDIF    ! IF (dhdt >= 0) 
     1762               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 
     1763               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 
     1764               ! Alan: this hml is never defined or used -- do we need it? 
     1765            ENDIF       ! IF (lconv) 
     1766 
     1767            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 
     1768            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     1769            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 
     1770            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     1771            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     1772         END DO 
     1773      END DO 
     1774 
     1775    END SUBROUTINE zdf_osm_pycnocline_thickness 
     1776 
     1777 
     1778   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 
     1779      !!---------------------------------------------------------------------- 
     1780      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  *** 
     1781      !! 
     1782      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 
     1783      !! 
     1784      !! ** Method  : 
     1785      !! 
     1786      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1787      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1788 
     1789 
     1790      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 
     1791      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == ! 
     1792      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy 
     1793 
     1794      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1795      INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     1796      REAL(wp)                         :: zc 
     1797      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value 
     1798      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH 
     1799      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 
     1800      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv 
     1801!!---------------------------------------------------------------------- 
     1802      ! 
     1803      !                                      !==  MLD used for MLE  ==! 
     1804 
     1805      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point 
     1806      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
     1807      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     1808      DO jk = nlb10, jpkm1 
     1809         DO jj = 1, jpj                ! Mixed layer level: w-level  
     1810            DO ji = 1, jpi 
     1811               ikt = mbkt(ji,jj) 
     1812               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 
     1813               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     1814            END DO 
     1815         END DO 
     1816      END DO 
     1817      DO jj = 1, jpj 
     1818         DO ji = 1, jpi 
     1819            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 
     1820            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1821         END DO 
     1822      END DO 
     1823      ! ensure mld_prof .ge. ibld 
     1824      ! 
     1825      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     1826      ! 
     1827      ztm(:,:) = 0._wp 
     1828      zsm(:,:) = 0._wp 
     1829      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
     1830         DO jj = 1, jpj 
     1831            DO ji = 1, jpi 
     1832               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     1833               ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 
     1834               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 
     1835            END DO 
     1836         END DO 
     1837      END DO 
     1838      ! average temperature and salinity. 
     1839      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1840      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     1841      ! calculate horizontal gradients at u & v points 
     1842 
     1843      DO jj = 2, jpjm1 
     1844         DO ji = 1, jpim1 
     1845            zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1846            zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     1847            zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 
     1848            ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 
     1849            ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 
     1850         END DO 
     1851      END DO 
     1852 
     1853      DO jj = 1, jpjm1 
     1854         DO ji = 2, jpim1 
     1855            zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1856            zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     1857            zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 
     1858            ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 
     1859            ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 
     1860         END DO 
     1861      END DO 
     1862 
     1863      CALL eos_rab(ztsm_midu, zmld_midu, zabu) 
     1864      CALL eos_rab(ztsm_midv, zmld_midv, zabv) 
     1865 
     1866      DO jj = 2, jpjm1 
     1867         DO ji = 1, jpim1 
     1868            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 
     1869         END DO 
     1870      END DO 
     1871      DO jj = 1, jpjm1 
     1872         DO ji = 2, jpim1 
     1873            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 
     1874         END DO 
     1875      END DO 
     1876 
     1877 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 
     1878  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 
     1879      !!---------------------------------------------------------------------- 
     1880      !!                  ***  ROUTINE zdf_osm_mle_parameters  *** 
     1881      !! 
     1882      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 
     1883      !! 
     1884      !! ** Method  : 
     1885      !! 
     1886      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     1887      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     1888 
     1889      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle 
     1890      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     1891      INTEGER  ::   ii, ij, ik  ! local integers 
     1892      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     1893      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle 
     1894 
     1895      mld_prof(:,:) = 4 
     1896      DO jk = 5, jpkm1 
     1897        DO jj = 2, jpjm1 
     1898          DO ji = 2, jpim1 
     1899            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1900          END DO 
     1901        END DO 
     1902      END DO 
     1903      ! DO jj = 2, jpjm1 
     1904      !    DO ji = 1, jpim1 
     1905      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 
     1906      !    END DO 
     1907      !  END DO 
     1908   ! Timestep mixed layer eddy depth. 
     1909      DO jj = 2, jpjm1 
     1910        DO ji = 2, jpim1 
     1911          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 
     1912          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 
     1913             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. 
     1914          ELSE 
     1915            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 
     1916             IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 
     1917               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 
     1918             ELSE 
     1919               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 
     1920             ENDIF 
     1921          ENDIF 
     1922          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 
     1923        END DO 
     1924      END DO 
     1925 
     1926      mld_prof = 4 
     1927      DO jk = 5, jpkm1 
     1928        DO jj = 2, jpjm1 
     1929          DO ji = 2, jpim1 
     1930            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
     1931          END DO 
     1932        END DO 
     1933      END DO 
     1934      DO jj = 2, jpjm1 
     1935         DO ji = 2, jpim1 
     1936            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 
     1937         END DO 
     1938       END DO 
     1939   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 
     1940 
     1941      DO jj = 2, jpjm1 
     1942        DO ji = 2, jpim1 
     1943          IF ( lconv(ji,jj) ) THEN 
     1944             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
     1945             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 
     1946             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 
     1947             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 
     1948             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 
     1949             zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 
     1950                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 
     1951             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 
     1952      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 
     1953             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1)  
     1954             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf 
     1955          ENDIF 
     1956        END DO 
     1957      END DO 
     1958END SUBROUTINE zdf_osm_mle_parameters 
     1959 
     1960END SUBROUTINE zdf_osm 
    10681961 
    10691962 
     
    10861979     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    10871980          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
    1088           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 
     1981          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
     1982          & ,ln_osm_mle 
     1983! Namelist for Fox-Kemper parametrization. 
     1984      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
     1985           & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 
     1986 
    10891987     !!---------------------------------------------------------------------- 
    10901988     ! 
     
    11022000        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 
    11032001        WRITE(numout,*) '~~~~~~~~~~~~' 
    1104         WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters' 
    1105         WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la 
     2002        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters' 
     2003        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la 
     2004        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    11062005        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
    11072006        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
     
    11282027     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    11292028 
    1130      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     2029 
     2030     IF( ln_osm_mle ) THEN 
     2031! Initialise Fox-Kemper parametrization 
     2032         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme 
     2033         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 
     2034903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 
     2035 
     2036         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 
     2037         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 
     2038904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 
     2039         IF(lwm) WRITE ( numond, namosm_mle ) 
     2040 
     2041         IF(lwp) THEN                     ! Namelist print 
     2042            WRITE(numout,*) 
     2043            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 
     2044            WRITE(numout,*) '~~~~~~~~~~~~~' 
     2045            WRITE(numout,*) '   Namelist namosm_mle : ' 
     2046            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle 
     2047            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce 
     2048            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm' 
     2049            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's' 
     2050            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg' 
     2051            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c 
     2052            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s' 
     2053            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's' 
     2054         ENDIF         ! 
     2055     ENDIF 
     2056      ! 
     2057      IF(lwp) THEN 
     2058         WRITE(numout,*) 
     2059         IF( ln_osm_mle ) THEN 
     2060            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 
     2061            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation' 
     2062            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation' 
     2063         ELSE 
     2064            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 
     2065         ENDIF 
     2066      ENDIF 
     2067      ! 
     2068      IF( ln_osm_mle ) THEN                ! MLE initialisation 
     2069         ! 
     2070         rb_c = grav * rn_osm_mle_rho_c /rau0        ! Mixed Layer buoyancy criteria 
     2071         IF(lwp) WRITE(numout,*) 
     2072         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
     2073         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 
     2074         ! 
     2075         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            ! 
     2076! 
     2077         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation 
     2078            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  ) 
     2079            ! 
     2080         ENDIF 
     2081         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
     2082         z1_t2 = 2.e-5 
     2083         do jj=1,jpj 
     2084            do ji = 1,jpi 
     2085               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 
     2086            end do 
     2087         end do 
     2088         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 
     2089         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  ) 
     2090         ! 
     2091      ENDIF 
     2092 
     2093     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 
     2094 
    11312095 
    11322096     IF( ln_zdfddm) THEN 
     
    13092273        END DO 
    13102274     END DO 
    1311      hbl = MAX(hbl,epsln) 
    1312      hbli(:,:) = hbl(:,:) 
    1313      DEALLOCATE( imld_rst ) 
     2275 
     2276     IF( ln_osm_mle ) hmle(:,:) = hbl(:,:)            ! Initialise MLE depth. 
     2277 
    13142278     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
    13152279     wn(:,:,:) = 0._wp 
Note: See TracChangeset for help on using the changeset viewer.