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

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

Changeset 14803 for NEMO


Ignore:
Timestamp:
2021-05-06T21:01:51+02:00 (3 years ago)
Author:
smueller
Message:

Upgrade of ten local arrays to module arrays (ticket #2353)

File:
1 edited

Legend:

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

    r14802 r14803  
    163163   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_bl   ! Velocity average (u) 
    164164   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_bl   ! Velocity average (v) 
    165    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_bl   ! Buoyancy 
     165   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_bl   ! Buoyancy average 
     166   ! 
     167   ! Difference between layer average and parameter at the base of the layer: BL 
     168   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dt_bl   ! Temperature difference 
     169   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_ds_bl   ! Salinity difference 
     170   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_du_bl   ! Velocity difference (u) 
     171   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dv_bl   ! Velocity difference (v) 
     172   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_db_bl   ! Buoyancy difference 
    166173   ! 
    167174   ! Layer averages: ML 
     
    170177   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_ml   ! Velocity average (u) 
    171178   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_ml   ! Velocity average (v) 
    172    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_ml   ! Buoyancy 
     179   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_ml   ! Buoyancy average 
     180   ! 
     181   ! Difference between layer average and parameter at the base of the layer: ML 
     182   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dt_ml   ! Temperature difference 
     183   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_ds_ml   ! Salinity difference 
     184   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_du_ml   ! Velocity difference (u) 
     185   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dv_ml   ! Velocity difference (v) 
     186   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_db_ml   ! Buoyancy difference 
    173187   ! 
    174188   ! Layer averages: MLE 
     
    177191   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_mle  ! Velocity average (u) 
    178192   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_mle  ! Velocity average (v) 
    179    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_mle  ! Buoyancy 
     193   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_mle  ! Buoyancy average 
    180194   ! 
    181195   !            ** Namelist  namzdf_osm  ** 
     
    270284      zdf_osm_alloc = zdf_osm_alloc + ierr 
    271285      ! 
     286      ALLOCATE( av_dt_bl(jpi,jpj), av_ds_bl(jpi,jpj), av_du_bl(jpi,jpj), av_dv_bl(jpi,jpj), av_db_bl(jpi,jpj), STAT=ierr) 
     287      zdf_osm_alloc = zdf_osm_alloc + ierr 
     288      ! 
    272289      ALLOCATE( av_t_ml(jpi,jpj), av_s_ml(jpi,jpj), av_u_ml(jpi,jpj), av_v_ml(jpi,jpj), av_b_ml(jpi,jpj), STAT=ierr) 
     290      zdf_osm_alloc = zdf_osm_alloc + ierr 
     291      ! 
     292      ALLOCATE( av_dt_ml(jpi,jpj), av_ds_ml(jpi,jpj), av_du_ml(jpi,jpj), av_dv_ml(jpi,jpj), av_db_ml(jpi,jpj), STAT=ierr) 
    273293      zdf_osm_alloc = zdf_osm_alloc + ierr 
    274294      ! 
     
    354374      REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization. 
    355375      ! 
    356       REAL(wp), DIMENSION(jpi,jpj)     ::   zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl   ! Difference between blayer average and parameter at base of blayer 
    357       REAL(wp), DIMENSION(jpi,jpj)     ::   zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml   ! Difference between mixed layer average and parameter at base of blayer 
    358376      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline 
    359377      REAL(wp), DIMENSION(jpi,jpj)     ::   zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
     
    400418      av_b_ml(:,:)  = pp_large ; av_t_mle(:,:) = pp_large ; av_s_mle(:,:) = pp_large 
    401419      av_u_mle(:,:) = pp_large ; av_v_mle(:,:) = pp_large ; av_b_mle(:,:) = pp_large 
    402       zdt_bl(:,:)   = pp_large ; zds_bl(:,:)   = pp_large ; zdu_bl(:,:)  = pp_large 
    403       zdv_bl(:,:)   = pp_large ; zdb_bl(:,:)   = pp_large ; zdt_ml(:,:)  = pp_large 
    404       zds_ml(:,:)   = pp_large ; zdu_ml(:,:)   = pp_large ; zdv_ml(:,:)  = pp_large 
    405       zdb_ml(:,:)  = pp_large 
     420      av_dt_bl(:,:) = pp_large ; av_ds_bl(:,:) = pp_large ; av_du_bl(:,:) = pp_large 
     421      av_dv_bl(:,:) = pp_large ; av_db_bl(:,:) = pp_large ; av_dt_ml(:,:) = pp_large 
     422      av_ds_ml(:,:) = pp_large ; av_du_ml(:,:) = pp_large ; av_dv_ml(:,:) = pp_large 
     423      av_db_ml(:,:) = pp_large 
    406424      ! 
    407425      zdbdz_pyc(:,:,:)    = pp_large 
     
    677695      ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 
    678696      jp_ext(:,:) = 1   ! ag 19/03 
    679       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,   & 
    680          &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl,  & 
    681          &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     697      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,             & 
     698         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,  & 
     699         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    682700      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    683       CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    & 
    684          &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml,   & 
    685          &                           zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
     701      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
     702         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
     703         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 
    686704#ifdef key_osm_debug 
    687705      IF(narea==nn_narea_db) THEN 
     
    690708            & 'After averaging, with old hbl (& jp_ext==2), hml: zt_bl=', av_t_bl(ji,jj),& 
    691709            & ' zs_bl=', av_s_bl(ji,jj),  ' zb_bl=', av_b_bl(ji,jj),& 
    692             & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj),  ' zdb_bl=', zdb_bl(ji,jj),& 
     710            & 'zdt_bl=', av_dt_bl(ji,jj), ' zds_bl=', av_ds_bl(ji,jj),  ' zdb_bl=', av_db_bl(ji,jj),& 
    693711            & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj),  ' zb_ml=', av_b_ml(ji,jj),& 
    694             & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj),  ' zdb_ml=', zdb_ml(ji,jj),& 
    695             & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 
    696             & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 
     712            & 'zdt_ml=', av_dt_ml(ji,jj), ' zds_ml=', av_ds_ml(ji,jj),  ' zdb_ml=', av_db_ml(ji,jj),& 
     713            & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', av_du_bl(ji,jj), ' zdv_bl=', av_dv_bl(ji,jj),& 
     714            & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', av_du_ml(ji,jj), ' zdv_ml=', av_dv_ml(ji,jj) 
    697715         FLUSH(narea+100) 
    698716      END IF 
     
    700718      ! Velocity components in frame aligned with surface stress 
    701719      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml  ) 
    702       CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 
     720      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml ) 
    703721      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl  ) 
    704       CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 
     722      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl ) 
    705723#ifdef key_osm_debug 
    706724      IF(narea==nn_narea_db) THEN 
     
    708726         WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & 
    709727            & 'After rotation, with old hbl (& jp_ext==2), hml:', & 
    710             & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 
    711             & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 
     728            & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', av_du_bl(ji,jj), ' zdv_bl=', av_dv_bl(ji,jj),& 
     729            & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', av_du_ml(ji,jj), ' zdv_ml=', av_dv_ml(ji,jj) 
    712730         FLUSH(narea+100) 
    713731      END IF 
     
    715733      ! 
    716734      ! Determine the state of the OSBL, stable/unstable, shear/no shear 
    717       CALL zdf_osm_osbl_state( Kmm,    zwb_ent, zwb_min, zshear, zhbl,     & 
    718          &                     zhml,   zdh,     zdb_bl,  zdu_bl, zdv_bl,   & 
    719          &                     zdb_ml, zdu_ml,  zdv_ml ) 
     735      CALL zdf_osm_osbl_state( Kmm, zwb_ent, zwb_min, zshear, zhbl,     & 
     736         &                     zhml, zdh ) 
    720737      ! 
    721738#ifdef key_osm_debug 
     
    756773         ! Calculate max vertical FK flux zwb_fk & set logical descriptors 
    757774         CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zhbl, zhmle, zwb_ent,   & 
    758             &                        zdb_bl, zdbds_mle ) 
     775            &                        zdbds_mle ) 
    759776         ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 
    760777         CALL zdf_osm_mle_parameters( Kmm, mld_prof, zmld, zhmle, zvel_mle,   & 
     
    780797         l_mle(:,:)  = .FALSE. 
    781798         DO_2D( 0, 0, 0, 0 ) 
    782             IF ( l_conv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
     799            IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
    783800         END_2D 
    784801      ENDIF   ! ln_osm_mle 
     
    816833      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 
    817834      jp_ext(:,:) = 1   ! ag 19/03 
    818       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,    & 
    819          &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl,   & 
    820          &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
     835      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,              & 
     836         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   & 
     837         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    821838      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    822       CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    & 
    823          &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml,   & 
    824          &                           zds_ml, zdb_ml, zdu_ml, zdv_ml )   ! ag 19/03 
     839      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
     840         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
     841         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )   ! ag 19/03 
    825842#ifdef key_osm_debug 
    826843      IF(narea==nn_narea_db) THEN 
     
    829846            & 'After averaging, with old hbl (&correct jp_ext), hml: zt_bl=', av_t_bl(ji,jj),& 
    830847            & ' zs_bl=', av_s_bl(ji,jj),  ' zb_bl=', av_b_bl(ji,jj),& 
    831             & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj),  ' zdb_bl=', zdb_bl(ji,jj),& 
     848            & 'zdt_bl=', av_dt_bl(ji,jj), ' zds_bl=', av_ds_bl(ji,jj),  ' zdb_bl=', av_db_bl(ji,jj),& 
    832849            & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj),  ' zb_ml=', av_b_ml(ji,jj),& 
    833             & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj),  ' zdb_ml=', zdb_ml(ji,jj),& 
    834             & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 
    835             & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 
     850            & 'zdt_ml=', av_dt_ml(ji,jj), ' zds_ml=', av_ds_ml(ji,jj),  ' zdb_ml=', av_db_ml(ji,jj),& 
     851            & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', av_du_bl(ji,jj), ' zdv_bl=', av_dv_bl(ji,jj),& 
     852            & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', av_du_ml(ji,jj), ' zdv_ml=', av_dv_ml(ji,jj) 
    836853         FLUSH(narea+100) 
    837854      END IF 
     
    839856      ! 
    840857      ! Rate of change of hbl 
    841       CALL zdf_osm_calculate_dhdt( zdhdt, zhbl, zdh, zwb_ent, zwb_min,               & 
    842          &                         zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk,   & 
    843          &                         zvel_mle ) 
     858      CALL zdf_osm_calculate_dhdt( zdhdt, zhbl, zdh, zwb_ent, zwb_min,   & 
     859         &                         zdbdz_bl_ext, zwb_fk_b, zwb_fk, zvel_mle ) 
    844860      ! Test if surface boundary layer coupled to bottom 
    845861      l_coup(:,:) = .FALSE.   ! ag 19/03 
     
    882898      ! Recalculate BL averages and differences using new BL depth 
    883899      jp_ext(:,:) = 1   ! ag 19/03 
    884       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,    & 
    885          &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl,   & 
    886          &                           zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    887       ! 
    888       CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zdb_bl,   & 
    889          &                               zhbl, zwb_ent, zdbdz_bl_ext, zwb_fk_b ) 
     900      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,              & 
     901         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   & 
     902         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
     903      ! 
     904      CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zhbl,   & 
     905         &                               zwb_ent, zdbdz_bl_ext, zwb_fk_b ) 
    890906      ! 
    891907      ! Reset l_pyc before calculating terms in the flux-gradient relationship 
    892908      DO_2D( 0, 0, 0, 0 ) 
    893          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR.   & 
     909         IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR.   & 
    894910            & nbld(ji,jj) - nmld(ji,jj) == 1   .OR. zdhdt(ji,jj) < 0.0_wp ) THEN   ! ag 19/03 
    895911            l_pyc(ji,jj) = .FALSE.   ! ag 19/03 
     
    920936      jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    921937      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    & 
    922          &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml,   & 
    923          &                           zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
     938         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
     939         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 
    924940      ! 
    925941      CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     
    930946            & 'After averaging, with new hbl (&correct jp_ext), hml: zt_bl=', av_t_bl(ji,jj),& 
    931947            & ' zs_bl=', av_s_bl(ji,jj),  ' zb_bl=', av_b_bl(ji,jj),& 
    932             & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj),  ' zdb_bl=', zdb_bl(ji,jj),& 
     948            & 'zdt_bl=', av_dt_bl(ji,jj), ' zds_bl=', av_ds_bl(ji,jj),  ' zdb_bl=', av_db_bl(ji,jj),& 
    933949            & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj),  ' zb_ml=', av_b_ml(ji,jj),& 
    934             & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj),  ' zdb_ml=', zdb_ml(ji,jj),& 
    935             & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 
    936             & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 
     950            & 'zdt_ml=', av_dt_ml(ji,jj), ' zds_ml=', av_ds_ml(ji,jj),  ' zdb_ml=', av_db_ml(ji,jj),& 
     951            & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', av_du_bl(ji,jj), ' zdv_bl=', av_dv_bl(ji,jj),& 
     952            & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', av_du_ml(ji,jj), ' zdv_ml=', av_dv_ml(ji,jj) 
    937953         FLUSH(narea+100) 
    938954      END IF 
     
    940956      ! Rotate mean currents and changes onto wind aligned co-ordinates 
    941957      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml ) 
    942       CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 
     958      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml ) 
    943959      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl ) 
    944       CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 
     960      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl ) 
    945961#ifdef key_osm_debug 
    946962      IF(narea==nn_narea_db) THEN 
     
    948964         WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & 
    949965            & 'After rotation, with new hbl (& correct jp_ext), hml:', & 
    950             & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 
    951             & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 
     966            & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', av_du_bl(ji,jj), ' zdv_bl=', av_dv_bl(ji,jj),& 
     967            & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', av_du_ml(ji,jj), ' zdv_ml=', av_dv_ml(ji,jj) 
    952968         FLUSH(narea+100) 
    953969      END IF 
     
    958974      jp_ext(:,:) = 1   ! ag 19/03 
    959975      CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, jp_ext, zdbdz_pyc, zalpha_pyc, zdh,    & 
    960          &                                       zdb_bl, zhbl, zdbdz_bl_ext, zhml, zdb_ml,   & 
    961          &                                       zdhdt ) 
     976         &                                       zhbl, zdbdz_bl_ext, zhml, zdhdt ) 
    962977#ifdef key_osm_debug 
    963978      IF(narea==nn_narea_db) THEN 
     
    9961011      ! Calculate non-gradient components of the flux-gradient relationships 
    9971012      ! -------------------------------------------------------------------- 
    998       CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                            & 
    999          &                    zdhdt, zshear, zdt_bl, zds_bl, zdb_bl,                   & 
    1000          &                    zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml,                  & 
    1001          &                    zdu_ml, zdv_ml, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc,   & 
     1013      CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                           & 
     1014         &                    zdhdt, zshear, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc,   & 
    10021015         &                    zalpha_pyc, zdiffut, zviscos ) 
    10031016      ! 
     
    11861199         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                     ! Boundary-layer depth 
    11871200         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )                  ! Boundary-layer max k 
    1188          IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )            ! dt at ml base 
    1189          IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )            ! ds at ml base 
    1190          IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )            ! db at ml base 
    1191          IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )            ! du at ml base 
    1192          IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )            ! dv at ml base 
     1201         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*av_dt_bl )          ! dt at ml base 
     1202         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*av_ds_bl )            ! ds at ml base 
     1203         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*av_db_bl )            ! db at ml base 
     1204         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*av_du_bl )            ! du at ml base 
     1205         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*av_dv_bl )            ! dv at ml base 
    11931206         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )                        ! Initial boundary-layer depth 
    11941207         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )                     ! Initial boundary-layer depth 
    1195          IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml )            ! dt at ml base 
    1196          IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml )            ! ds at ml base 
    1197          IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml )            ! db at ml base 
     1208         IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*av_dt_ml )            ! dt at ml base 
     1209         IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*av_ds_ml )            ! ds at ml base 
     1210         IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*av_db_ml )            ! db at ml base 
    11981211         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )         ! Stokes drift penetration depth 
    11991212         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*sustke )            ! Stokes drift magnitude at T-points 
     
    14261439 
    14271440   SUBROUTINE zdf_osm_osbl_state( Kmm, pwb_ent, pwb_min, pshear, phbl,   & 
    1428       &                           phml, pdh, pdb_bl, pdu_bl, pdv_bl,     & 
    1429       &                           pdb_ml, pdu_ml, pdv_ml ) 
     1441      &                           phml, pdh ) 
    14301442      !!--------------------------------------------------------------------- 
    14311443      !!                 ***  ROUTINE zdf_osm_osbl_state  *** 
     
    14461458      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phml      ! ML depth 
    14471459      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdh       ! Pycnocline depth 
    1448       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdb_bl    ! Buoyancy diff. between BL average and basal value 
    1449       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdu_bl    ! Velocity diff. (u) between BL average and basal value 
    1450       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdv_bl    ! Velocity diff. (v) between BL average and basal value 
    1451       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdb_ml    ! Buoyancy diff. between mixed-layer average and basal value 
    1452       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdu_ml    ! Velocity diff. (u) between mixed-layer average and basal value 
    1453       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdv_ml    ! Velocity diff. (v) between mixed-layer average and basal value 
    14541460      ! 
    14551461      ! Local variables 
     
    15021508      DO_2D( 0, 0, 0, 0 ) 
    15031509         IF ( l_conv(ji,jj) ) THEN 
    1504             IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
    1505                zri_p(ji,jj) = MAX (  SQRT( pdb_bl(ji,jj) * pdh(ji,jj) / MAX( pdu_bl(ji,jj)**2 + pdv_bl(ji,jj)**2, 1e-8_wp ) ) *   & 
    1506                   &                  ( phbl(ji,jj) / pdh(ji,jj) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 /         & 
    1507                   &                  MAX( zekman(ji,jj), 1.0e-6_wp )  , 5.0_wp ) 
     1510            IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
     1511               zri_p(ji,jj) = MAX (  SQRT( av_db_bl(ji,jj) * pdh(ji,jj) / MAX( av_du_bl(ji,jj)**2 + av_dv_bl(ji,jj)**2,       & 
     1512                  &                                                          1e-8_wp ) ) * ( phbl(ji,jj) / pdh(ji,jj) ) *   & 
     1513                  &                  ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 /                                  & 
     1514                  &                  MAX( zekman(ji,jj), 1.0e-6_wp ), 5.0_wp ) 
    15081515               IF ( ff_t(ji,jj) >= 0.0_wp ) THEN   ! Northern hemisphere 
    1509                   zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 +   & 
    1510                      &                                          MAX( -1.0_wp * pdv_ml(ji,jj), 1e-5_wp)**2 ) 
     1516                  zri_b(ji,jj) = av_db_ml(ji,jj) * pdh(ji,jj) / ( MAX( av_du_ml(ji,jj), 1e-5_wp )**2 +   & 
     1517                     &                                          MAX( -1.0_wp * av_dv_ml(ji,jj), 1e-5_wp)**2 ) 
    15111518               ELSE                                ! Southern hemisphere 
    1512                   zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 +   & 
    1513                      &                                          MAX(           pdv_ml(ji,jj), 1e-5_wp)**2 ) 
     1519                  zri_b(ji,jj) = av_db_ml(ji,jj) * pdh(ji,jj) / ( MAX( av_du_ml(ji,jj), 1e-5_wp )**2 +   & 
     1520                     &                                          MAX(           av_dv_ml(ji,jj), 1e-5_wp)**2 ) 
    15141521               END IF 
    15151522               pshear(ji,jj) = pp_a_shr * zekman(ji,jj) *                                                   & 
    1516                   &            ( MAX( sustar(ji,jj)**2 * pdu_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) +          & 
     1523                  &            ( MAX( sustar(ji,jj)**2 * av_du_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) +          & 
    15171524                  &              pp_b_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) *   & 
    1518                   &                            pdv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) ) 
     1525                  &                            av_dv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) ) 
    15191526#ifdef key_osm_debug 
    15201527               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     
    15561563               !               l_shear(ji,jj) = .FALSE. 
    15571564               !            ENDIF 
    1558             ELSE   ! pdb_bl test, note pshear set to zero 
     1565            ELSE   ! av_db_bl test, note pshear set to zero 
    15591566               n_ddh(ji,jj) = 2 
    15601567               l_shear(ji,jj) = .FALSE. 
     
    16051612                  ! Developing shear layer, additional shear production possible. 
    16061613 
    1607                   !              pshear_u = MAX( zustar(ji,jj)**2 * MAX( pdu_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp ) 
     1614                  !              pshear_u = MAX( zustar(ji,jj)**2 * MAX( av_du_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp ) 
    16081615                  !              pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1.d0 )**2 ) 
    16091616                  !              pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u ) 
     
    16871694   END SUBROUTINE zdf_osm_external_gradients 
    16881695 
    1689    SUBROUTINE zdf_osm_calculate_dhdt( pdhdt, phbl, pdh, pwb_ent, pwb_min,               & 
    1690       &                               pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk,   & 
    1691       &                               pvel_mle ) 
     1696   SUBROUTINE zdf_osm_calculate_dhdt( pdhdt, phbl, pdh, pwb_ent, pwb_min,   & 
     1697      &                               pdbdz_bl_ext, pwb_fk_b, pwb_fk, pvel_mle ) 
    16921698      !!--------------------------------------------------------------------- 
    16931699      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     
    17031709      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    17041710      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_min 
    1705       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
    17061711      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    1707       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
    17081712      REAL(wp), DIMENSION(:,:),    INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
    17091713      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
     
    17401744                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening, 
    17411745                     !                                                                 !    entrainment > restratification 
    1742                      IF ( pdb_bl(ji,jj) > 1e-15_wp ) THEN 
    1743                         zgamma_b_nd = MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) * pdh(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1746                     IF ( av_db_bl(ji,jj) > 1e-15_wp ) THEN 
     1747                        zgamma_b_nd = MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) * pdh(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    17441748                        zpsi = ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
    17451749                           &   ( swb0(ji,jj) - MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) ) * pdh(ji,jj) / phbl(ji,jj) 
     
    17481752                        zpsi = pp_alpha_b * MAX( zpsi, 0.0_wp ) 
    17491753                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   & 
    1750                            &                      ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) +   & 
    1751                            &            zpsi / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1754                           &                      ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) +   & 
     1755                           &            zpsi / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    17521756#ifdef key_osm_debug 
    17531757                        IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     
    17591763                        IF ( n_ddh(ji,jj) == 1 ) THEN 
    17601764                           IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN 
    1761                               zari = MIN( 1.5_wp * pdb_bl(ji,jj) /                                                   & 
     1765                              zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                   & 
    17621766                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                     & 
    1763                                  &                               pdb_bl(ji,jj)**2 / MAX( 4.5_wp * svstr(ji,jj)**2,   & 
     1767                                 &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * svstr(ji,jj)**2,   & 
    17641768                                 &                                                       1e-12_wp ) ) ), 0.2_wp ) 
    17651769                           ELSE 
    1766                               zari = MIN( 1.5_wp * pdb_bl(ji,jj) /                                                    & 
     1770                              zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                    & 
    17671771                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                      & 
    1768                                  &                               pdb_bl(ji,jj)**2 / MAX( 4.5_wp * swstrc(ji,jj)**2,   & 
     1772                                 &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * swstrc(ji,jj)**2,   & 
    17691773                                 &                                                       1e-12_wp ) ) ), 0.2_wp ) 
    17701774                           ENDIF 
    17711775                           ! Relaxation to dh_ref = zari * hbl 
    17721776                           zddhdt = -1.0_wp * pp_ddh_2 * ( 1.0_wp - pdh(ji,jj) / ( zari * phbl(ji,jj) ) ) * pwb_ent(ji,jj) /   & 
    1773                               &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1777                              &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    17741778#ifdef key_osm_debug 
    17751779                           IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     
    17801784                        ELSE IF ( n_ddh(ji,jj) == 0 ) THEN   ! Growing shear layer 
    17811785                           zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
    1782                               &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1786                              &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    17831787                           zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt 
    17841788                        ELSE 
     
    17861790                        ENDIF   ! n_ddh 
    17871791                        pdhdt(ji,jj) = pdhdt(ji,jj) + pp_alpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
    1788                            &                            pdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) /   & 
    1789                            &                            ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
    1790                      ELSE   ! pdb_bl >0 
     1792                           &                            av_db_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) /   & 
     1793                           &                            ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
     1794                     ELSE   ! av_db_bl >0 
    17911795                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1e-15_wp ) 
    17921796                     ENDIF 
     
    17991803                     &                                                         rn_Dt / hbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
    18001804                     &       MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln ) 
    1801                   pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1805                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    18021806                  ! added ajgn 23 July as temporay fix 
    18031807               ENDIF   ! ln_osm_mle 
     
    18091813                  zpert = 2.0_wp * ( 1.0_wp + 0.0_wp * 2.0_wp * svstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * svstr(ji,jj)**2 / hbl(ji,jj) 
    18101814               ELSE 
    1811                   zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), pdb_bl(ji,jj) ) 
     1815                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), av_db_bl(ji,jj) ) 
    18121816               ENDIF 
    18131817               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX( zpert, epsln ) 
     
    18311835                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening, 
    18321836                     !                                                                 !    entrainment > restratification 
    1833                      IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
     1837                     IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
    18341838                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   & 
    1835                            &            ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1839                           &            ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    18361840                     ELSE 
    18371841                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) / MAX( zvel_max, 1e-15_wp ) 
     
    18421846               ELSE   ! Fox-Kemper not used 
    18431847                  zvel_max = -1.0_wp * pwb_ent(ji,jj) / MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln ) 
    1844                   pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     1848                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    18451849                  ! added ajgn 23 July as temporay fix 
    18461850               ENDIF  ! ln_osm_mle 
     
    18531857                  zpert = 2.0_wp * svstr(ji,jj)**2 / hbl(ji,jj) 
    18541858               ELSE 
    1855                   zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), pdb_bl(ji,jj) ) 
     1859                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), av_db_bl(ji,jj) ) 
    18561860               ENDIF 
    18571861               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX(zpert, epsln) 
     
    20202024   END SUBROUTINE zdf_osm_timestep_hbl 
    20212025 
    2022    SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, pdh, phml, pdhdt, pdb_bl,   & 
    2023       &                                     phbl, pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 
     2026   SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, pdh, phml, pdhdt, phbl,   & 
     2027      &                                     pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 
    20242028      !!--------------------------------------------------------------------- 
    20252029      !!            ***  ROUTINE zdf_osm_pycnocline_thickness  *** 
     
    20382042      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phml           ! ML depth 
    20392043      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2040       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
    20412044      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    20422045      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     
    20612064            IF ( l_conv(ji,jj) ) THEN 
    20622065               ! 
    2063                IF ( pdb_bl(ji,jj) > 1e-15_wp ) THEN 
     2066               IF ( av_db_bl(ji,jj) > 1e-15_wp ) THEN 
    20642067                  IF ( n_ddh(ji,jj) == 0 ) THEN 
    20652068                     zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
    20662069                     ! ddhdt for pycnocline determined in osm_calculate_dhdt 
    20672070                     zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
    2068                         &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15 ) ) 
     2071                        &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15 ) ) 
    20692072                     zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt 
    20702073                     ! Maximum limit for how thick the shear layer can grow relative to the thickness of the boundary layer 
     
    20762079                        ztmp = swstrc(ji,jj) 
    20772080                     END IF 
    2078                      zari = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +        & 
    2079                         &                                                   pdb_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2,   & 
     2081                     zari = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +        & 
     2082                        &                                                   av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2,   & 
    20802083                        &                                                                           1e-12_wp ) ) ), 0.2_wp ) 
    2081                      ztau = MAX( pdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( pp_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ),   & 
    2082                         &        2.0_wp * rn_Dt ) 
     2084                     ztau = MAX( av_db_bl(ji,jj) * ( zari * hbl(ji,jj) ) /   & 
     2085                        &        ( pp_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ), 2.0_wp * rn_Dt ) 
    20832086                     dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) +   & 
    20842087                        &        zari * phbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 
     
    20972100               IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN   ! Probably shouldn't include wm here 
    20982101                  ! Boundary layer deepening 
    2099                   IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
     2102                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
    21002103                     ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions 
    2101                      zari    = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( pdb_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp, 0.2_wp ) 
     2104                     zari    = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( av_db_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp, 0.2_wp ) 
    21022105                     zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj) 
    21032106                  ELSE 
     
    21192122                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening. Note wb_fk_b is zero if 
    21202123                     !                                                                 !    ln_osm_mle=F 
    2121                      IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
     2124                     IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
    21222125                        IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln) )**3 <= 0.5_wp ) THEN   ! Near neutral stability 
    21232126                           ztmp = svstr(ji,jj) 
     
    21252128                           ztmp = swstrc(ji,jj) 
    21262129                        END IF 
    2127                         zari = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) *                           & 
     2130                        zari = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                           & 
    21282131                           &                                   ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp) +   & 
    2129                            &                                     pdb_bl(ji,jj)**2 / MAX(4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 
     2132                           &                                     av_db_bl(ji,jj)**2 / MAX(4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 
    21302133                     ELSE 
    21312134                        zari = 0.2_wp 
     
    21372140                  zdh_ref = zari * hbl(ji,jj) 
    21382141               ELSE   ! ln_osm_mle 
    2139                   IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
     2142                  IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
    21402143                     IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln ) )**3 <= 0.5_wp ) THEN   ! Near neutral stability 
    21412144                        ztmp = svstr(ji,jj) 
     
    21432146                        ztmp = swstrc(ji,jj) 
    21442147                     END IF 
    2145                      zari    = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) *                            & 
     2148                     zari    = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                            & 
    21462149                        &                                      ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +   & 
    2147                         &                                        pdb_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 
     2150                        &                                        av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 
    21482151                  ELSE 
    21492152                     zari    = 0.2_wp 
     
    21612164               IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN   ! Probably shouldn't include wm here 
    21622165                  ! Boundary layer deepening 
    2163                   IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
     2166                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
    21642167                     ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    2165                      zari    = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( pdb_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp , 0.2_wp ) 
     2168                     zari    = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( av_db_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp , 0.2_wp ) 
    21662169                     zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj) 
    21672170                  ELSE 
     
    21982201   END SUBROUTINE zdf_osm_pycnocline_thickness 
    21992202 
    2200    SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, pdbdz, palpha, pdh,            & 
    2201       &                                             pdb_bl, phbl, pdbdz_bl_ext, phml, pdb_ml,   & 
    2202       &                                             pdhdt ) 
     2203   SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, pdbdz, palpha, pdh,   & 
     2204      &                                             phbl, pdbdz_bl_ext, phml, pdhdt ) 
    22032205      !!--------------------------------------------------------------------- 
    22042206      !!       ***  ROUTINE zdf_osm_pycnocline_buoyancy_profiles  *** 
     
    22142216      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   palpha 
    22152217      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline thickness 
    2216       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
    22172218      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    22182219      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    22192220      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
    2220       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
    22212221      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! Rates of change of hbl 
    22222222      ! 
     
    22422242                  zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 
    22432243                  palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / pp_gamma_b ) ) *   & 
    2244                      &                                pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / pdb_ml(ji,jj) ) /                & 
     2244                     &                                pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / av_db_ml(ji,jj) ) /                & 
    22452245                     &            ( 0.723_wp + SQRT( 3.14159_wp / pp_gamma_b ) ) 
    22462246                  palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) 
     
    22502250                  ! can be removed                                                          ! 
    22512251                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    2252                   ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
    2253                   ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
    2254                   zbgrad = palpha(ji,jj) * pdb_ml(ji,jj) * ztmp + pdbdz_bl_ext(ji,jj) 
    2255                   zgamma_b_nd = pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / MAX( pdb_ml(ji,jj), epsln ) 
     2252                  ! ztgrad = zalpha * av_dt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
     2253                  ! zsgrad = zalpha * av_ds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
     2254                  zbgrad = palpha(ji,jj) * av_db_ml(ji,jj) * ztmp + pdbdz_bl_ext(ji,jj) 
     2255                  zgamma_b_nd = pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / MAX( av_db_ml(ji,jj), epsln ) 
    22562256                  DO jk = 2, nbld(ji,jj) 
    22572257                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) * ztmp 
    22582258                     IF ( znd <= zzeta_m ) THEN 
    2259                         ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 
     2259                        ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * av_dt_ml(ji,jj) * ztmp * & 
    22602260                        ! &        EXP( -6.0 * ( znd -zzeta_m )**2 ) 
    2261                         ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 
     2261                        ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * av_ds_ml(ji,jj) * ztmp * & 
    22622262                        ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
    2263                         pdbdz(ji,jj,jk) = pdbdz_bl_ext(ji,jj) + palpha(ji,jj) * pdb_ml(ji,jj) * ztmp * & 
     2263                        pdbdz(ji,jj,jk) = pdbdz_bl_ext(ji,jj) + palpha(ji,jj) * av_db_ml(ji,jj) * ztmp * & 
    22642264                           & EXP( -6.0_wp * ( znd -zzeta_m )**2 ) 
    22652265                     ELSE 
     
    22822282               ! If pycnocline profile only defined when depth steady of increasing. 
    22832283               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
    2284                   IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
     2284                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
    22852285                     IF ( shol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline 
    22862286                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 
    2287                         zbgrad = pdb_bl(ji,jj) * ztmp 
     2287                        zbgrad = av_db_bl(ji,jj) * ztmp 
    22882288                        DO jk = 2, nbld(ji,jj) 
    22892289                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
     
    22922292                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    22932293                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 
    2294                         zbgrad = pdb_bl(ji,jj) * ztmp 
     2294                        zbgrad = av_db_bl(ji,jj) * ztmp 
    22952295                        DO jk = 2, nbld(ji,jj) 
    22962296                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
     
    23062306                     END IF 
    23072307#endif 
    2308                   END IF      ! IF (pdb_bl> 0.) 
     2308                  END IF      ! IF (av_db_bl> 0.) 
    23092309               END IF         ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    23102310               ! 
     
    25822582   END SUBROUTINE zdf_osm_diffusivity_viscosity 
    25832583 
    2584    SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,                            & 
    2585       &                          pdhdt, pshear, pdt_bl, pds_bl, pdb_bl,                   & 
    2586       &                          pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml,                  & 
    2587       &                          pdu_ml, pdv_ml, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc,   & 
     2584   SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,                           & 
     2585      &                          pdhdt, pshear, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc,   & 
    25882586      &                          palpha_pyc, pdiffut, pviscos ) 
    25892587      !!--------------------------------------------------------------------- 
     
    26022600      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    26032601      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pshear         ! Shear production 
    2604       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdt_bl         ! Temperature diff. between BL average and basal value 
    2605       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pds_bl         ! Salinity diff. between BL average and basal value 
    2606       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
    2607       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdu_bl         ! Velocity diff. (u) between BL average and basal value 
    2608       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdv_bl         ! Velocity diff. (v) between BL average and basal value 
    2609       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdt_ml         ! Temperature diff. between mixed-layer average and basal value 
    2610       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pds_ml         ! Salinity diff. between mixed-layer average and basal value 
    2611       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
    2612       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdu_ml         ! Velocity diff. (u) between mixed-layer average and basal value 
    2613       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdv_ml         ! Velocity diff. (v) between mixed-layer average and basal value 
    26142602      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
    26152603      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
     
    27882776               &                ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.5_wp ) 
    27892777            zwth_ent(ji,jj)     = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *   & 
    2790                &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj) 
     2778               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dt_ml(ji,jj) 
    27912779            zws_ent(ji,jj)      = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *   & 
    2792                &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 
     2780               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_ds_ml(ji,jj) 
    27932781            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 
    2794                zbuoy_pyc_sc        = 2.0_wp * MAX( pdb_ml(ji,jj), 0.0_wp ) / pdh(ji,jj) 
     2782               zbuoy_pyc_sc        = 2.0_wp * MAX( av_db_ml(ji,jj), 0.0_wp ) / pdh(ji,jj) 
    27952783               zdelta_pyc          = ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird /   & 
    27962784                  &                       SQRT( MAX( zbuoy_pyc_sc, ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 
    2797                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) ) *   & 
     2785               zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * av_dt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   & 
    27982786                  &                     zdelta_pyc**2 / pdh(ji,jj) 
    2799                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) ) *   & 
     2787               zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * av_ds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   & 
    28002788                  &                     zdelta_pyc**2 / pdh(ji,jj) 
    28012789               zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 
     
    28832871               ! Place holding code. Parametrization needs checking for these conditions. 
    28842872               zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
    2885                zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
    2886                zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
     2873               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_du_ml(ji,jj) 
     2874               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dv_ml(ji,jj) 
    28872875            ELSE 
    28882876               zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 
    2889                zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
    2890                zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
     2877               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_du_ml(ji,jj) 
     2878               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dv_ml(ji,jj) 
    28912879            ENDIF 
    28922880            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * suw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj) 
     
    29412929         zsc_ws_1(:,:)  = sws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 
    29422930         WHERE ( l_pyc(A2D(0)) )   ! Pycnocline scales 
    2943             zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pdt_ml(A2D(0)) 
    2944             zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pds_ml(A2D(0)) 
     2931            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_dt_ml(A2D(0)) 
     2932            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_ds_ml(A2D(0)) 
    29452933         END WHERE 
    29462934      ELSEWHERE 
     
    30773065         IF ( l_conv (ji,jj) ) THEN 
    30783066            ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
    3079             !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     3067            !                  zugrad = 0.7 * av_du_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
    30803068            !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
    30813069            !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
    30823070            !Alan is this right? 
    3083             !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     3071            !                  zvgrad = ( 0.7 * av_dv_ml(ji,jj) + & 
    30843072            !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
    30853073            !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     
    30993087               ! Pycnocline profile only defined when depth steady of increasing. 
    31003088               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
    3101                   IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
     3089                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
    31023090                     IF ( shol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline 
    31033091                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 
    3104                         ztgrad = pdt_bl(ji,jj) * ztmp 
    3105                         zsgrad = pds_bl(ji,jj) * ztmp 
    3106                         zbgrad = pdb_bl(ji,jj) * ztmp 
     3092                        ztgrad = av_dt_bl(ji,jj) * ztmp 
     3093                        zsgrad = av_ds_bl(ji,jj) * ztmp 
     3094                        zbgrad = av_db_bl(ji,jj) * ztmp 
    31073095                        IF ( jk <= nbld(ji,jj) ) THEN 
    31083096                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
     
    31183106                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    31193107                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 
    3120                         ztgrad = pdt_bl(ji,jj) * ztmp 
    3121                         zsgrad = pds_bl(ji,jj) * ztmp 
    3122                         zbgrad = pdb_bl(ji,jj) * ztmp 
     3108                        ztgrad = av_dt_bl(ji,jj) * ztmp 
     3109                        zsgrad = av_ds_bl(ji,jj) * ztmp 
     3110                        zbgrad = av_db_bl(ji,jj) * ztmp 
    31233111                        IF ( jk <= nbld(ji,jj) ) THEN 
    31243112                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
     
    31333121                        END IF 
    31343122                     ENDIF   ! IF (shol >=0.5) 
    3135                   ENDIF      ! IF (zdb_bl> 0.) 
     3123                  ENDIF      ! IF (av_db_bl> 0.) 
    31363124               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    31373125            END IF 
     
    31453133         IF ( .NOT. l_conv (ji,jj) ) THEN 
    31463134            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    3147                zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj) 
    3148                zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj) 
     3135               zugrad = 3.25_wp * av_du_bl(ji,jj) / phbl(ji,jj) 
     3136               zvgrad = 2.75_wp * av_dv_bl(ji,jj) / phbl(ji,jj) 
    31493137               IF ( jk <= nbld(ji,jj) ) THEN 
    31503138                  znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
     
    33093297 
    33103298   SUBROUTINE zdf_osm_osbl_state_fk( Kmm, pwb_fk, phbl, phmle, pwb_ent,   & 
    3311       &                              pdb_bl, pdbds_mle ) 
     3299      &                              pdbds_mle ) 
    33123300      !!--------------------------------------------------------------------- 
    33133301      !!               ***  ROUTINE zdf_osm_osbl_state_fk  *** 
     
    33333321      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phmle       ! MLE depth 
    33343322      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux 
    3335       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdb_bl      ! Difference between blayer average and parameter at base of blayer 
    33363323      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    33373324      ! 
     
    34033390                  ENDIF  ! znd_param > 100 
    34043391                  ! 
    3405                   IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
     3392                  IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) THEN 
    34063393                     l_pyc(ji,jj) = .FALSE. 
    34073394                  ELSE 
     
    34093396                  ENDIF 
    34103397               ELSE   ! MLE layer restricted to OSBL or just below 
    3411                   IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN   ! Weak stratification MLE layer can grow 
     3398                  IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) THEN   ! Weak stratification MLE layer can grow 
    34123399                     l_pyc(ji,jj)  = .FALSE. 
    34133400                     l_flux(ji,jj) = .TRUE. 
     
    34173404                     l_flux(ji,jj) = .FALSE. 
    34183405                     l_mle(ji,jj)  = .FALSE. 
    3419                   END IF   ! pdb_bl < rn_mle_thresh_bl and 
     3406                  END IF   ! av_db_bl < rn_mle_thresh_bl and 
    34203407               END IF   ! phmle > 1.2 phbl 
    34213408            ELSE 
     
    34233410               l_flux(ji,jj) = .FALSE. 
    34243411               l_mle(ji,jj)  = .FALSE. 
    3425                IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
     3412               IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
    34263413            END IF   !  -2.0 * pwb_fk(ji,jj) / pwb_ent > 0.5 
    34273414         ELSE   ! Stable Boundary Layer 
Note: See TracChangeset for help on using the changeset viewer.