Changeset 12216


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

Alan's changes

File:
1 edited

Legend:

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

    r12178 r12216  
    2525   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 
    2626   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 
    27    !!            (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out). 
     27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 
    2828   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 
    2929   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 
    3030   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 
     31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, 
     32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers 
     33   !!                  (b) The stable OSBL depth parametrization changed. 
     34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. 
     35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 
    3136   !!---------------------------------------------------------------------- 
    3237 
     
    4045   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD) 
    4146   !!   dyn_osm       : compute and add to u & v trensd the non-local flux 
     47   !! 
     48   !! Subroutines in revised code. 
    4249   !!---------------------------------------------------------------------- 
    4350   USE oce            ! ocean dynamics and active tracers 
     
    97104 
    98105   !                                    !!! ** General constants  ** 
    99    REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number 
     106   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero 
     107   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw 
    100108   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3 
    101109   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3 
     
    166174      REAL(wp) ::   zbeta, zthermal                                  ! 
    167175      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 
    168       REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       ! 
     176      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      ! 
     177 
    169178      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density 
    170179      INTEGER  ::   jm                          ! dummy loop indices 
     
    196205      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 
    197206      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer 
    198       LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 
     207      LOGICAL, DIMENSION(jpi,jpj)  :: lconv    ! unstable/stable bl 
    199208 
    200209      ! mixed-layer variables 
     
    238247      ! Temporary variables 
    239248      INTEGER :: inhml 
    240       INTEGER :: i_lconv_alloc 
    241249      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 
    242250      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables 
     
    248256      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    249257 
     258      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
     259      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 
     260      REAL(wp) :: zzeta_s = 0._wp 
     261      REAL(wp) :: zzeta_v = 0.46 
     262      REAL(wp) :: zabsstke 
     263 
    250264      ! For debugging 
    251265      INTEGER :: ikt 
    252266      !!-------------------------------------------------------------------- 
    253267      ! 
    254       ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc ) 
    255       IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv') 
    256  
    257268      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    258269      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
     
    268279      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp 
    269280      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
     281 
    270282      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    271283      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
     
    287299      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
    288300 
     301      zdhdt_2(:,:) = 0._wp 
    289302      ! hbl = MAX(hbl,epsln) 
    290303      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    350363              ! Use wind speed wndm included in sbc_oce module 
    351364              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    352               dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
     365              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    353366           END DO 
    354367        END DO 
     
    362375              ! It could represent the effects of the spread of wave directions 
    363376              ! around the mean wind. The effect of this adjustment needs to be tested. 
    364               zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
    365                    &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
    366               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
     377              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     378              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     379              dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
    367380           END DO 
    368381        END DO 
     
    375388           ! Langmuir velocity scale (zwstrl), at T-point 
    376389           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
    377            ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
    378            ! Intended as a possible test. Doesn't affect LES results for entrainment, 
    379            !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
    380            zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
    381                 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
    382                 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
    383            ! define La this way so effects of Stokes penetration depth on velocity scale are included 
    384            zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
     390           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 
     391           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 
    385392           ! Velocity scale that tends to zustar for large Langmuir numbers 
    386393           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
     
    389396           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    390397           ! Note zustke and zwstrl are not amended. 
    391            IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    392398           ! 
    393399           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
     
    406412     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 
    407413     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408      ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
    410       ibld(:,:) = 3 
    411       DO jk = 4, jpkm1 
     414     ! BL must be always 4 levels deep. 
     415      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
     416      ibld(:,:) = 4 
     417      DO jk = 5, jpkm1 
    412418         DO jj = 2, jpjm1 
    413419            DO ji = 2, jpim1 
     
    419425      END DO 
    420426 
    421       DO jj = 2, jpjm1                                 !  Vertical slab 
     427      DO jj = 2, jpjm1 
    422428         DO ji = 2, jpim1 
    423                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    424                zbeta    = rab_n(ji,jj,1,jp_sal) 
    425                zt   = 0._wp 
    426                zs   = 0._wp 
    427                zu   = 0._wp 
    428                zv   = 0._wp 
    429                ! average over depth of boundary layer 
    430                zthick=0._wp 
    431                DO jm = 2, ibld(ji,jj) 
    432                   zthick=zthick+e3t_n(ji,jj,jm) 
    433                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    434                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    435                   zu   = zu  + e3t_n(ji,jj,jm) & 
    436                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    437                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    438                   zv   = zv  + e3t_n(ji,jj,jm) & 
    439                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    440                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    441                END DO 
    442                zt_bl(ji,jj) = zt / zthick 
    443                zs_bl(ji,jj) = zs / zthick 
    444                zu_bl(ji,jj) = zu / zthick 
    445                zv_bl(ji,jj) = zv / zthick 
    446                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    447                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    448                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    449                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    450                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    451                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    452                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    453                IF ( lconv(ji,jj) ) THEN    ! Convective 
    454                       zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    455                            &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    456  
    457                       zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    458                            &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    460 !                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    461 !                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    462  
    463 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    464 !                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    465                       zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    466                ELSE                        ! Stable 
    467                       zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
    468                            &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
    469                            & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
    470                            &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
    471                       zzdhdt = zzdhdt + zwbav(ji,jj) 
    472                       IF ( zzdhdt < 0._wp ) THEN 
    473                       ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    474                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    475                       ELSE 
    476                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    477                               &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    478                       ENDIF 
    479                       zzdhdt = 2.0 * zzdhdt / zpert 
    480                ENDIF 
    481                zdhdt(ji,jj) = zzdhdt 
    482            END DO 
     429            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     430            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 
     431            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     432         END DO 
    483433      END DO 
    484434 
     
    495445            DO ji = 2, jpim1 
    496446               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    497                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     447                  ibld(ji,jj) = jk 
    498448               ENDIF 
    499449            END DO 
     
    504454! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    505455! 
    506       DO jj = 2, jpjm1 
    507          DO ji = 2, jpim1 
    508             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     456      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 
     457      ! Alan: do we need zb_ml? 
     458      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    509459! 
    510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    511460! 
    512                zhbl_s = hbl(ji,jj) 
    513                jm = imld(ji,jj) 
    514                zthermal = rab_n(ji,jj,1,jp_tem) 
    515                zbeta = rab_n(ji,jj,1,jp_sal) 
    516                IF ( lconv(ji,jj) ) THEN 
    517 !unstable 
    518                   zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    519                        &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    520  
    521                   DO jk = imld(ji,jj), ibld(ji,jj) 
    522                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
    523                           & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524  
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
    526                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527  
    528                      IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
    529                   END DO 
    530                   hbl(ji,jj) = zhbl_s 
    531                   ibld(ji,jj) = jm 
    532                   hbli(ji,jj) = hbl(ji,jj) 
    533                ELSE 
    534 ! stable 
    535                   DO jk = imld(ji,jj), ibld(ji,jj) 
    536                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          & 
    537                           &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 
    538                           & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    539  
    540                      zhbl_s = zhbl_s +  (                                                                                & 
    541                           &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
    542                           &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
    543                           &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    544                           &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    545                           &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here 
    546  
    547                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    548                      IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
    549                   END DO 
    550                   hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 
    551                   ibld(ji,jj) = MAX(jm, 3 ) 
    552                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    553                ENDIF   ! IF ( lconv ) 
    554             ELSE 
    555 ! change zero or one model level. 
    556                hbl(ji,jj) = zhbl_t(ji,jj) 
    557                IF ( lconv(ji,jj) ) THEN 
    558                   hbli(ji,jj) = hbl(ji,jj) 
    559                ELSE 
    560                   hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 
    561                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    562                ENDIF 
    563             ENDIF 
    564             zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    565          END DO 
    566       END DO 
     461      CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
    567462      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    568  
    569 ! Recalculate averages over boundary layer after depth updated 
    570      ! Consider later  combining this into the loop above and looking for columns 
    571      ! where the index for base of the boundary layer have changed 
    572       DO jj = 2, jpjm1                                 !  Vertical slab 
    573          DO ji = 2, jpim1 
    574                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    575                zbeta    = rab_n(ji,jj,1,jp_sal) 
    576                zt   = 0._wp 
    577                zs   = 0._wp 
    578                zu   = 0._wp 
    579                zv   = 0._wp 
    580                ! average over depth of boundary layer 
    581                zthick=0._wp 
    582                DO jm = 2, ibld(ji,jj) 
    583                   zthick=zthick+e3t_n(ji,jj,jm) 
    584                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    585                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    586                   zu   = zu  + e3t_n(ji,jj,jm) & 
    587                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    588                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    589                   zv   = zv  + e3t_n(ji,jj,jm) & 
    590                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    591                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    592                END DO 
    593                zt_bl(ji,jj) = zt / zthick 
    594                zs_bl(ji,jj) = zs / zthick 
    595                zu_bl(ji,jj) = zu / zthick 
    596                zv_bl(ji,jj) = zv / zthick 
    597                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    598                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    599                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    600                       &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    601                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    602                       &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    603                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    604                zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    605                IF ( lconv(ji,jj) ) THEN 
    606                   IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    607                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    608                            zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
    609                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    610                      ELSE                                                     ! unstable 
    611                            zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
    612                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    613                      ENDIF 
    614                      IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
    615                         zari = 0.2 
    616                         zwb_ent(ji,jj) = 0._wp 
    617                      ENDIF 
    618                      inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    619                      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    620                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    621                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    622                   ELSE  ! IF (zdb_bl) 
    623                      imld(ji,jj) = ibld(ji,jj) - 1 
    624                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    625                      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    626                   ENDIF 
    627                ELSE   ! IF (lconv) 
    628                   IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    629                   ! boundary layer deepening 
    630                      IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    631                   ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    632                         zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    633                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    634                         inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    635                         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    636                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    637                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    638                      ELSE 
    639                         imld(ji,jj) = ibld(ji,jj) - 1 
    640                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    641                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    642                      ENDIF ! IF (zdb_bl > 0.0) 
    643                   ELSE     ! IF(dhdt >= 0) 
    644                   ! boundary layer collapsing. 
    645                      imld(ji,jj) = ibld(ji,jj) 
    646                      zhml(ji,jj) = zhbl(ji,jj) 
    647                      zdh(ji,jj) = 0._wp 
    648                   ENDIF    ! IF (dhdt >= 0) 
    649                ENDIF       ! IF (lconv) 
    650          END DO 
    651       END DO 
    652  
    653       ! Average over the depth of the mixed layer in the convective boundary layer 
    654       ! Also calculate entrainment fluxes for temperature and salinity 
    655       DO jj = 2, jpjm1                                 !  Vertical slab 
    656          DO ji = 2, jpim1 
    657             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    658             zbeta    = rab_n(ji,jj,1,jp_sal) 
    659             IF ( lconv(ji,jj) ) THEN 
    660                zt   = 0._wp 
    661                zs   = 0._wp 
    662                zu   = 0._wp 
    663                zv   = 0._wp 
    664                ! average over depth of boundary layer 
    665                zthick=0._wp 
    666                DO jm = 2, imld(ji,jj) 
    667                   zthick=zthick+e3t_n(ji,jj,jm) 
    668                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    669                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    670                   zu   = zu  + e3t_n(ji,jj,jm) & 
    671                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    672                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    673                   zv   = zv  + e3t_n(ji,jj,jm) & 
    674                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    675                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    676                END DO 
    677                zt_ml(ji,jj) = zt / zthick 
    678                zs_ml(ji,jj) = zs / zthick 
    679                zu_ml(ji,jj) = zu / zthick 
    680                zv_ml(ji,jj) = zv / zthick 
    681                zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    682                zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    683                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    684                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    685                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    686                      &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    687                zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    688             ELSE 
    689             ! stable, if entraining calulate average below interface layer. 
    690                IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    691                   zt   = 0._wp 
    692                   zs   = 0._wp 
    693                   zu   = 0._wp 
    694                   zv   = 0._wp 
    695                   ! average over depth of boundary layer 
    696                   zthick=0._wp 
    697                   DO jm = 2, imld(ji,jj) 
    698                      zthick=zthick+e3t_n(ji,jj,jm) 
    699                      zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    700                      zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    701                      zu   = zu  + e3t_n(ji,jj,jm) & 
    702                         &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    703                         &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    704                      zv   = zv  + e3t_n(ji,jj,jm) & 
    705                         &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    706                         &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    707                   END DO 
    708                   zt_ml(ji,jj) = zt / zthick 
    709                   zs_ml(ji,jj) = zs / zthick 
    710                   zu_ml(ji,jj) = zu / zthick 
    711                   zv_ml(ji,jj) = zv / zthick 
    712                   zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    713                   zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    714                   zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    715                         &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    716                   zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    717                         &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    718                   zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    719                ENDIF 
    720             ENDIF 
    721          END DO 
    722       END DO 
    723     ! 
     463! 
     464    ! Average over the depth of the mixed layer in the convective boundary layer 
     465    ! Alan: do we need zb_ml? 
     466    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    724467    ! rotate mean currents and changes onto wind align co-ordinates 
    725468    ! 
    726  
    727       DO jj = 2, jpjm1 
    728          DO ji = 2, jpim1 
    729             ztemp = zu_ml(ji,jj) 
    730             zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
    731             zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    732             ztemp = zdu_ml(ji,jj) 
    733             zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
    734             zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    735     ! 
    736             ztemp = zu_bl(ji,jj) 
    737             zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
    738             zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    739             ztemp = zdu_bl(ji,jj) 
    740             zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
    741             zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    742          END DO 
    743       END DO 
    744  
     469    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
     470    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    745471     zuw_bse = 0._wp 
    746472     zvw_bse = 0._wp 
     473     zwth_ent = 0._wp 
     474     zws_ent = 0._wp 
    747475     DO jj = 2, jpjm1 
    748476        DO ji = 2, jpim1 
    749  
    750            IF ( lconv(ji,jj) ) THEN 
    751               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    752                  zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    753                  zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     477           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 
     478              IF ( lconv(ji,jj) ) THEN 
     479             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 
     480                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 
     481                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 
     482            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 
     483                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 
     484                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 
     485                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     486                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     487                 ENDIF 
     488              ELSE 
     489                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
     490                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    754491              ENDIF 
    755            ELSE 
    756               zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    757               zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
    758492           ENDIF 
    759493        END DO 
     
    764498      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    765499 
    766        DO jj = 2, jpjm1 
    767           DO ji = 2, jpim1 
    768           ! 
    769              IF ( lconv (ji,jj) ) THEN 
    770              ! Unstable conditions 
    771                 IF( zdb_bl(ji,jj) > 0._wp ) THEN 
    772                 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
    773                    ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
    774                    zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
    775                    zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    776                    DO jk = 2 , ibld(ji,jj) 
    777                       znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    778                       zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    779                       zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    780                       zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    781                    END DO 
    782                 ENDIF 
    783              ELSE 
    784              ! stable conditions 
    785              ! if pycnocline profile only defined when depth steady of increasing. 
    786                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
    787                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    788                      IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    789                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
    790                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
    791                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    792                          DO jk = 2, ibld(ji,jj) 
    793                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    794                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    795                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    796                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    797                          END DO 
    798                      ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    799                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
    800                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
    801                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    802                          DO jk = 2, ibld(ji,jj) 
    803                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    804                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    805                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    806                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    807                          END DO 
    808                       ENDIF ! IF (zhol >=0.5) 
    809                    ENDIF    ! IF (zdb_bl> 0.) 
    810                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
    811              ENDIF          ! IF (lconv) 
    812             ! 
    813           END DO 
    814        END DO 
    815 ! 
    816        DO jj = 2, jpjm1 
    817           DO ji = 2, jpim1 
    818           ! 
    819              IF ( lconv (ji,jj) ) THEN 
    820              ! Unstable conditions 
    821                  zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
    822                & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
    823                 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
    824               & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    825                 DO jk = 2 , ibld(ji,jj)-1 
    826                    znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    827                    zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    828                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    829                 END DO 
    830              ELSE 
    831              ! stable conditions 
    832                 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    833                 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    834                 DO jk = 2, ibld(ji,jj) 
    835                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    836                    IF ( znd < 1.0 ) THEN 
    837                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    838                    ELSE 
    839                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    840                    ENDIF 
    841                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
    842                 END DO 
    843              ENDIF 
    844             ! 
    845           END DO 
    846        END DO 
     500      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 
     501      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 
     502      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    847503       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    848504       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    849505       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    850506 
    851       ! WHERE ( lconv ) 
    852       !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 
    853       !     zvisml_sc = zdifml_sc 
    854       !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 
    855       !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 
    856       !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 
    857       !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 
    858       !  ELSEWHERE 
    859       !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    860       !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    861       !  ENDWHERE 
    862507       DO jj = 2, jpjm1 
    863508          DO ji = 2, jpim1 
     
    868513               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    869514               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    870                zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
     515               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 
    871516             ELSE 
    872517               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    873518               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    874             END IF 
    875         END DO 
    876     END DO 
     519             END IF 
     520          END DO 
     521       END DO 
    877522! 
    878523       DO jj = 2, jpjm1 
     
    889534                ! pycnocline - if present linear profile 
    890535                IF ( zdh(ji,jj) > 0._wp ) THEN 
     536                   zgamma_b = 6.0 
    891537                   DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892538                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    893539                       ! 
    894                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     540                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    895541                       ! 
    896                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
     542                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 
    897543                   END DO 
     544                   IF ( ibld_ext == 0 ) THEN 
     545                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     546                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     547                   ELSE 
     548                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     549                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 
     550                   ENDIF 
    898551                ENDIF 
    899552                ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
     
    908561                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    909562                END DO 
     563 
     564                IF ( ibld_ext == 0 ) THEN 
     565                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     566                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     567                ELSE 
     568                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     569                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 
     570                ENDIF 
    910571             ENDIF   ! end if ( lconv ) 
    911 ! 
     572             ! 
    912573          END DO  ! end of ji loop 
    913574       END DO     ! end of jj loop 
     
    952613       END DO     ! end of jj loop 
    953614 
    954  
    955615! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 
    956616       WHERE ( lconv ) 
    957           zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ) 
    958           zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln ) 
    959           zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
     617          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 ) 
     618          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
     619          zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 
    960620       ELSEWHERE 
    961621          zsc_uw_1 = zustar**2 
    962           zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 
     622          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    963623       ENDWHERE 
    964  
     624       IF(ln_dia_osm) THEN 
     625          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     626          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     627       END IF 
    965628       DO jj = 2, jpjm1 
    966629          DO ji = 2, jpim1 
     
    1007670                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1008671                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    1009                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     672                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
    1010673                   ! non-gradient buoyancy terms 
    1011674                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
     
    1020683          END DO   ! ji loop 
    1021684       END DO      ! jj loop 
    1022  
    1023685 
    1024686       WHERE ( lconv ) 
     
    1051713       END DO           ! jj loop 
    1052714 
     715       IF(ln_dia_osm) THEN 
     716          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
     717          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     718       END IF 
    1053719! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    1054720 
     
    1089755       END DO      ! jj loop 
    1090756 
    1091  
    1092757       WHERE ( lconv ) 
    1093758          zsc_uw_1 = zustar**2 
     
    1134799          END DO 
    1135800       END DO 
     801 
     802       IF(ln_dia_osm) THEN 
     803          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
     804          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
     805          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     806          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     807          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     808          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     809       END IF 
    1136810! 
    1137811! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
     
    1165839      END DO 
    1166840 
     841       IF(ln_dia_osm) THEN 
     842          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     843          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     844       END IF 
    1167845      ! pynocline contributions 
    1168846       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
     
    1170848       DO jj = 2, jpjm1 
    1171849          DO ji = 2, jpim1 
    1172              DO jk= 2, ibld(ji,jj) 
    1173                 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1174                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1175                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1176                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1177                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
    1178                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    1179              END DO 
    1180            END DO 
     850             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 
     851                DO jk= 2, ibld(ji,jj) 
     852                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     853                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     854                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     855                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     856                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
     857                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     858                END DO 
     859             END IF 
     860          END DO 
    1181861       END DO 
    1182862 
     
    1185865       DO jj=2, jpjm1 
    1186866          DO ji = 2, jpim1 
    1187              IF ( lconv(ji,jj) ) THEN 
     867            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 
    1188868               DO jk = 1, imld(ji,jj) - 1 
    1189869                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1190                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1191                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
     870                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
     871                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1192872                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1193873                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
     
    1195875               DO jk = imld(ji,jj), ibld(ji,jj) 
    1196876                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    1197                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1198                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
     877                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
     878                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1199879                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1200880                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1201881                END DO 
    1202882             ENDIF 
    1203              ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
    1204              ghams(ji,jj,ibld(ji,jj)) = 0._wp 
    1205              ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
    1206              ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
     883 
     884             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     885             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     886             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
     887             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1207888          END DO       ! ji loop 
    1208889       END DO          ! jj loop 
    1209890 
    1210  
     891       IF(ln_dia_osm) THEN 
     892          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
     893          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
     894          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 
     895          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 
     896          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
     897          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
     898          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
     899       END IF 
    1211900       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1212901       ! Need to put in code for contributions that are applied explicitly to 
     
    1287976 
    1288977       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
    1289        CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. ) 
     978       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 
    1290979 
    1291980       ! GN 25/8: need to change tmask --> wmask 
     
    13191008        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    13201009        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   & 
    1321          &                  ghamu, 'U', 1. , ghamv, 'V', 1. ) 
    1322  
    1323        IF(ln_dia_osm) THEN 
     1010         &                  ghamu, 'U', -1. , ghamv, 'V', -1. ) 
     1011 
     1012      IF(ln_dia_osm) THEN 
    13241013         SELECT CASE (nn_osm_wave) 
    13251014         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1). 
     
    13301019         ! Stokes drift read in from sbcwave  (=2). 
    13311020         CASE(2) 
    1332             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    1333             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
     1021            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
     1022            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     1023            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period 
     1024            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height 
     1025            IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) )                  ! wave mean period from NP spectrum 
     1026            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum 
     1027            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10 
    13341028            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
    13351029                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
     
    13421036         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0> 
    13431037         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth 
    1344          IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth 
     1038         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k 
     1039         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base 
     1040         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base 
     1041         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base 
     1042         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base 
     1043         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base 
     1044         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth 
     1045         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth 
    13451046         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth 
    13461047         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points 
     
    13481049         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    13491050         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
     1051         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale 
     1052         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir # 
    13501053         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    13511054         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    13521055         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    13531056         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    1354          IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine 
     1057         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine 
     1058         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine 
    13551059         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    13561060         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine 
     
    13781082     INTEGER  ::   ios            ! local integer 
    13791083     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     1084     REAL z1_t2 
    13801085     !! 
    13811086     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
     
    14011106        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
    14021107        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    1403         WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes 
     1108        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
    14041109        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave 
    14051110        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave 
     
    15361241     REAL(wp) ::   zN2_c           ! local scalar 
    15371242     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    1538      INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top) 
     1243     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 
    15391244     !!---------------------------------------------------------------------- 
    15401245     ! 
     
    15511256           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    15521257        END IF 
     1258 
    15531259        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    1554         id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. ) 
     1260        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. ) 
    15551261        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return 
    15561262           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 
    1557            CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios  ) 
    1558            WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file' 
     1263           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  ) 
     1264           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file' 
    15591265           RETURN 
    1560         ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate 
     1266        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate 
    15611267           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 
    15621268        END IF 
     
    15701276         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
    15711277         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    1572          CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
     1278         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios ) 
    15731279        RETURN 
    15741280     END IF 
     
    15781284     !!----------------------------------------------------------------------------- 
    15791285     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 
    1580      ALLOCATE( imld_rst(jpi,jpj) ) 
    15811286     ! w-level of the mixing and mixed layers 
    15821287     CALL eos_rab( tsn, rab_n ) 
     
    15991304     DO jj = 1, jpj 
    16001305        DO ji = 1, jpi 
    1601            iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
     1306           iiki = MAX(4,imld_rst(ji,jj)) 
     1307           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth 
     1308           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth 
    16031309        END DO 
    16041310     END DO 
     
    16071313     DEALLOCATE( imld_rst ) 
    16081314     WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 
     1315     wn(:,:,:) = 0._wp 
     1316     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
    16091317   END SUBROUTINE osm_rst 
    16101318 
     
    16341342      ENDIF 
    16351343 
    1636       ! add non-local temperature and salinity flux 
    16371344      DO jk = 1, jpkm1 
    16381345         DO jj = 2, jpjm1 
     
    16481355      END DO 
    16491356 
    1650  
    1651       ! save the non-local tracer flux trends for diagnostic 
     1357      ! save the non-local tracer flux trends for diagnostics 
    16521358      IF( l_trdtra )   THEN 
    16531359         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    16541360         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    1655 !!bug gm jpttdzdf ==> jpttosm 
    1656          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    1657          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     1361 
     1362         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 
     1363         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 
    16581364         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    16591365      ENDIF 
     
    17231429 
    17241430   !!====================================================================== 
     1431 
    17251432END MODULE zdfosm 
Note: See TracChangeset for help on using the changeset viewer.