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

Changeset 14816 for NEMO/branches/2021


Ignore:
Timestamp:
2021-05-07T18:51:36+02:00 (4 years ago)
Author:
smueller
Message:

Slight rearrangement of the subroutine structure and removal of the halo regions from the majority of arrays in module zdfosm (ticket #2353)

File:
1 edited

Legend:

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

    r14803 r14816  
    4949   !!      zdf_osm_timestep_hbl                 : hbl timestep 
    5050   !!      zdf_osm_pycnocline_thickness         : calculate thickness of pycnocline 
    51    !!      zdf_osm_pycnocline_buoyancy_profiles : calculate pycnocline buoyancy profiles 
    5251   !!      zdf_osm_diffusivity_viscosity        : compute eddy diffusivity and viscosity profiles 
    5352   !!      zdf_osm_fgr_terms                    : compute flux-gradient relationship terms 
     53   !!         zdf_osm_pycnocline_buoyancy_profiles : calculate pycnocline buoyancy profiles 
    5454   !!      zdf_osm_zmld_horizontal_gradients    : calculate horizontal buoyancy gradients for use with Fox-Kemper parametrization 
    5555   !!      zdf_osm_osbl_state_fk                : determine state of OSBL and MLE layers 
     
    193193   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_mle  ! Buoyancy average 
    194194   ! 
     195   ! Diagnostic output 
     196   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   osmdia2d   ! Auxiliary array for diagnostic output 
     197   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   osmdia3d   ! Auxiliary array for diagnostic output 
     198   ! 
    195199   !            ** Namelist  namzdf_osm  ** 
    196200   LOGICAL  ::   ln_use_osm_la                      ! Use namelist rn_osm_la 
     
    264268      zdf_osm_alloc = zdf_osm_alloc + ierr 
    265269      ! 
    266       ALLOCATE( etmean(jpi,jpj,jpk), dh(jpi,jpj), r1_ft(jpi,jpj), STAT=ierr ) 
     270      ALLOCATE( etmean(A2D(0),jpk), dh(jpi,jpj), r1_ft(A2D(0)), STAT=ierr ) 
    267271      zdf_osm_alloc = zdf_osm_alloc + ierr 
    268272      ! 
    269       ALLOCATE( nbld(jpi,jpj), nmld(jpi,jpj), STAT=ierr ) 
     273      ALLOCATE( nbld(jpi,jpj), nmld(A2D(0)), STAT=ierr ) 
    270274      zdf_osm_alloc = zdf_osm_alloc + ierr 
    271275      ! 
    272       ALLOCATE( n_ddh(jpi,jpj), STAT=ierr ) 
     276      ALLOCATE( n_ddh(A2D(0)), STAT=ierr ) 
    273277      zdf_osm_alloc = zdf_osm_alloc + ierr 
    274278      ! 
    275       ALLOCATE( l_conv(jpi,jpj), l_shear(jpi,jpj), l_coup(jpi,jpj), l_pyc(jpi,jpj), l_flux(jpi,jpj), l_mle(jpi,jpj), STAT=ierr ) 
     279      ALLOCATE( l_conv(A2D(0)), l_shear(A2D(0)), l_coup(A2D(0)), l_pyc(A2D(0)), l_flux(A2D(0)), l_mle(A2D(0)), STAT=ierr ) 
    276280      zdf_osm_alloc = zdf_osm_alloc + ierr 
    277281      ! 
    278       ALLOCATE( swth0(jpi,jpj),     sws0(jpi,jpj),   swb0(jpi,jpj),  suw0(jpi,jpj),  sustar(jpi,jpj), scos_wind(jpi,jpj),   & 
    279          &      ssin_wind(jpi,jpj), swthav(jpi,jpj), swsav(jpi,jpj), swbav(jpi,jpj), sustke(jpi,jpj), dstokes(jpi,jpj),     & 
    280          &      swstrl(jpi,jpj),    swstrc(jpi,jpj), sla(jpi,jpj),   svstr(jpi,jpj), shol(jpi,jpj),   STAT=ierr ) 
     282      ALLOCATE( swth0(A2D(0)),     sws0(A2D(0)),   swb0(A2D(0)),  suw0(A2D(0)),  sustar(A2D(0)), scos_wind(A2D(0)),   & 
     283         &      ssin_wind(A2D(0)), swthav(A2D(0)), swsav(A2D(0)), swbav(A2D(0)), sustke(A2D(0)), dstokes(A2D(0)),     & 
     284         &      swstrl(A2D(0)),    swstrc(A2D(0)), sla(A2D(0)),   svstr(A2D(0)), shol(A2D(0)),   STAT=ierr ) 
    281285      zdf_osm_alloc = zdf_osm_alloc + ierr 
    282286      ! 
    283       ALLOCATE( av_t_bl(jpi,jpj), av_s_bl(jpi,jpj), av_u_bl(jpi,jpj), av_v_bl(jpi,jpj), av_b_bl(jpi,jpj), STAT=ierr) 
     287      ALLOCATE( av_t_bl(A2D(0)), av_s_bl(A2D(0)), av_u_bl(A2D(0)), av_v_bl(A2D(0)), av_b_bl(A2D(0)), STAT=ierr) 
    284288      zdf_osm_alloc = zdf_osm_alloc + ierr 
    285289      ! 
    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) 
     290      ALLOCATE( av_dt_bl(A2D(0)), av_ds_bl(A2D(0)), av_du_bl(A2D(0)), av_dv_bl(A2D(0)), av_db_bl(A2D(0)), STAT=ierr) 
    287291      zdf_osm_alloc = zdf_osm_alloc + ierr 
    288292      ! 
    289       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) 
     293      ALLOCATE( av_t_ml(A2D(0)), av_s_ml(A2D(0)), av_u_ml(A2D(0)), av_v_ml(A2D(0)), av_b_ml(A2D(0)), STAT=ierr) 
    290294      zdf_osm_alloc = zdf_osm_alloc + ierr 
    291295      ! 
    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) 
     296      ALLOCATE( av_dt_ml(A2D(0)), av_ds_ml(A2D(0)), av_du_ml(A2D(0)), av_dv_ml(A2D(0)), av_db_ml(A2D(0)), STAT=ierr) 
    293297      zdf_osm_alloc = zdf_osm_alloc + ierr 
    294298      ! 
    295       ALLOCATE( av_t_mle(jpi,jpj), av_s_mle(jpi,jpj), av_u_mle(jpi,jpj), av_v_mle(jpi,jpj), av_b_mle(jpi,jpj), STAT=ierr) 
     299      ALLOCATE( av_t_mle(A2D(0)), av_s_mle(A2D(0)), av_u_mle(A2D(0)), av_v_mle(A2D(0)), av_b_mle(A2D(0)), STAT=ierr) 
    296300      zdf_osm_alloc = zdf_osm_alloc + ierr 
     301      ! 
     302      IF ( ln_dia_osm ) THEN 
     303         ALLOCATE( osmdia2d(jpi,jpj), osmdia3d(jpi,jpj,jpk), STAT=ierr ) 
     304         zdf_osm_alloc = zdf_osm_alloc + ierr 
     305      END IF 
    297306      ! 
    298307      CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 
     
    353362      REAL(wp)                     ::   zvw0        ! Surface v-momentum flux 
    354363      REAL(wp), DIMENSION(A2D(0))  ::   zwb0tot     ! Total surface buoyancy flux including insolation 
    355       REAL(wp), DIMENSION(jpi,jpj) ::   zwb_ent     ! Buoyancy entrainment flux 
    356       REAL(wp), DIMENSION(jpi,jpj) ::   zwb_min 
    357       REAL(wp), DIMENSION(jpi,jpj) ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL 
    358       REAL(wp), DIMENSION(jpi,jpj) ::   zwb_fk      ! Max MLE buoyancy flux 
    359       REAL(wp), DIMENSION(jpi,jpj) ::   zdiff_mle   ! Extra MLE vertical diff 
    360       REAL(wp), DIMENSION(jpi,jpj) ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK 
     364      REAL(wp), DIMENSION(A2D(0)) ::   zwb_ent     ! Buoyancy entrainment flux 
     365      REAL(wp), DIMENSION(A2D(0)) ::   zwb_min 
     366      REAL(wp), DIMENSION(A2D(0)) ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL 
     367      REAL(wp), DIMENSION(A2D(0)) ::   zwb_fk      ! Max MLE buoyancy flux 
     368      REAL(wp), DIMENSION(A2D(0)) ::   zdiff_mle   ! Extra MLE vertical diff 
     369      REAL(wp), DIMENSION(A2D(0)) ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK 
    361370      ! 
    362371      ! mixed-layer variables 
    363       INTEGER,  DIMENSION(jpi,jpj) ::   jp_ext   ! Offset for external level 
    364       ! 
    365       REAL(wp), DIMENSION(jpi,jpj) ::   zhbl   ! BL depth - grid 
    366       REAL(wp), DIMENSION(jpi,jpj) ::   zhml   ! ML depth - grid 
    367       ! 
    368       REAL(wp), DIMENSION(jpi,jpj) ::   zhmle   ! MLE depth - grid 
     372      INTEGER,  DIMENSION(A2D(0)) ::   jp_ext   ! Offset for external level 
     373      ! 
     374      REAL(wp), DIMENSION(A2D(0)) ::   zhbl   ! BL depth - grid 
     375      REAL(wp), DIMENSION(A2D(0)) ::   zhml   ! ML depth - grid 
     376      ! 
     377      REAL(wp), DIMENSION(A2D(0)) ::   zhmle   ! MLE depth - grid 
    369378      REAL(wp), DIMENSION(jpi,jpj) ::   zmld   ! ML depth on grid 
    370379      ! 
    371       REAL(wp), DIMENSION(jpi,jpj) ::   zdh   ! Pycnocline depth - grid 
     380      REAL(wp), DIMENSION(A2D(0)) ::   zdh   ! Pycnocline depth - grid 
    372381      REAL(wp), DIMENSION(A2D(0))  ::   zdhdt   ! BL depth tendency 
    373382      REAL(wp), DIMENSION(A2D(0))  ::   zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext   ! External temperature/salinity and buoyancy gradients 
    374       REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization. 
    375       ! 
    376       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline 
    377       REAL(wp), DIMENSION(jpi,jpj)     ::   zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
     383      REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization 
     384      ! 
     385      REAL(wp), DIMENSION(A2D(0)) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    378386      ! Flux-gradient relationship variables 
    379       REAL(wp), DIMENSION(jpi,jpj) ::   zshear   ! Shear production 
    380       ! 
    381       REAL(wp), DIMENSION(jpi,jpj) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep 
     387      REAL(wp), DIMENSION(A2D(0)) ::   zshear   ! Shear production 
     388      ! 
     389      REAL(wp), DIMENSION(A2D(0)) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep 
    382390      ! 
    383391      ! For calculating Ri#-dependent mixing 
     
    387395      ! 
    388396      ! Temporary variables 
    389       REAL(wp)                         ::   znd   ! Temporary non-dimensional depth 
    390       REAL(wp)                         ::   zz0, zz1, zfac 
    391       REAL(wp)                         ::   zus_x, zus_y   ! Temporary Stokes drift 
    392       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zviscos   ! Viscosity 
    393       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdiffut   ! t-diffusivity 
    394       REAL(wp), DIMENSION(jpi,jpj)     ::   zalpha_pyc 
    395       REAL(wp)                         ::   zabsstke 
    396       REAL(wp)                         ::   zsqrtpi, z_two_thirds, zthickness 
    397       REAL(wp)                         ::   z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc 
     397      REAL(wp)                        ::   znd   ! Temporary non-dimensional depth 
     398      REAL(wp)                        ::   zz0, zz1, zfac 
     399      REAL(wp)                        ::   zus_x, zus_y   ! Temporary Stokes drift 
     400      REAL(wp), DIMENSION(A2D(0),jpk) ::   zviscos   ! Viscosity 
     401      REAL(wp), DIMENSION(A2D(0),jpk) ::   zdiffut   ! t-diffusivity 
     402      REAL(wp)                        ::   zabsstke 
     403      REAL(wp)                        ::   zsqrtpi, z_two_thirds, zthickness 
     404      REAL(wp)                        ::   z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc 
    398405      ! 
    399406      ! For debugging 
     
    402409      IF( ln_timing ) CALL timing_start('zdf_osm') 
    403410      ! 
    404       nbld(:,:)   = 0        ; nmld(:,:)  = 0 
    405       sustar(:,:) = pp_large 
    406       swstrl(:,:) = pp_large ; svstr(:,:)      = pp_large ; swstrc(:,:)     = pp_large ; suw0(:,:)      = pp_large 
    407       swth0(:,:)  = pp_large ; sws0(:,:)       = pp_large ; swb0(:,:)       = pp_large 
    408       swthav(:,:) = pp_large ; swsav(:,:)      = pp_large ; swbav(:,:)      = pp_large 
    409       sustke(:,:) = pp_large ; sla(:,:)        = pp_large ; scos_wind(:,:)  = pp_large ; ssin_wind(:,:) = pp_large 
    410       shol(:,:)   = pp_large ; zalpha_pyc(:,:) = pp_large 
    411       l_pyc(:,:)  = .FALSE.  ; l_flux(:,:)     = .FALSE.  ; l_mle(:,:)      = .FALSE. 
     411      nbld(:,:)   = 0 
     412      nmld(:,:)   = 0 
     413      sustke(:,:) = pp_large 
     414      l_pyc(:,:)  = .FALSE. 
     415      l_flux(:,:) = .FALSE. 
     416      l_mle(:,:)  = .FALSE. 
    412417      ! Mixed layer 
    413418      ! No initialization of zhbl or zhml (or zdh?) 
    414419      zhbl(:,:)     = pp_large ; zhml(:,:)     = pp_large ; zdh(:,:)      = pp_large 
    415       av_t_bl(:,:)  = pp_large ; av_s_bl(:,:)  = pp_large ; av_u_bl(:,:)  = pp_large 
    416       av_v_bl(:,:)  = pp_large ; av_b_bl(:,:)  = pp_large ; av_t_ml(:,:)  = pp_large 
    417       av_s_ml(:,:)  = pp_large ; av_u_ml(:,:)  = pp_large ; av_v_ml(:,:)  = pp_large 
    418       av_b_ml(:,:)  = pp_large ; av_t_mle(:,:) = pp_large ; av_s_mle(:,:) = pp_large 
    419       av_u_mle(:,:) = pp_large ; av_v_mle(:,:) = pp_large ; av_b_mle(:,:) = 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 
    424       ! 
    425       zdbdz_pyc(:,:,:)    = pp_large 
    426       zdbdz_pyc(A2D(0),:) = 0.0_wp 
    427420      ! 
    428421      IF ( ln_osm_mle ) THEN   ! Only initialise arrays if needed 
    429422         zdtdx(:,:)  = pp_large ; zdtdy(:,:)    = pp_large ; zdsdx(:,:)     = pp_large 
    430423         zdsdy(:,:)  = pp_large ; dbdx_mle(:,:) = pp_large ; dbdy_mle(:,:)  = pp_large 
    431          zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large ; zdiff_mle(:,:) = pp_large 
     424         zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large 
    432425         zhmle(:,:)  = pp_large ; zmld(:,:)     = pp_large 
    433426      ENDIF 
    434427      zhbl_t(:,:)   = pp_large 
    435428      ! 
    436       zdiffut(:,:,:)    = pp_large ; zviscos(:,:,:)    = pp_large 
    437       zdiffut(A2D(0),:) = 0.0_wp   ; zviscos(A2D(0),:) = 0.0_wp 
     429      zdiffut(:,:,:) = 0.0_wp 
     430      zviscos(:,:,:) = 0.0_wp 
     431      ! 
    438432      ghamt(:,:,:)      = pp_large ; ghams(:,:,:)      = pp_large 
    439433      ghamt(A2D(0),:)   = 0.0_wp   ; ghams(A2D(0),:)   = 0.0_wp 
    440434      ghamu(:,:,:)      = pp_large ; ghamv(:,:,:)      = pp_large 
    441435      ghamu(A2D(0),:)   = 0.0_wp   ; ghamv(A2D(0),:)   = 0.0_wp 
    442       zdiff_mle(A2D(0)) = 0.0_wp 
     436      ! 
     437      zdiff_mle(:,:) = 0.0_wp 
    443438      ! 
    444439#ifdef key_osm_debug 
     
    695690      ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 
    696691      jp_ext(:,:) = 1   ! ag 19/03 
    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,  & 
     692      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
     693         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   & 
    699694         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    700       jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     695      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    701696      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
    702697         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
     
    751746            IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
    752747         END_3D 
    753          CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof, av_t_mle, av_s_mle,   & 
     748         CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D(0)), av_t_mle, av_s_mle,   & 
    754749            &                           av_b_mle, av_u_mle, av_v_mle ) 
    755750         ! 
     
    802797      ! 
    803798      !! External gradient below BL needed both with and w/o FK 
    804       CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
     799      CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
    805800      ! 
    806801      ! Test if pycnocline well resolved 
     
    833828      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 
    834829      jp_ext(:,:) = 1   ! ag 19/03 
    835       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,              & 
     830      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
    836831         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   & 
    837832         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    838       jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     833      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
    839834      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
    840835         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
     
    880875      END_2D 
    881876      ! 
    882       nmld(:,:) = nbld(:,:)           ! use nmld to hold previous blayer index 
     877      nmld(:,:) = nbld(A2D(0))           ! use nmld to hold previous blayer index 
    883878      nbld(:,:) = 4 
    884879      ! 
     
    898893      ! Recalculate BL averages and differences using new BL depth 
    899894      jp_ext(:,:) = 1   ! ag 19/03 
    900       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl,              & 
     895      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
    901896         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   & 
    902897         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
     
    928923      END_2D 
    929924      ! 
    930       dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating 
     925      dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D(0))/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating 
    931926      !                                                       !    flux-gradient terms 
    932927      ! 
     
    934929      !      jp_ext = nbld - nmld + 1 
    935930      ! Recalculate ML averages and differences using new ML depth 
    936       jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03 
     931      jp_ext(:,:) = nbld(A2D(0)) - nmld(A2D(0)) + jp_ext(:,:) + 1   ! ag 19/03 
    937932      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    & 
    938933         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   & 
    939934         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 
    940935      ! 
    941       CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     936      CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    942937#ifdef key_osm_debug 
    943938      IF(narea==nn_narea_db) THEN 
     
    969964      END IF 
    970965#endif 
    971       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    972       !  Pycnocline gradients for scalars and velocity 
    973       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    974       jp_ext(:,:) = 1   ! ag 19/03 
    975       CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, jp_ext, zdbdz_pyc, zalpha_pyc, zdh,    & 
    976          &                                       zhbl, zdbdz_bl_ext, zhml, zdhdt ) 
    977 #ifdef key_osm_debug 
    978       IF(narea==nn_narea_db) THEN 
    979          ji=iloc_db; jj=jloc_db 
    980          jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
    981          WRITE(narea+100,'(a,l7,/,3(a,g11.3),/)') & 
    982             & 'After pycnocline profiles BL  lpyc=', l_pyc(ji,jj),& 
    983             & 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj), & 
    984             & 'Pycnocline: zalpha_pyc=', zalpha_pyc(ji,jj) 
    985          !       WRITE(narea+100,'(a,*(g11.3))') ' zdtdz_pyc[imld-1..ibld+2] =', ( zdtdz_pyc(ji,jj,jk), jk=jl,jm ) 
    986          !       WRITE(narea+100,'(a,*(g11.3))') ' zdsdz_pyc[imld-1..ibld+2] =', ( zdsdz_pyc(ji,jj,jk), jk=jl,jm ) 
    987          WRITE(narea+100,'(a,*(g11.3))') ' zdbdz_pyc[imld-1..ibld+2] =', ( zdbdz_pyc(ji,jj,jk), jk=jl,jm ) 
    988          !       WRITE(narea+100,'(a,*(g11.3))') ' zdudz_pyc[imld-1..ibld+2] =', ( zdudz_pyc(ji,jj,jk), jk=jl,jm ) 
    989          !       WRITE(narea+100,'(a,*(g11.3))') ' zdvdz_pyc[imld-1..ibld+2] =', ( zdvdz_pyc(ji,jj,jk), jk=jl,jm ) 
    990          WRITE(narea+100,*) 
    991          FLUSH(narea+100) 
    992       END IF 
    993 #endif 
    994966      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    995967      ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
     
    1011983      ! Calculate non-gradient components of the flux-gradient relationships 
    1012984      ! -------------------------------------------------------------------- 
    1013       CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                           & 
    1014          &                    zdhdt, zshear, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc,   & 
    1015          &                    zalpha_pyc, zdiffut, zviscos ) 
     985      jp_ext(:,:) = 1   ! ag 19/03 
     986      CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                              & 
     987         &                    zdhdt, zshear, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext,   & 
     988         &                    zdiffut, zviscos ) 
    1016989      ! 
    1017990      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    1022995      ! 
    1023996      ! Rotate non-gradient velocity terms back to model reference frame 
    1024       CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE., 2, nbld ) 
     997      CALL zdf_osm_velocity_rotation( ghamu(A2D(0),:), ghamv(A2D(0),:), .FALSE., 2, nbld(A2D(0)) ) 
    1025998      ! 
    1026999      ! KPP-style Ri# mixing 
     
    11651138#endif 
    11661139      ! 
    1167       IF(ln_dia_osm) THEN 
     1140      IF ( ln_dia_osm ) THEN 
    11681141         SELECT CASE (nn_osm_wave) 
    11691142            ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1) 
    11701143         CASE(0:1) 
    1171             IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*sustke*scos_wind )          ! x surface Stokes drift 
    1172             IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*sustke*ssin_wind )          ! y surface Stokes drift 
    1173             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(:,:,1) *   & 
    1174                &                                                                       sustar**2 * sustke ) 
     1144            IF ( iom_use("us_x") ) THEN                                                           ! x surface Stokes drift 
     1145               osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * scos_wind 
     1146               CALL iom_put( "us_x", osmdia2d ) 
     1147            END IF 
     1148            IF ( iom_use("us_y") ) THEN                                                           ! y surface Stokes drift 
     1149               osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * ssin_wind 
     1150               CALL iom_put( "us_y", osmdia2d ) 
     1151            END IF 
     1152            IF ( iom_use("wind_wave_abs_power") ) THEN 
     1153               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke 
     1154               CALL iom_put( "wind_wave_abs_power", osmdia2d ) 
     1155            END IF 
    11751156            ! Stokes drift read in from sbcwave  (=2). 
    11761157         CASE(2:3) 
     
    11851166               &                                             tmask(:,:,1) )                       !    height from NP spectrum 
    11861167            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                      ! U_10 
    1187             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power",                          & 
    1188                &                                                1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 *   & 
    1189                &                                                SQRT( ut0sd**2 + vt0sd**2 ) ) 
     1168            IF ( iom_use("wind_wave_abs_power") ) THEN 
     1169               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * SQRT( ut0sd**2 + vt0sd**2 ) 
     1170               CALL iom_put( "wind_wave_abs_power", osmdia2d ) 
     1171            END IF 
    11901172         END SELECT 
    11911173         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )                      ! <Tw_NL> 
     
    11931175         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )                      ! <uw_NL> 
    11941176         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )                      ! <vw_NL> 
    1195          IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*swth0 )               ! <Tw_0> 
    1196          IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*sws0 )                  ! <Sw_0> 
    1197          IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*swb0 )                  ! <Sw_0> 
    1198          IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*swth0 )               ! Upward BL-avged turb buoyancy flux 
     1177         IF ( iom_use("zwth0") ) THEN                                                      ! <Tw_0> 
     1178            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwth0",     osmdia2d ) 
     1179         END IF 
     1180         IF ( iom_use("zws0") ) THEN                                                       ! <Sw_0> 
     1181            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sws0;      CALL iom_put( "zws0",      osmdia2d ) 
     1182         END IF 
     1183         IF ( iom_use("zwb0") ) THEN                                                       ! <Sw_0> 
     1184            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swb0;      CALL iom_put( "zwb0",      osmdia2d ) 
     1185         END IF 
     1186         IF ( iom_use("zwbav") ) THEN                                                      ! Upward BL-avged turb buoyancy flux 
     1187            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwbav",     osmdia2d ) 
     1188         END IF 
    11991189         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                     ! Boundary-layer depth 
    12001190         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )                  ! Boundary-layer max k 
    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 
     1191         IF ( iom_use("zdt_bl") ) THEN                                                     ! dt at ml base 
     1192            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_bl;  CALL iom_put( "zdt_bl", osmdia2d ) 
     1193         END IF 
     1194         IF ( iom_use("zds_bl") ) THEN                                                     ! ds at ml base 
     1195            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_bl;  CALL iom_put( "zds_bl", osmdia2d ) 
     1196         END IF 
     1197         IF ( iom_use("zdb_bl") ) THEN                                                     ! db at ml base 
     1198            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_bl;  CALL iom_put( "zdb_bl", osmdia2d ) 
     1199         END IF 
     1200         IF ( iom_use("zdu_bl") ) THEN                                                     ! du at ml base 
     1201            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_du_bl;  CALL iom_put( "zdu_bl", osmdia2d ) 
     1202         END IF 
     1203         IF ( iom_use("zdv_bl") ) THEN                                                     ! dv at ml base 
     1204            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dv_bl;  CALL iom_put( "zdv_bl", osmdia2d ) 
     1205         END IF 
    12061206         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )                        ! Initial boundary-layer depth 
    12071207         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )                     ! Initial boundary-layer depth 
    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 
    1211          IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )         ! Stokes drift penetration depth 
    1212          IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*sustke )            ! Stokes drift magnitude at T-points 
    1213          IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*swstrc )            ! Convective velocity scale 
    1214          IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*swstrl )            ! Langmuir velocity scale 
    1215          IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*sustar )            ! Friction velocity scale 
    1216          IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*svstr )               ! Mixed velocity scale 
    1217          IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*sla )                     ! Langmuir # 
    1218          IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.0_wp * rho0 *   &   ! BL depth internal to zdf_osm routine 
    1219             &                                                     tmask(:,:,1) * sustar**3 )  
    1220          IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * sustke ) 
    1221          IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )                  ! BL depth internal to zdf_osm routine 
    1222          IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )                  ! ML depth internal to zdf_osm routine 
    1223          IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*nmld )                  ! Index for ML depth internal to zdf_osm routine 
    1224          IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext )            ! =1 if pycnocline resolved internal to zdf_osm routine 
    1225          IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*n_ddh )               ! Index forpyc thicknessh internal to zdf_osm routine 
    1226          IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear )            ! Shear production of TKE internal to zdf_osm routine 
    1227          IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                     ! Pyc thicknessh internal to zdf_osm routine 
    1228          IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*shol )                  ! ML depth internal to zdf_osm routine 
    1229          IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )         ! Upward turb buoyancy entrainment flux 
    1230          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*av_t_ml )             ! Average T in ML 
     1208         IF ( iom_use("zdt_ml") ) THEN                                                     ! dt at ml base 
     1209            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_ml;  CALL iom_put( "zdt_ml", osmdia2d ) 
     1210         END IF 
     1211         IF ( iom_use("zds_ml") ) THEN                                                     ! ds at ml base 
     1212            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_ml;  CALL iom_put( "zds_ml", osmdia2d ) 
     1213         END IF 
     1214         IF ( iom_use("zdb_ml") ) THEN                                                     ! db at ml base 
     1215            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_ml;  CALL iom_put( "zdb_ml", osmdia2d ) 
     1216         END IF 
     1217         IF ( iom_use("dstokes") ) THEN                                                    ! Stokes drift penetration depth 
     1218            osmdia2d(A2D(0)) = tmask(A2D(0),1) * dstokes;   CALL iom_put( "dstokes",   osmdia2d ) 
     1219         END IF 
     1220         IF ( iom_use("zustke") ) THEN                                                     ! Stokes drift magnitude at T-points 
     1221            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke;    CALL iom_put( "zustke",    osmdia2d ) 
     1222         END IF 
     1223         IF ( iom_use("zwstrc") ) THEN                                                     ! Convective velocity scale 
     1224            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrc;    CALL iom_put( "zwstrc",    osmdia2d ) 
     1225         END IF 
     1226         IF ( iom_use("zwstrl") ) THEN                                                     ! Langmuir velocity scale 
     1227            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrl;    CALL iom_put( "zwstrl",    osmdia2d ) 
     1228         END IF 
     1229         IF ( iom_use("zustar") ) THEN                                                     ! Friction velocity scale 
     1230            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustar;    CALL iom_put( "zustar",    osmdia2d ) 
     1231         END IF 
     1232         IF ( iom_use("zvstr") ) THEN                                                      ! Mixed velocity scale 
     1233            osmdia2d(A2D(0)) = tmask(A2D(0),1) * svstr;     CALL iom_put( "zvstr",     osmdia2d ) 
     1234         END IF 
     1235         IF ( iom_use("zla") ) THEN                                                        ! Langmuir # 
     1236            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sla;       CALL iom_put( "zla",       osmdia2d ) 
     1237         END IF 
     1238         IF ( iom_use("wind_power") ) THEN                                                 ! BL depth internal to zdf_osm routine 
     1239            osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**3 
     1240            CALL iom_put( "wind_power", osmdia2d ) 
     1241         END IF 
     1242         IF ( iom_use("wind_wave_power") ) THEN 
     1243            osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke 
     1244            CALL iom_put( "wind_wave_power", osmdia2d ) 
     1245         END IF 
     1246         IF ( iom_use("zhbl") ) THEN                                                       ! BL depth internal to zdf_osm routine 
     1247            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhbl;      CALL iom_put( "zhbl",      osmdia2d ) 
     1248         END IF 
     1249         IF ( iom_use("zhml") ) THEN                                                       ! ML depth internal to zdf_osm routine 
     1250            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhml;      CALL iom_put( "zhml",      osmdia2d ) 
     1251         END IF 
     1252         IF ( iom_use("imld") ) THEN                                                       ! Index for ML depth internal to zdf_osm routine 
     1253            osmdia2d(A2D(0)) = tmask(A2D(0),1) * nmld;      CALL iom_put( "imld",      osmdia2d ) 
     1254         END IF 
     1255         IF ( iom_use("jp_ext") ) THEN                                                     ! =1 if pycnocline resolved internal to zdf_osm routine 
     1256            osmdia2d(A2D(0)) = tmask(A2D(0),1) * jp_ext;    CALL iom_put( "jp_ext",    osmdia2d ) 
     1257         END IF 
     1258         IF ( iom_use("j_ddh") ) THEN                                                      ! Index forpyc thicknessh internal to zdf_osm routine 
     1259            osmdia2d(A2D(0)) = tmask(A2D(0),1) * n_ddh;     CALL iom_put( "j_ddh",     osmdia2d ) 
     1260         END IF 
     1261         IF ( iom_use("zshear") ) THEN                                                     ! Shear production of TKE internal to zdf_osm routine 
     1262            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zshear;    CALL iom_put( "zshear",    osmdia2d ) 
     1263         END IF 
     1264         IF ( iom_use("zdh") ) THEN                                                        ! Pyc thicknessh internal to zdf_osm routine 
     1265            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdh;       CALL iom_put( "zdh",       osmdia2d ) 
     1266         END IF 
     1267         IF ( iom_use("zhol") ) THEN                                                       ! ML depth internal to zdf_osm routine 
     1268            osmdia2d(A2D(0)) = tmask(A2D(0),1) * shol;      CALL iom_put( "zhol",      osmdia2d ) 
     1269         END IF 
     1270         IF ( iom_use("zwb_ent") ) THEN                                                    ! Upward turb buoyancy entrainment flux 
     1271            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_ent;   CALL iom_put( "zwb_ent",   osmdia2d ) 
     1272         END IF 
     1273         IF ( iom_use("zt_ml") ) THEN                                                      ! Average T in ML 
     1274            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_t_ml;   CALL iom_put( "zt_ml",     osmdia2d ) 
     1275         END IF 
    12311276         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )                  ! FK layer depth 
    12321277         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )                  ! FK target layer depth 
    1233          IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )            ! FK b flux 
    1234          IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )      ! FK b flux averaged over ML 
     1278         IF ( iom_use("zwb_fk") ) THEN                                                     ! FK b flux 
     1279            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk;    CALL iom_put( "zwb_fk",    osmdia2d ) 
     1280         END IF 
     1281         IF ( iom_use("zwb_fk_b") ) THEN                                                   ! FK b flux averaged over ML 
     1282            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk_b;  CALL iom_put( "zwb_fk_b",  osmdia2d ) 
     1283         END IF 
    12351284         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )      ! FK layer max k 
    12361285         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )               ! FK dtdx at u-pt 
     
    12401289         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )      ! FK dbdx at u-pt 
    12411290         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )      ! FK dbdy at v-pt 
    1242          IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )   ! FK diff in MLE at t-pt 
    1243          IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )     ! FK diff in MLE at t-pt 
    1244  
     1291         IF ( iom_use("zdiff_mle") ) THEN                                                  ! FK diff in MLE at t-pt 
     1292            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d ) 
     1293         END IF 
     1294         IF ( iom_use("zvel_mle") ) THEN                                                   ! FK diff in MLE at t-pt 
     1295            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zvel_mle",  osmdia2d ) 
     1296         END IF 
    12451297      END IF 
    12461298      IF( ln_timing ) CALL timing_stop('zdf_osm') 
     
    12481300   END SUBROUTINE zdf_osm 
    12491301 
    1250    SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt, ps,    & 
     1302   SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt, ps,   & 
    12511303      &                                 pb, pu, pv, kp_ext, pdt,   & 
    12521304      &                                 pds, pdb, pdu, pdv ) 
     
    12621314      !!              knlev+kp_ext 
    12631315      !!---------------------------------------------------------------------- 
    1264       INTEGER,                      INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices 
    1265       INTEGER,  DIMENSION(jpi,jpj), INTENT(in   )           ::   knlev      ! Number of levels to average over. 
    1266       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity 
    1267       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out)           ::   pb         ! Average buoyancy 
    1268       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out)           ::   pu, pv     ! Average current components 
    1269       INTEGER,  DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets 
    1270       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity, 
    1271       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy, 
    1272       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pdu, pdv   ! velocity components and the OSBL 
    1273       ! 
    1274       INTEGER                      ::   jk, jkflt, jkmax, ji, jj   ! Loop indices 
    1275       INTEGER                      ::   ibld_ext                   ! External-layer index 
    1276       REAL(wp), DIMENSION(jpi,jpj) ::   zthick                     ! Layer thickness 
    1277       REAL(wp)                     ::   zthermal, zbeta            ! Thermal/haline expansion/contraction coefficients 
     1316      INTEGER,                     INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices 
     1317      INTEGER,  DIMENSION(A2D(0)), INTENT(in   )           ::   knlev      ! Number of levels to average over. 
     1318      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity 
     1319      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pb         ! Average buoyancy 
     1320      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pu, pv     ! Average current components 
     1321      INTEGER,  DIMENSION(A2D(0)), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets 
     1322      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity, 
     1323      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy, 
     1324      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdu, pdv   ! velocity components and the OSBL 
     1325      ! 
     1326      INTEGER                     ::   jk, jkflt, jkmax, ji, jj   ! Loop indices 
     1327      INTEGER                     ::   ibld_ext                   ! External-layer index 
     1328      REAL(wp), DIMENSION(A2D(0)) ::   zthick                     ! Layer thickness 
     1329      REAL(wp)                    ::   zthermal, zbeta            ! Thermal/haline expansion/contraction coefficients 
    12781330      !!---------------------------------------------------------------------- 
    12791331      ! 
     
    12811333      ! 
    12821334      ! Averages over depth of boundary layer 
    1283       pt(A2D(0))     = 0.0_wp 
    1284       ps(A2D(0))     = 0.0_wp 
    1285       pu(A2D(0))     = 0.0_wp 
    1286       pv(A2D(0))     = 0.0_wp 
     1335      pt(:,:)     = 0.0_wp 
     1336      ps(:,:)     = 0.0_wp 
     1337      pu(:,:)     = 0.0_wp 
     1338      pv(:,:)     = 0.0_wp 
    12871339      zthick(:,:) = epsln 
    12881340      jkflt = jpk 
     
    13671419      !! 
    13681420      !!----------------------------------------------------------------------       
    1369       REAL(wp),           INTENT(inout), DIMENSION(jpi,jpj) ::   pu, pv           ! Components of current 
    1370       LOGICAL,  OPTIONAL, INTENT(in   )                     ::   fwd              ! Forward (default) or reverse rotation 
     1421      REAL(wp),           INTENT(inout), DIMENSION(A2D(0)) ::   pu, pv           ! Components of current 
     1422      LOGICAL,  OPTIONAL, INTENT(in   )                    ::   fwd              ! Forward (default) or reverse rotation 
    13711423      ! 
    13721424      INTEGER  ::   ji, jj       ! Loop indices 
     
    14011453      !! 
    14021454      !!----------------------------------------------------------------------       
    1403       REAL(wp),           INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pu, pv           ! Components of current 
    1404       LOGICAL,  OPTIONAL, INTENT(in   )                         ::   fwd              ! Forward (default) or reverse rotation 
    1405       INTEGER,  OPTIONAL, INTENT(in   )                         ::   ktop             ! Minimum depth index 
    1406       INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(jpi,jpj)     ::   knlev            ! Array of maximum depth indices 
     1455      REAL(wp),           INTENT(inout), DIMENSION(A2D(0),jpk) ::   pu, pv           ! Components of current 
     1456      LOGICAL,  OPTIONAL, INTENT(in   )                        ::   fwd              ! Forward (default) or reverse rotation 
     1457      INTEGER,  OPTIONAL, INTENT(in   )                        ::   ktop             ! Minimum depth index 
     1458      INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D(0))     ::   knlev            ! Array of maximum depth indices 
    14071459      !  
    14081460      INTEGER  ::   ji, jj, jk, jktop, jkmax   ! Loop indices 
     
    14511503      !! 
    14521504      !!---------------------------------------------------------------------- 
    1453       INTEGER,                  INTENT(in   ) ::   Kmm       ! Ocean time-level index 
    1454       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base 
    1455       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_min   !    of well-mixed layer 
    1456       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline 
    1457       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phbl      ! BL depth 
    1458       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phml      ! ML depth 
    1459       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdh       ! Pycnocline depth 
     1505      INTEGER,                     INTENT(in   ) ::   Kmm       ! Ocean time-level index 
     1506      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base 
     1507      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_min   !    of well-mixed layer 
     1508      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline 
     1509      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl      ! BL depth 
     1510      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phml      ! ML depth 
     1511      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh       ! Pycnocline depth 
    14601512      ! 
    14611513      ! Local variables 
     
    14771529      IF( ln_timing ) CALL timing_start('zdf_osm_os') 
    14781530      ! 
    1479       ! Initialise INTENT(  out) arrays 
     1531      ! Initialise arrays 
    14801532      l_conv(:,:)  = .FALSE. 
    14811533      l_shear(:,:) = .FALSE. 
    14821534      n_ddh(:,:)   = 1 
     1535      ! Initialise INTENT(  out) arrays 
    14831536      pwb_ent(:,:) = pp_large 
    14841537      pwb_min(:,:) = pp_large 
    1485       pshear(:,:)  = pp_large 
    14861538      ! 
    14871539      ! Determins stability and set flag l_conv 
     
    14941546      END_2D 
    14951547      ! 
    1496       pshear(A2D(0)) = 0._wp 
    1497       zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(0)) ) * phbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) ) 
     1548      pshear(:,:) = 0.0_wp 
     1549      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(0)) ) * phbl(:,:) / MAX( sustar(A2D(0)), 1.e-8 ) ) 
    14981550      ! 
    14991551#ifdef key_osm_debug 
     
    16581710      !! 
    16591711      !!----------------------------------------------------------------------    
    1660       INTEGER,                      INTENT(in   ) ::   Kmm                   ! Ocean time-level index 
    1661       INTEGER,  DIMENSION(jpi,jpj), INTENT(in   ) ::   kbase                 ! OSBL base layer index 
    1662       REAL(wp), DIMENSION(A2D(0)),  INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy 
     1712      INTEGER,                     INTENT(in   ) ::   Kmm                   ! Ocean time-level index 
     1713      INTEGER,  DIMENSION(A2D(0)), INTENT(in   ) ::   kbase                 ! OSBL base layer index 
     1714      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy 
    16631715      ! 
    16641716      ! Local variables 
     
    17051757      !!---------------------------------------------------------------------- 
    17061758      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl 
    1707       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    1708       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
    1709       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    1710       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_min 
     1759      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth 
     1760      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh            ! Pycnocline depth 
     1761      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     1762      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_min 
    17111763      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    1712       REAL(wp), DIMENSION(:,:),    INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
    1713       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
    1714       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK 
     1764      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
     1765      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
     1766      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK 
    17151767      ! 
    17161768      ! Local variables 
     
    17851837                           zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
    17861838                              &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) 
    1787                            zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt 
     1839                           zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8_wp ) ) * zddhdt 
    17881840                        ELSE 
    17891841                           zddhdt = 0.0_wp 
     
    18981950      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index 
    18991951      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl 
    1900       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phbl       ! BL depth 
    1901       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl_t     ! BL depth 
    1902       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux 
    1903       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL 
     1952      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phbl       ! BL depth 
     1953      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl_t     ! BL depth 
     1954      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux 
     1955      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL 
    19041956      ! 
    19051957      ! Local variables 
     
    20392091      !!---------------------------------------------------------------------- 
    20402092      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    2041       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   pdh            ! Pycnocline thickness 
    2042       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phml           ! ML depth 
     2093      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdh            ! Pycnocline thickness 
     2094      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phml           ! ML depth 
    20432095      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2044       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    2045       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     2096      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth 
     2097      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    20462098      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    2047       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
     2099      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
    20482100 
    20492101      ! 
     
    22112263      !! 
    22122264      !!---------------------------------------------------------------------- 
    2213       INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    2214       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kp_ext         ! External-level offsets 
    2215       REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pdbdz          ! Gradients in the pycnocline 
    2216       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   palpha 
    2217       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline thickness 
    2218       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    2219       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    2220       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
    2221       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! Rates of change of hbl 
     2265      INTEGER,                          INTENT(in   ) ::   Kmm            ! Ocean time-level index 
     2266      INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! External-level offsets 
     2267      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline 
     2268      REAL(wp), DIMENSION(A2D(0)),      INTENT(  out) ::   palpha 
     2269      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline thickness 
     2270      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
     2271      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     2272      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
     2273      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl 
    22222274      ! 
    22232275      ! Local variables 
     
    22292281      ! 
    22302282      REAL(wp), PARAMETER ::   pp_gamma_b = 2.25_wp 
     2283      REAL(wp), PARAMETER ::   pp_large   = -1e10_wp 
    22312284      ! 
    22322285      IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 
     2286      ! 
     2287      pdbdz(:,:,:) = pp_large 
     2288      palpha(:,:)  = pp_large 
    22332289      ! 
    22342290      DO_2D( 0, 0, 0, 0 ) 
     
    22882344                        DO jk = 2, nbld(ji,jj) 
    22892345                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
    2290                            pdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
     2346                           pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
    22912347                        END DO 
    22922348                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     
    22952351                        DO jk = 2, nbld(ji,jj) 
    22962352                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
    2297                            pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
     2353                           pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
    22982354                        END DO 
    22992355                     END IF   ! IF (shol >=0.5) 
     
    23352391      !! 
    23362392      !!---------------------------------------------------------------------- 
    2337       INTEGER,                     INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
    2338       REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
    2339       REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pviscos        ! Viscosity 
    2340       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    2341       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
    2342       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
    2343       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2344       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pshear         ! Shear production 
    2345       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    2346       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_min 
     2393      INTEGER,                          INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
     2394      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
     2395      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pviscos        ! Viscosity 
     2396      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
     2397      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
     2398      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2399      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2400      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production 
     2401      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     2402      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_min 
    23472403      ! 
    23482404      ! Local variables 
     
    23762432            zvel_sc_ml  = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird 
    23772433            zstab_fac   = ( phml(ji,jj) / zvel_sc_ml *   & 
    2378                &            ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10(-shol(ji,jj) ) ) )**1.25_wp ) )**2 
     2434               &            ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.25_wp ) )**2 
    23792435            ! 
    23802436            zdifml_sc(ji,jj) = pp_dif_ml * phml(ji,jj) * zvel_sc_ml 
     
    24972553#endif 
    24982554         ELSE 
    2499             zdifml_sc(ji,jj) = svstr(ji,jj) * phbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     2555            zdifml_sc(ji,jj) = svstr(ji,jj) * phbl(ji,jj) * MAX( EXP ( -1.0_wp * ( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
    25002556            zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 
    25012557#ifdef key_osm_debug 
     
    25092565      ! 
    25102566      DO_2D( 0, 0, 0, 0 ) 
    2511       IF ( l_conv(ji,jj) ) THEN 
    2512          DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity 
    2513             zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
    2514             pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
    2515             pviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_v_sc(ji,jj) * zznd_ml ) *   & 
    2516                &                ( 1.0_wp - 0.5_wp * zznd_ml**2 ) 
    2517          END DO 
     2567         IF ( l_conv(ji,jj) ) THEN 
     2568            DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity 
     2569               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
     2570               pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
     2571               pviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_v_sc(ji,jj) * zznd_ml ) *   & 
     2572                  &                ( 1.0_wp - 0.5_wp * zznd_ml**2 ) 
     2573            END DO 
     2574            ! 
     2575            ! Coupling to bottom 
     2576            ! 
     2577            IF ( l_coup(ji,jj) ) THEN                                                         ! ag 19/03 
     2578               DO jk = mbkt(ji,jj), nmld(ji,jj), -1                                           ! ag 19/03 
     2579                  zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )   ! ag 19/03 
     2580                  pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03 
     2581                  pdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2    ! ag 19/03 
     2582               END DO                                                                         ! ag 19/03 
     2583            ENDIF                                                                             ! ag 19/03 
     2584            ! Pycnocline 
     2585            IF ( l_pyc(ji,jj) ) THEN  
     2586               ! Diffusivity and viscosity profiles in the pycnocline given by 
     2587               ! cubic polynomial. Note, if l_pyc TRUE can't be coupled to seabed. 
     2588               za_cubic = 0.5_wp 
     2589               zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
     2590               zd_d_cubic = ( pdh(ji,jj) * zdifml_sc(ji,jj) / phml(ji,jj) * SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) *   & 
     2591                  &           ( 2.5_wp * zbeta_d_sc(ji,jj) - 1.0_wp ) - 0.85_wp * zdifpyc_s_sc(ji,jj) ) /            & 
     2592                  &           MAX( zdifpyc_n_sc(ji,jj), 1.0e-8_wp ) 
     2593               zd_d_cubic = zd_d_cubic - zb_d_cubic - 2.0_wp * ( 1.0_wp - za_cubic  - zb_d_cubic ) 
     2594               zc_d_cubic = 1.0_wp - za_cubic - zb_d_cubic - zd_d_cubic 
     2595               zb_v_cubic = -1.75_wp * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
     2596               zd_v_cubic = ( 0.5_wp * zvisml_sc(ji,jj) * pdh(ji,jj) / phml(ji,jj) - 0.85_wp * zvispyc_s_sc(ji,jj) ) /   & 
     2597                  &           MAX( zvispyc_n_sc(ji,jj), 1.0e-8_wp ) 
     2598               zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic ) 
     2599               zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic 
     2600               DO jk = nmld(ji,jj) , nbld(ji,jj) 
     2601                  zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp ) 
     2602                  ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ) 
     2603                  ! 
     2604                  pdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) *   & 
     2605                     &                ( za_cubic + zb_d_cubic * zznd_pyc + zc_d_cubic * zznd_pyc**2 + zd_d_cubic * zznd_pyc**3 ) 
     2606                  ! 
     2607                  pdiffut(ji,jj,jk) = pdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ztmp 
     2608                  pviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) *   & 
     2609                     &                ( za_cubic + zb_v_cubic * zznd_pyc + zc_v_cubic * zznd_pyc**2 + zd_v_cubic * zznd_pyc**3 ) 
     2610                  pviscos(ji,jj,jk) = pviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ztmp 
     2611               END DO 
     2612   !                  IF ( pdhdt(ji,jj) > 0._wp ) THEN 
     2613   !                     zdiffut(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
     2614   !                     zviscos(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
     2615   !                  ELSE 
     2616   !                     zdiffut(ji,jj,nbld(ji,jj)) = 0._wp 
     2617   !                     zviscos(ji,jj,nbld(ji,jj)) = 0._wp 
     2618   !                  ENDIF 
     2619            ENDIF 
     2620         ELSE 
     2621            ! Stable conditions 
     2622            DO jk = 2, nbld(ji,jj) 
     2623               zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
     2624               pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp 
     2625               pviscos(ji,jj,jk) = 0.375_wp * zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml ) * ( 1.0_wp - zznd_ml**2 ) 
     2626            END DO 
     2627            ! 
     2628            IF ( pdhdt(ji,jj) > 0.0_wp ) THEN 
     2629               pdiffut(ji,jj,nbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, nbld(ji,jj), Kmm) 
     2630               pviscos(ji,jj,nbld(ji,jj)) = pdiffut(ji,jj,nbld(ji,jj)) 
     2631            ENDIF 
     2632         ENDIF   ! End if ( l_conv ) 
    25182633         ! 
    2519          ! Coupling to bottom 
    2520          ! 
    2521          IF ( l_coup(ji,jj) ) THEN                                                         ! ag 19/03 
    2522             DO jk = mbkt(ji,jj), nmld(ji,jj), -1                                           ! ag 19/03 
    2523                zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )   ! ag 19/03 
    2524                pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03 
    2525                pdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2    ! ag 19/03 
    2526             END DO                                                                         ! ag 19/03 
    2527          ENDIF                                                                             ! ag 19/03 
    2528          ! Pycnocline 
    2529          IF ( l_pyc(ji,jj) ) THEN  
    2530             ! Diffusivity and viscosity profiles in the pycnocline given by 
    2531             ! cubic polynomial. Note, if l_pyc TRUE can't be coupled to seabed. 
    2532             za_cubic = 0.5_wp 
    2533             zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
    2534             zd_d_cubic = ( pdh(ji,jj) * zdifml_sc(ji,jj) / phml(ji,jj) * SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) *   & 
    2535                &           ( 2.5_wp * zbeta_d_sc(ji,jj) - 1.0_wp ) - 0.85_wp * zdifpyc_s_sc(ji,jj) ) /            & 
    2536                &           MAX( zdifpyc_n_sc(ji,jj), 1.0e-8_wp ) 
    2537             zd_d_cubic = zd_d_cubic - zb_d_cubic - 2.0_wp * ( 1.0_wp - za_cubic  - zb_d_cubic ) 
    2538             zc_d_cubic = 1.0_wp - za_cubic - zb_d_cubic - zd_d_cubic 
    2539             zb_v_cubic = -1.75_wp * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
    2540             zd_v_cubic = ( 0.5_wp * zvisml_sc(ji,jj) * pdh(ji,jj) / phml(ji,jj) - 0.85_wp * zvispyc_s_sc(ji,jj) ) /   & 
    2541                &           MAX( zvispyc_n_sc(ji,jj), 1.0e-8_wp ) 
    2542             zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic ) 
    2543             zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic 
    2544             DO jk = nmld(ji,jj) , nbld(ji,jj) 
    2545                zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp ) 
    2546                ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ) 
    2547                ! 
    2548                pdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) *   & 
    2549                   &                ( za_cubic + zb_d_cubic * zznd_pyc + zc_d_cubic * zznd_pyc**2 + zd_d_cubic * zznd_pyc**3 ) 
    2550                ! 
    2551                pdiffut(ji,jj,jk) = pdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ztmp 
    2552                pviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) *   & 
    2553                   &                ( za_cubic + zb_v_cubic * zznd_pyc + zc_v_cubic * zznd_pyc**2 + zd_v_cubic * zznd_pyc**3 ) 
    2554                pviscos(ji,jj,jk) = pviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ztmp 
    2555             END DO 
    2556 !                  IF ( pdhdt(ji,jj) > 0._wp ) THEN 
    2557 !                     zdiffut(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
    2558 !                     zviscos(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 ) 
    2559 !                  ELSE 
    2560 !                     zdiffut(ji,jj,nbld(ji,jj)) = 0._wp 
    2561 !                     zviscos(ji,jj,nbld(ji,jj)) = 0._wp 
    2562 !                  ENDIF 
    2563          ENDIF 
    2564       ELSE 
    2565          ! Stable conditions 
    2566          DO jk = 2, nbld(ji,jj) 
    2567             zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
    2568             pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp 
    2569             pviscos(ji,jj,jk) = 0.375_wp * zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml ) * ( 1.0_wp - zznd_ml**2 ) 
    2570          END DO 
    2571          ! 
    2572          IF ( pdhdt(ji,jj) > 0.0_wp ) THEN 
    2573             pdiffut(ji,jj,nbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, nbld(ji,jj), Kmm) 
    2574             pviscos(ji,jj,nbld(ji,jj)) = pdiffut(ji,jj,nbld(ji,jj)) 
    2575          ENDIF 
    2576       ENDIF   ! End if ( l_conv ) 
    2577       ! 
    25782634      END_2D 
    25792635      IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) )   ! BBL-coupling velocity scale 
     
    25822638   END SUBROUTINE zdf_osm_diffusivity_viscosity 
    25832639 
    2584    SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,                           & 
    2585       &                          pdhdt, pshear, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc,   & 
    2586       &                          palpha_pyc, pdiffut, pviscos ) 
     2640   SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,                              & 
     2641      &                          pdhdt, pshear, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext,   & 
     2642      &                          pdiffut, pviscos ) 
    25872643      !!--------------------------------------------------------------------- 
    25882644      !!                 ***  ROUTINE zdf_osm_fgr_terms *** 
     
    25932649      !! 
    25942650      !!---------------------------------------------------------------------- 
    2595       INTEGER,                     INTENT(in   ) ::   Kmm            ! Time-level index 
    2596       INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kp_ext         ! Offset for external level 
    2597       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
    2598       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
    2599       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
    2600       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2601       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pshear         ! Shear production 
    2602       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
    2603       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
    2604       REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pdbdz_pyc      ! Pycnocline buoyancy gradients 
    2605       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   palpha_pyc     ! 
    2606       REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pdiffut        ! t-diffusivity 
    2607       REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pviscos        ! Viscosity 
    2608       ! 
     2651      INTEGER,                          INTENT(in   ) ::   Kmm            ! Time-level index 
     2652      INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! Offset for external level 
     2653      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
     2654      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
     2655      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2656      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2657      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production 
     2658      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
     2659      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
     2660      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     2661      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity 
     2662      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity 
     2663      ! 
     2664      REAL(wp), DIMENSION(A2D(0))     ::   zalpha_pyc   ! 
     2665      REAL(wp), DIMENSION(A2D(0),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline 
    26092666      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles 
    26102667      ! 
     
    26392696      ! 
    26402697      IF( ln_timing ) CALL timing_start('zdf_osm_ft') 
     2698      ! 
     2699      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     2700      !  Pycnocline gradients for scalars and velocity 
     2701      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     2702      CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, zdbdz_pyc, zalpha_pyc, pdh,    & 
     2703         &                                       phbl, pdbdz_bl_ext, phml, pdhdt ) 
     2704#ifdef key_osm_debug 
     2705      IF(narea==nn_narea_db) THEN 
     2706         ji=iloc_db; jj=jloc_db 
     2707         jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) 
     2708         WRITE(narea+100,'(a,l7,/,3(a,g11.3),/)') & 
     2709            & 'After pycnocline profiles BL  lpyc=', l_pyc(ji,jj),& 
     2710            & 'sub-BL strat: zdtdz_bl_ext=', pdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', pdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', pdbdz_bl_ext(ji,jj), & 
     2711            & 'Pycnocline: zalpha_pyc=', zalpha_pyc(ji,jj) 
     2712         !       WRITE(narea+100,'(a,*(g11.3))') ' zdtdz_pyc[imld-1..ibld+2] =', ( zdtdz_pyc(ji,jj,jk), jk=jl,jm ) 
     2713         !       WRITE(narea+100,'(a,*(g11.3))') ' zdsdz_pyc[imld-1..ibld+2] =', ( zdsdz_pyc(ji,jj,jk), jk=jl,jm ) 
     2714         WRITE(narea+100,'(a,*(g11.3))') ' zdbdz_pyc[imld-1..ibld+2] =', ( zdbdz_pyc(ji,jj,jk), jk=jl,jm ) 
     2715         !       WRITE(narea+100,'(a,*(g11.3))') ' zdudz_pyc[imld-1..ibld+2] =', ( zdudz_pyc(ji,jj,jk), jk=jl,jm ) 
     2716         !       WRITE(narea+100,'(a,*(g11.3))') ' zdvdz_pyc[imld-1..ibld+2] =', ( zdvdz_pyc(ji,jj,jk), jk=jl,jm ) 
     2717         WRITE(narea+100,*) 
     2718         FLUSH(narea+100) 
     2719      END IF 
     2720#endif 
    26412721      ! 
    26422722      ! Auxiliary indices 
     
    27832863               zdelta_pyc          = ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird /   & 
    27842864                  &                       SQRT( MAX( zbuoy_pyc_sc, ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 
    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) ) *   & 
     2865               zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( zalpha_pyc(ji,jj) * av_dt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   & 
    27862866                  &                     zdelta_pyc**2 / pdh(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) ) *   & 
     2867               zws_pyc_sc_1(ji,jj) = 0.325_wp * ( zalpha_pyc(ji,jj) * av_ds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   & 
    27882868                  &                     zdelta_pyc**2 / pdh(ji,jj) 
    27892869               zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 
     
    28032883            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
    28042884            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) -                                                                                & 
    2805                &              0.045_wp * ( ( zwth_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
     2885               &              0.045_wp * ( ( zwth_ent(ji,jj) * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
    28062886               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 
    28072887            ghams(ji,jj,jk) = ghams(ji,jj,jk) -                                                                                & 
    2808                &              0.045_wp * ( ( zws_ent(ji,jj)  * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
     2888               &              0.045_wp * ( ( zws_ent(ji,jj)  * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
    28092889               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 
    28102890#ifdef key_osm_debug 
     
    28352915      END_3D 
    28362916      ! 
    2837       IF(ln_dia_osm) THEN 
     2917      IF ( ln_dia_osm ) THEN 
    28382918         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! Upward turb. temperature entrainment flux 
    28392919         IF ( iom_use("zws_ent")  ) CALL iom_put( "zws_ent",  tmask(:,:,1)*zws_ent  )   ! Upward turb. salinity    entrainment flux 
     
    28902970            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) *                 & 
    28912971               &                                ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) *   & 
    2892                &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
     2972               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    28932973            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse(ji,jj) *                 & 
    28942974               &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   & 
    2895                &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
     2975               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    28962976         END IF   ! l_conv .AND. l_pyc 
    28972977      END_3D 
     
    29172997#endif 
    29182998 
    2919       IF(ln_dia_osm) THEN 
     2999      IF ( ln_dia_osm ) THEN 
    29203000         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0",    wmask*ghamu           ) 
    29213001         IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     
    30273107#endif 
    30283108      ! 
    3029       IF(ln_dia_osm) THEN 
     3109      IF ( ln_dia_osm ) THEN 
    30303110         IF ( iom_use("ghamu_f") )    CALL iom_put( "ghamu_f",    wmask       *ghamu    ) 
    30313111         IF ( iom_use("ghamv_f") )    CALL iom_put( "ghamv_f",    wmask       *ghamv    ) 
     
    31573237         IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
    31583238      END IF 
    3159       IF(ln_dia_osm) THEN 
     3239      IF ( ln_dia_osm ) THEN 
    31603240         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
    31613241         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     
    31853265#endif 
    31863266      ! 
    3187       IF(ln_dia_osm) THEN 
    3188          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu   ) 
    3189          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv   ) 
    3190          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*pviscos ) 
     3267      IF ( ln_dia_osm ) THEN 
     3268         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     3269         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     3270         IF ( iom_use("zviscos") ) THEN 
     3271            osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pviscos; CALL iom_put( "zviscos", osmdia3d ) 
     3272         END IF 
    31913273      END IF 
    31923274      ! 
     
    32093291      !! 
    32103292      !!---------------------------------------------------------------------- 
    3211       INTEGER,                  INTENT(in   ) ::   Kmm          ! Time-level index 
    3212       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == ! 
    3213       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization 
    3214       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization 
    3215       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization 
    3216       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization 
    3217       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdbdx_mle    ! MLE horiz gradients at u points 
    3218       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdbdy_mle    ! MLE horiz gradients at v points 
    3219       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient 
     3293      INTEGER,                      INTENT(in   ) ::   Kmm          ! Time-level index 
     3294      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == ! 
     3295      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization 
     3296      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization 
     3297      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization 
     3298      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization 
     3299      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdx_mle    ! MLE horiz gradients at u points 
     3300      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdy_mle    ! MLE horiz gradients at v points 
     3301      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient 
    32203302      ! 
    32213303      ! Local variables 
     
    33163398      !!----------------------------------------------------------------------       
    33173399      ! Outputs 
    3318       INTEGER,                  INTENT(in   ) ::   Kmm         ! Time-level index 
    3319       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pwb_fk 
    3320       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phbl        ! BL depth 
    3321       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phmle       ! MLE depth 
    3322       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux 
    3323       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
     3400      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index 
     3401      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pwb_fk 
     3402      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth 
     3403      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phmle       ! MLE depth 
     3404      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux 
     3405      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    33243406      ! 
    33253407      ! Local variables 
     
    34503532      INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kmld_prof 
    34513533      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == ! 
    3452       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phmle       ! MLE depth 
    3453       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK 
    3454       REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff 
    3455       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    3456       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl        ! BL depth 
     3534      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phmle       ! MLE depth 
     3535      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK 
     3536      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff 
     3537      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
     3538      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth 
    34573539      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation 
    34583540      ! 
     
    36703752         ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
    36713753         z1_t2 = 2e-5_wp 
    3672          DO_2D( 1, 1, 1, 1 ) 
     3754         DO_2D( 0, 0, 0, 0 ) 
    36733755            r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 ) 
    36743756         END_2D 
     
    37483830      ghamv(:,:,:) = 0.0_wp 
    37493831      ! 
     3832      IF ( ln_dia_osm ) osmdia2d(:,:) = 0.0_wp   ! Initialise auxiliary array for diagnostic output 
     3833      ! 
    37503834      IF( ln_timing ) CALL timing_stop('zdf_osm_init') 
    37513835      ! 
Note: See TracChangeset for help on using the changeset viewer.