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 14316 – NEMO

Changeset 14316


Ignore:
Timestamp:
2021-01-19T19:57:12+01:00 (2 years ago)
Author:
smueller
Message:

Upgrade of a section of subroutine zdf_osm to a streamlined module subroutine of module zdfosm (ticket #2353)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90

    r14305 r14316  
    277277      INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 
    278278 
    279       REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients 
    280       REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear 
    281  
    282279      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid 
    283280      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid 
     
    299296      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 
    300297!      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
    301       REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
    302       REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
    303       REAL(wp) ::   zdtdz_pyc    ! Parametrized gradient of temperature in pycnocline 
    304       REAL(wp) ::   zdsdz_pyc    ! Parametrised gradient of salinity in pycnocline 
    305298      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline 
    306       REAL(wp) ::   zdudz_pyc    ! u-shear across the pycnocline 
    307       REAL(wp) ::   zdvdz_pyc    ! v-shear across the pycnocline 
    308299      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
    309300      ! Flux-gradient relationship variables 
    310301      REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 
    311302 
    312       REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale. 
    313  
    314       REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline. 
    315       REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. 
    316       REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ 
    317       REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. 
    318303      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep 
    319304 
     
    334319      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 
    335320      REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 
    336       REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. 
    337       REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc 
    338       REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max 
    339321      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme 
    340322      REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 
     
    344326      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 
    345327      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 
    346       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles 
    347       INTEGER  ::   istat                                                     ! Memory allocation status 
    348328 
    349329      ! For debugging 
     
    354334      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    355335      zustar(:,:) = 0.0_wp 
    356       zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp 
     336      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp 
    357337      zwth0(:,:)  = 0.0_wp ; zws0(:,:)  = 0.0_wp ; zwb0(:,:)      = 0.0_wp 
    358338      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp 
     
    372352      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp 
    373353      zdb_ml(:,:)  = 0._wp 
    374       zwth_ent = 0._wp ; zws_ent = 0._wp 
    375354      ! 
    376355      zdbdz_pyc(:,:,:) = 0.0_wp 
     
    386365      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
    387366 
    388       ! Flux-Gradient arrays. 
    389       zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp 
    390       zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp 
    391       zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp 
     367      zhbl_t(:,:)     = 0._wp 
    392368 
    393369      zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp 
     
    706682       CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 
    707683 
    708        ! 
    709        ! calculate non-gradient components of the flux-gradient relationships 
    710        ! 
    711 ! Stokes term in scalar flux, flux-gradient relationship 
    712        WHERE ( lconv ) 
    713           zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln) 
    714           ! 
    715           zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    716        ELSEWHERE 
    717           zsc_wth_1 = 2.0 * zwthav 
    718           ! 
    719           zsc_ws_1 = 2.0 * zwsav 
    720        ENDWHERE 
    721  
    722  
    723        DO_2D( 0, 0, 0, 0 ) 
    724          IF ( lconv(ji,jj) ) THEN 
    725            DO jk = 2, imld(ji,jj) 
    726               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    727               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
    728               ! 
    729               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
    730            END DO ! end jk loop 
    731          ELSE     ! else for if (lconv) 
    732  ! Stable conditions 
    733             DO jk = 2, ibld(ji,jj) 
    734                zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    735                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
    736                     &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
    737                ! 
    738                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
    739                     &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
    740             END DO 
    741          ENDIF               ! endif for check on lconv 
    742  
    743        END_2D 
    744  
    745 ! 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) 
    746        WHERE ( lconv ) 
    747           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 ) 
    748           zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
    749           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 ) 
    750        ELSEWHERE 
    751           zsc_uw_1 = zustar**2 
    752           zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
    753        ENDWHERE 
    754        IF(ln_dia_osm) THEN 
    755           IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
    756           IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
    757        END IF 
    758        DO_2D( 0, 0, 0, 0 ) 
    759           IF ( lconv(ji,jj) ) THEN 
    760              DO jk = 2, imld(ji,jj) 
    761                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    762                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   & 
    763                      &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) & 
    764                      &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) 
    765 ! 
    766                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       & 
    767                      &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 
    768              END DO   ! end jk loop 
    769           ELSE 
    770 ! Stable conditions 
    771              DO jk = 2, ibld(ji,jj) ! corrected to ibld 
    772                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    773                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       & 
    774                      &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 
    775                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 
    776              END DO   ! end jk loop 
    777           ENDIF 
    778        END_2D 
    779  
    780 ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] 
    781  
    782        WHERE ( lconv ) 
    783           zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    784           zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    785        ELSEWHERE 
    786           zsc_wth_1 = 0._wp 
    787           zsc_ws_1 = 0._wp 
    788        ENDWHERE 
    789  
    790        DO_2D( 0, 0, 0, 0 ) 
    791           IF (lconv(ji,jj) ) THEN 
    792              DO jk = 2, imld(ji,jj) 
    793                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    794                 ! calculate turbulent length scale 
    795                 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    796                      &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) ) 
    797                 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    798                      &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    799                 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 
    800                 ! non-gradient buoyancy terms 
    801                 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 ) 
    802                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    803              END DO 
    804  
    805              IF ( lpyc(ji,jj) ) THEN 
    806                ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
    807                ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) 
    808                zwth_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) 
    809                zws_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 
    810 ! Cubic profile used for buoyancy term 
    811                za_cubic = 0.755 * ztau_sc_u(ji,jj) 
    812                zb_cubic = 0.25 * ztau_sc_u(ji,jj) 
    813                DO jk = 2, ibld(ji,jj) 
    814                  zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
    815                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
    816  
    817                  ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 
    818                END DO 
    819 ! 
    820                zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 
    821                zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) ) 
    822 ! 
    823                zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
    824 ! 
    825                zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 
    826 ! 
    827                zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
    828                DO jk = 2, ibld(ji,jj) 
    829                  zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
    830                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
    831 ! 
    832                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
    833                END DO 
    834             ENDIF ! End of pycnocline 
    835           ELSE ! lconv test - stable conditions 
    836              DO jk = 2, ibld(ji,jj) 
    837                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
    838                 ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj) 
    839              END DO 
    840           ENDIF 
    841        END_2D 
    842  
    843        WHERE ( lconv ) 
    844           zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    845           zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0) 
    846           zsc_vw_1 = 0._wp 
    847        ELSEWHERE 
    848          zsc_uw_1 = 0._wp 
    849          zsc_vw_1 = 0._wp 
    850        ENDWHERE 
    851  
    852        DO_2D( 0, 0, 0, 0 ) 
    853           IF ( lconv(ji,jj) ) THEN 
    854              DO jk = 2 , imld(ji,jj) 
    855                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    856                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     & 
    857                      &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   & 
    858                      &                                          * zsc_uw_2(ji,jj)                                    ) 
    859                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
    860              END DO  ! jk loop 
    861           ELSE 
    862           ! stable conditions 
    863              DO jk = 2, ibld(ji,jj) 
    864                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 
    865                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
    866              END DO 
    867           ENDIF 
    868        END_2D 
    869  
    870        DO_2D( 0, 0, 0, 0 ) 
    871         IF ( lpyc(ji,jj) ) THEN 
    872           IF ( j_ddh(ji,jj) == 0 ) THEN 
    873 ! Place holding code. Parametrization needs checking for these conditions. 
    874             zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
    875             zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
    876             zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
    877           ELSE 
    878             zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
    879             zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
    880             zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
    881           ENDIF 
    882           zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 
    883           zc_cubic = zuw_bse - zd_cubic 
    884 ! need ztau_sc_u to be available. Change to array. 
    885           DO jk = imld(ji,jj), ibld(ji,jj) 
    886              zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
    887              ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    888           END DO 
    889           zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 
    890           zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 
    891           zc_cubic = zvw_bse - zd_cubic 
    892           DO jk = imld(ji,jj), ibld(ji,jj) 
    893             zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj) 
    894             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    895           END DO 
    896         ENDIF  ! lpyc 
    897        END_2D 
    898  
    899        IF(ln_dia_osm) THEN 
    900           IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 
    901           IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
    902        END IF 
    903 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    904  
    905        DO_2D( 1, 0, 1, 0 ) 
    906  
    907          IF ( lconv(ji,jj) ) THEN 
    908            zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 
    909            zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 
    910            IF ( lpyc(ji,jj) ) THEN 
    911 ! Pycnocline scales 
    912               zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 
    913               zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 
    914             ENDIF 
    915          ELSE 
    916            zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 
    917            zsc_ws_1(ji,jj) = zws0(ji,jj) 
    918          ENDIF 
    919        END_2D 
    920  
    921        DO_2D( 0, 0, 0, 0 ) 
    922          IF ( lconv(ji,jj) ) THEN 
    923             DO jk = 2, imld(ji,jj) 
    924                zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    925                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                & 
    926                     &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
    927                     &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
    928                     &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) ) 
    929                ! 
    930                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  & 
    931                     &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
    932                     &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
    933                     &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) ) 
    934             END DO 
    935 ! 
    936             IF ( lpyc(ji,jj) ) THEN 
    937 ! pycnocline 
    938               DO jk = imld(ji,jj), ibld(ji,jj) 
    939                 zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 
    940                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 
    941                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 
    942               END DO 
    943            ENDIF 
    944          ELSE 
    945             IF( zdhdt(ji,jj) > 0. ) THEN 
    946               DO jk = 2, ibld(ji,jj) 
    947                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    948                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    949                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    950               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
    951                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    952                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
    953               END DO 
    954             ENDIF 
    955          ENDIF 
    956        END_2D 
    957  
    958        WHERE ( lconv ) 
    959           zsc_uw_1 = zustar**2 
    960           zsc_vw_1 = ff_t * zustke * zhml 
    961        ELSEWHERE 
    962           zsc_uw_1 = zustar**2 
    963           zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1 
    964           zsc_vw_1 = ff_t * zustke * zhbl 
    965           zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1 
    966        ENDWHERE 
    967  
    968        DO_2D( 0, 0, 0, 0 ) 
    969           IF ( lconv(ji,jj) ) THEN 
    970             DO jk = 2, imld(ji,jj) 
    971                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    972                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    973                ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
    974                     & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 
    975                ! 
    976                ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    977                     & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 
    978             END DO 
    979           ELSE 
    980             DO jk = 2, ibld(ji,jj) 
    981                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    982                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    983                IF ( zznd_d <= 2.0 ) THEN 
    984                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 
    985                        &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 
    986                   ! 
    987                ELSE 
    988                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
    989                        & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 
    990                   ! 
    991                ENDIF 
    992  
    993                ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    994                     & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 
    995                ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    996                     & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 
    997             END DO 
    998           ENDIF 
    999        END_2D 
    1000  
    1001        IF(ln_dia_osm) THEN 
    1002           IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 
    1003           IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 
    1004           IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
    1005           IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
    1006           IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
    1007           IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
    1008        END IF 
    1009 ! 
    1010 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
    1011  
    1012  
    1013  ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 
    1014  
    1015       DO_2D( 0, 0, 0, 0 ) 
    1016          IF ( .not. lconv(ji,jj) ) THEN 
    1017             DO jk = 2, ibld(ji,jj) 
    1018                znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 
    1019                IF ( znd >= 0.0 ) THEN 
    1020                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    1021                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    1022                ELSE 
    1023                   ghamu(ji,jj,jk) = 0._wp 
    1024                   ghamv(ji,jj,jk) = 0._wp 
    1025                ENDIF 
    1026             END DO 
    1027          ENDIF 
    1028       END_2D 
    1029  
    1030       ! Pynocline contributions 
    1031       ! 
    1032       IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays used for output of pycnocline gradient/shear profiles 
    1033          ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) 
    1034          CALL mpp_sum( 'zdfosm', istat ) 
    1035          IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 
    1036          z3ddz_pyc_1(:,:,:) = 0.0_wp 
    1037          z3ddz_pyc_2(:,:,:) = 0.0_wp 
    1038       END IF 
    1039       DO_2D( 0, 0, 0, 0 ) 
    1040          IF ( lconv (ji,jj) ) THEN 
    1041             ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
    1042             !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
    1043             !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
    1044             !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
    1045             !Alan is this right? 
    1046             !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
    1047             !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
    1048             !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
    1049             !                       &      )/ (zdh(ji,jj)  + epsln ) 
    1050             !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
    1051             !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
    1052             !                     IF ( znd <= 0.0 ) THEN 
    1053             !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
    1054             !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
    1055             !                     ELSE 
    1056             !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
    1057             !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
    1058             !                     ENDIF 
    1059             !                  END DO 
    1060          ELSE   ! Stable conditions 
    1061             IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    1062                ! If pycnocline profile only defined when depth steady of increasing. 
    1063                IF ( zdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
    1064                   IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 
    1065                      IF ( zhol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline 
    1066                         ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 
    1067                         ztgrad = zdt_bl(ji,jj) * ztmp 
    1068                         zsgrad = zds_bl(ji,jj) * ztmp 
    1069                         zbgrad = zdb_bl(ji,jj) * ztmp 
    1070                         DO jk = 2, ibld(ji,jj) 
    1071                            znd = gdepw(ji,jj,jk,Kmm) * ztmp 
    1072                            zdtdz_pyc =  ztgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 
    1073                            zdsdz_pyc =  zsgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 
    1074                            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 
    1075                            ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 
    1076                            IF ( ln_dia_pyc_scl ) THEN 
    1077                               z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 
    1078                               z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 
    1079                            END IF 
    1080                         END DO 
    1081                      ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    1082                         ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 
    1083                         ztgrad = zdt_bl(ji,jj) * ztmp 
    1084                         zsgrad = zds_bl(ji,jj) * ztmp 
    1085                         zbgrad = zdb_bl(ji,jj) * ztmp 
    1086                         DO jk = 2, ibld(ji,jj) 
    1087                            znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 
    1088                            zdtdz_pyc =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    1089                            zdsdz_pyc =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    1090                            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 
    1091                            ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 
    1092                            IF ( ln_dia_pyc_scl ) THEN 
    1093                               z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 
    1094                               z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 
    1095                            END IF 
    1096                         END DO 
    1097                      ENDIF   ! IF (zhol >=0.5) 
    1098                   ENDIF      ! IF (zdb_bl> 0.) 
    1099                ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    1100             END IF 
    1101          END IF 
    1102       END_2D 
    1103       IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles 
    1104          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 
    1105          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
    1106          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * zdbdz_pyc(:,:,:)   ) 
    1107       END IF 
    1108       DO_2D( 0, 0, 0, 0 ) 
    1109          IF ( .NOT. lconv (ji,jj) ) THEN 
    1110             IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    1111                zugrad = 3.25_wp * zdu_bl(ji,jj) / zhbl(ji,jj) 
    1112                zvgrad = 2.75_wp * zdv_bl(ji,jj) / zhbl(ji,jj) 
    1113                DO jk = 2, ibld(ji,jj) 
    1114                   znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    1115                   IF ( znd < 1.0 ) THEN 
    1116                      zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) 
    1117                   ELSE 
    1118                      zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) 
    1119                   ENDIF 
    1120                   zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) 
    1121                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc 
    1122                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc 
    1123                   IF ( ln_dia_pyc_shr ) THEN 
    1124                      z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc 
    1125                      z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc 
    1126                   END IF 
    1127                END DO 
    1128             END IF 
    1129          END IF 
    1130       END_2D 
    1131       IF ( ln_dia_pyc_shr ) THEN   ! Output of pycnocline shear profiles 
    1132          IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 
    1133          IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
    1134       END IF 
    1135       IF(ln_dia_osm) THEN 
    1136          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
    1137          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
    1138       END IF 
    1139       IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles 
    1140          DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) 
    1141       END IF 
    1142  
    1143        DO_2D( 0, 0, 0, 0 ) 
    1144           ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1145           ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1146           ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1147           ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 
    1148        END_2D 
    1149  
    1150        IF(ln_dia_osm) THEN 
    1151           IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
    1152           IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
    1153           IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
    1154        END IF 
     684      ! 
     685      ! Calculate non-gradient components of the flux-gradient relationships 
     686      ! -------------------------------------------------------------------- 
     687      CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, ibld_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear,   & 
     688         &                    zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla,     & 
     689         &                    zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml,                & 
     690         &                    zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 
    1155691 
    1156692       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    1329865         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine 
    1330866         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux 
    1331          IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux 
    1332867         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux 
    1333          IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux 
    1334868         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML 
    1335869 
     
    1374908! 
    1375909      REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 
     910      REAL(wp) ::   za_cubic, zb_cubic, zc_cubic, zd_cubic   ! Coefficients in cubic polynomial specifying diffusivity in pycnocline 
    1376911 
    1377912      REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 
     
    18831418      END_2D 
    18841419      ! 
     1420      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles 
     1421         IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) ) 
     1422      END IF 
     1423      ! 
    18851424      IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 
    18861425      ! 
     
    25422081   END SUBROUTINE zdf_osm_vertical_average 
    25432082 
     2083   SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, kbld_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear,   & 
     2084      &                          pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla,       & 
     2085      &                          pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml,                  & 
     2086      &                          pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) 
     2087      !!--------------------------------------------------------------------- 
     2088      !!                 ***  ROUTINE zdf_osm_fgr_terms *** 
     2089      !! 
     2090      !! ** Purpose : Compute non-gradient terms in flux-gradient relationship 
     2091      !! 
     2092      !! ** Method  : 
     2093      !! 
     2094      !!---------------------------------------------------------------------- 
     2095      INTEGER,                    INTENT(in   ) ::   Kmm            ! Time-level index 
     2096      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kbld           ! BL base layer 
     2097      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer 
     2098      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kp_ext         ! Offset for external level 
     2099      INTEGER,                    INTENT(in   ) ::   kbld_ext       ! Offset for external level 
     2100      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags 
     2101      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
     2102      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   k_ddh          ! Type of shear layer 
     2103      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phbl           ! BL depth 
     2104      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phml           ! ML depth 
     2105      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2106      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2107      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phol           ! Stability parameter for boundary layer 
     2108      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pshear         ! Shear production 
     2109      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pustar         ! Friction velocity 
     2110      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwstrl         ! Langmuir velocity scale 
     2111      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pvstr          ! Velocity scale (approaches zustar for large Langmuir number) 
     2112      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwstrc         ! Convective velocity scale 
     2113      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   puw0           ! Surface u-momentum flux 
     2114      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwth0          ! Surface heat flux 
     2115      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pws0           ! Surface freshwater flux 
     2116      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb0           ! Surface buoyancy flux 
     2117      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwthav         ! BL average heat flux 
     2118      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwsav          ! BL average freshwater flux 
     2119      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwbav          ! BL average buoyancy flux 
     2120      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pustke         ! Surface Stokes drift 
     2121      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pla            ! Langmuir number 
     2122      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_bl         ! Temperature diff. between BL average and basal value 
     2123      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_bl         ! Salinity diff. between BL average and basal value 
     2124      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
     2125      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdu_bl         ! Velocity diff. (u) between BL average and basal value 
     2126      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdv_bl         ! Velocity diff. (u) between BL average and basal value 
     2127      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_ml         ! Temperature diff. between mixed-layer average and basal value 
     2128      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_ml         ! Salinity diff. between mixed-layer average and basal value 
     2129      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
     2130      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdu_ml         ! Velocity diff. (u) between mixed-layer average and basal value 
     2131      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdv_ml         ! Velocity diff. (v) between mixed-layer average and basal value 
     2132      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
     2133      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
     2134      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     2135      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdbdz_pyc      ! Pycnocline buoyancy gradients 
     2136      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   palpha_pyc     ! 
     2137      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdiffut        ! t-diffusivity 
     2138      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pviscos        ! Viscosity 
     2139      ! 
     2140      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles 
     2141      ! 
     2142      INTEGER                     ::   ji, jj, jk, jkm_bld, jkf_mld, jkm_mld   ! Loop indices 
     2143      INTEGER                     ::   istat                                   ! Memory allocation status 
     2144      REAL(wp)                    ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths 
     2145      REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales 
     2146      REAL(wp), DIMENSION(A2D(0)) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales 
     2147      REAL(wp), DIMENSION(A2D(0)) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales 
     2148      REAL(wp), DIMENSION(A2D(0)) ::   ztau_sc_u                               ! Dissipation timescale at base of WML 
     2149      REAL(wp)                    ::   zbuoy_pyc_sc, zdelta_pyc                ! 
     2150      REAL(wp)                    ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale 
     2151      REAL(wp), DIMENSION(A2D(0)) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying 
     2152      REAL(wp), DIMENSION(A2D(0)) ::   zc_cubic, zd_cubic                      ! diffusivity in pycnocline 
     2153      REAL(wp), DIMENSION(A2D(0)) ::   zwt_pyc_sc_1, zws_pyc_sc_1              ! 
     2154      REAL(wp), DIMENSION(A2D(0)) ::   zzeta_pyc                               ! 
     2155      REAL(wp)                    ::   zomega, zvw_max                         ! 
     2156      REAL(wp)                    ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes 
     2157      REAL(wp), DIMENSION(A2D(0)) ::   zwth_ent,zws_ent                        ! at the top of the pycnocline 
     2158      REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term 
     2159      REAL(wp)                    ::   ztmp                                    ! 
     2160      REAL(wp)                    ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline gradients 
     2161      REAL(wp)                    ::   zugrad, zvgrad                          ! Variables for calculating pycnocline shear 
     2162      REAL(wp)                    ::   zdtdz_pyc                               ! Parametrized gradient of temperature in pycnocline 
     2163      REAL(wp)                    ::   zdsdz_pyc                               ! Parametrised gradient of salinity in pycnocline 
     2164      REAL(wp)                    ::   zdudz_pyc                               ! u-shear across the pycnocline 
     2165      REAL(wp)                    ::   zdvdz_pyc                               ! v-shear across the pycnocline 
     2166      !!---------------------------------------------------------------------- 
     2167      ! 
     2168      IF( ln_timing ) CALL timing_start('zdf_osm_ft') 
     2169      ! 
     2170      ! Auxiliary indices 
     2171      ! ----------------- 
     2172      jkm_bld = 0 
     2173      jkf_mld = jpk 
     2174      jkm_mld = 0 
     2175      DO_2D( 0, 0, 0, 0 ) 
     2176         IF ( kbld(ji,jj) > jkm_bld ) jkm_bld = kbld(ji,jj) 
     2177         IF ( kbld(ji,jj) < jkf_mld ) jkf_mld = kbld(ji,jj) 
     2178         IF ( kmld(ji,jj) > jkm_mld ) jkm_mld = kmld(ji,jj) 
     2179      END_2D 
     2180      ! 
     2181      ! Stokes term in scalar flux, flux-gradient relationship 
     2182      ! ------------------------------------------------------ 
     2183      WHERE ( ldconv(A2D(0)) ) 
     2184         zsc_wth_1(:,:) = pwstrl(A2D(0))**3 * pwth0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
     2185         zsc_ws_1(:,:)  = pwstrl(A2D(0))**3 * pws0(A2D(0))  / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
     2186      ELSEWHERE 
     2187         zsc_wth_1(:,:) = 2.0_wp * pwthav(A2D(0)) 
     2188         zsc_ws_1(:,:)  = 2.0_wp * pwsav(A2D(0)) 
     2189      ENDWHERE 
     2190      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2191         IF ( ldconv(ji,jj) ) THEN 
     2192            IF ( jk <= kmld(ji,jj) ) THEN 
     2193               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2194               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   & 
     2195                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 
     2196               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   & 
     2197                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 
     2198            END IF 
     2199         ELSE   ! Stable conditions 
     2200            IF ( jk <= kbld(ji,jj) ) THEN 
     2201               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2202               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   & 
     2203                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 
     2204               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   & 
     2205                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 
     2206            END IF 
     2207         END IF   ! Check on ldconv 
     2208      END_3D 
     2209      ! 
     2210      IF ( ln_dia_osm ) THEN 
     2211         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 
     2212         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 
     2213      END IF 
     2214      ! 
     2215      ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use 
     2216      ! zvstr since term needs to go to zero as zwstrl goes to zero) 
     2217      ! --------------------------------------------------------------------- 
     2218      WHERE ( ldconv(A2D(0)) ) 
     2219         zsc_uw_1(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) /   & 
     2220            &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * pla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 
     2221         zsc_uw_2(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) /   & 
     2222            &                                  MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 
     2223         zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
     2224            &            ( ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 
     2225      ELSEWHERE 
     2226         zsc_uw_1(:,:) = pustar(A2D(0))**2 
     2227         zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
     2228            &            ( pvstr(A2D(0))**2 + epsln ) 
     2229      ENDWHERE 
     2230      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2231         IF ( ldconv(ji,jj) ) THEN 
     2232            IF ( jk <= kmld(ji,jj) ) THEN 
     2233               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2234               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp   * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) +     & 
     2235                  &                                  0.00125_wp * EXP( -1.0_wp * zznd_d ) * zsc_uw_2(ji,jj) ) *   & 
     2236                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) 
     2237               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65_wp *  0.15_wp * EXP( -1.0_wp * zznd_d ) *                 & 
     2238                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_vw_1(ji,jj) 
     2239            END IF 
     2240         ELSE   ! Stable conditions 
     2241            IF ( jk <= kbld(ji,jj) ) THEN   ! Corrected to ibld 
     2242               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2243               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) *             & 
     2244                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_uw_1(ji,jj) 
     2245            END IF 
     2246         END IF 
     2247      END_3D 
     2248      ! 
     2249      ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio 
     2250      ! (X0.3) and pressure (X0.5)] 
     2251      ! ---------------------------------------------------------------------- 
     2252      WHERE ( ldconv(A2D(0)) ) 
     2253         zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   & 
     2254            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
     2255         zsc_ws_1(:,:)  = pwbav(A2D(0)) * pws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   & 
     2256            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
     2257      ELSEWHERE 
     2258         zsc_wth_1(:,:) = 0.0_wp 
     2259         zsc_ws_1(:,:)  = 0.0_wp 
     2260      ENDWHERE 
     2261      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2262         IF ( ldconv(ji,jj) ) THEN 
     2263            IF ( jk <= kmld(ji,jj) ) THEN 
     2264               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
     2265               ! Calculate turbulent length scale 
     2266               zl_c   = 0.9_wp * ( 1.0_wp - EXP( -7.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         & 
     2267                  &     ( 1.0_wp - EXP( -15.0_wp * ( 1.1_wp - zznd_ml ) ) ) 
     2268               zl_l   = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         & 
     2269                  &     ( 1.0_wp - EXP( -5.0_wp  * ( 1.0_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 
     2270               zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) 
     2271               ! Non-gradient buoyancy terms 
     2272               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.5_wp * zsc_wth_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 
     2273               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.5_wp *  zsc_ws_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 
     2274            END IF 
     2275         ELSE   ! Stable conditions 
     2276            IF ( jk <= kbld(ji,jj) ) THEN 
     2277               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
     2278               ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj) 
     2279            END IF 
     2280         END IF 
     2281      END_3D 
     2282      DO_2D( 0, 0, 0, 0 ) 
     2283         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 
     2284            ztau_sc_u(ji,jj)    = phml(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *                             & 
     2285               &                ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**1.5_wp ) 
     2286            zwth_ent(ji,jj)     = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *   & 
     2287               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj) 
     2288            zws_ent(ji,jj)      = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *   & 
     2289               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 
     2290            zbuoy_pyc_sc        = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 
     2291            zdelta_pyc          = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird /   & 
     2292               &                  SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 
     2293            zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   & 
     2294               &                  zdelta_pyc**2 / pdh(ji,jj) 
     2295            zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   & 
     2296               &                  zdelta_pyc**2 / pdh(ji,jj) 
     2297            zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 
     2298         END IF 
     2299      END_2D 
     2300      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2301         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2302            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
     2303            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) -                                                                                & 
     2304               &              0.045_wp * ( ( zwth_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
     2305               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 
     2306            ghams(ji,jj,jk) = ghams(ji,jj,jk) -                                                                                & 
     2307               &              0.045_wp * ( ( zws_ent(ji,jj)  * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 & 
     2308               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 
     2309            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              & 
     2310               &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
     2311               &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
     2312            ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp  * zws_pyc_sc_1(ji,jj) *                              & 
     2313               &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        & 
     2314               &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 
     2315         END IF   ! End of pycnocline 
     2316      END_3D 
     2317      ! 
     2318      IF(ln_dia_osm) THEN 
     2319         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! Upward turb. temperature entrainment flux 
     2320         IF ( iom_use("zws_ent")  ) CALL iom_put( "zws_ent",  tmask(:,:,1)*zws_ent  )   ! Upward turb. salinity    entrainment flux 
     2321      END IF 
     2322      ! 
     2323      zsc_vw_1(:,:) = 0.0_wp 
     2324      WHERE ( ldconv(A2D(0)) ) 
     2325         zsc_uw_1(:,:) = -1.0_wp * pwb0(A2D(0)) * pustar(A2D(0))**2 * phml(A2D(0)) /   & 
     2326            &            ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 
     2327         zsc_uw_2(:,:) =           pwb0(A2D(0)) * pustke(A2D(0))    * phml(A2D(0)) /   & 
     2328            &            ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp ) 
     2329      ELSEWHERE 
     2330         zsc_uw_1(:,:) = 0.0_wp 
     2331      ENDWHERE 
     2332      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2333         IF ( ldconv(ji,jj) ) THEN 
     2334            IF ( jk <= kmld(ji,jj) ) THEN 
     2335               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2336               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp *   & 
     2337                  &                                ( zsc_uw_1(ji,jj) + 0.125_wp * EXP( -0.5_wp * zznd_d ) *       & 
     2338                  &                                  (   1.0_wp - EXP( -0.5_wp * zznd_d ) ) * zsc_uw_2(ji,jj) ) 
     2339               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
     2340            END IF 
     2341         ELSE   ! Stable conditions 
     2342            IF ( jk <= kbld(ji,jj) ) THEN 
     2343               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 
     2344               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
     2345            END IF 
     2346         ENDIF 
     2347      END_3D 
     2348      ! 
     2349      DO_2D( 0, 0, 0, 0 ) 
     2350         IF ( ldpyc(ji,jj) ) THEN 
     2351            IF ( k_ddh(ji,jj) == 0 ) THEN 
     2352               ! Place holding code. Parametrization needs checking for these conditions. 
     2353               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird 
     2354               zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
     2355               zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
     2356            ELSE 
     2357               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird 
     2358               zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 
     2359               zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 
     2360            ENDIF 
     2361            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse 
     2362            za_cubic(ji,jj) = zuw_bse - zb_cubic(ji,jj) 
     2363            zvw_max = 0.7_wp * ff_t(ji,jj) * ( pustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * pustar(ji,jj) * phml(ji,jj) ) 
     2364            zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse 
     2365            zc_cubic(ji,jj) = zvw_bse - zd_cubic(ji,jj) 
     2366         END IF 
     2367      END_2D 
     2368      DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array. 
     2369         IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2370            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
     2371            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 *                                      & 
     2372               &                                ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) *   & 
     2373               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
     2374            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 *                                      & 
     2375               &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   & 
     2376               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 
     2377         END IF   ! ldpyc 
     2378      END_3D 
     2379      ! 
     2380      IF(ln_dia_osm) THEN 
     2381         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0",    wmask*ghamu           ) 
     2382         IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 
     2383      END IF 
     2384      ! 
     2385      ! Transport term in flux-gradient relationship [note : includes ROI ratio 
     2386      ! (X0.3) ] 
     2387      ! ----------------------------------------------------------------------- 
     2388      WHERE ( ldconv(A2D(0)) ) 
     2389         zsc_wth_1(:,:) = pwth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) 
     2390         zsc_ws_1(:,:)  = pws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) 
     2391         WHERE ( ldpyc(A2D(0)) )   ! Pycnocline scales 
     2392            zsc_wth_pyc(:,:) = -0.2_wp * pwb0(A2D(0)) * pdt_bl(A2D(0)) / pdb_bl(A2D(0)) 
     2393            zsc_ws_pyc(:,:)  = -0.2_wp * pwb0(A2D(0)) * pds_bl(A2D(0)) / pdb_bl(A2D(0)) 
     2394         END WHERE 
     2395      ELSEWHERE 
     2396         zsc_wth_1(:,:) = 2.0 * pwthav(A2D(0)) 
     2397         zsc_ws_1(:,:)  =       pws0(A2D(0)) 
     2398      END WHERE 
     2399      DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) 
     2400         IF ( ldconv(ji,jj) ) THEN 
     2401            IF ( ( jk > 1 ) .AND. ( jk <= kmld(ji,jj) ) ) THEN 
     2402               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
     2403               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   & 
     2404                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) 
     2405               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj)  * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   & 
     2406                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) 
     2407            END IF 
     2408            IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN   ! Pycnocline 
     2409               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
     2410               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) 
     2411               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj)  * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) 
     2412            END IF 
     2413         ELSE 
     2414            IF( pdhdt(ji,jj) > 0. ) THEN 
     2415               IF ( ( jk > 1 ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2416                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2417                  znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
     2418                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) +   & 
     2419                                      7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_wth_1(ji,jj) 
     2420                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) +   & 
     2421                                      7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_ws_1(ji,jj) 
     2422               END IF 
     2423            ENDIF 
     2424         ENDIF 
     2425      END_3D 
     2426      ! 
     2427      WHERE ( ldconv(A2D(0)) ) 
     2428         zsc_uw_1(:,:) = pustar(A2D(0))**2 
     2429         zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) 
     2430      ELSEWHERE 
     2431         zsc_uw_1(:,:) = pustar(A2D(0))**2 
     2432         zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) *   & 
     2433            &            zsc_uw_1(:,:) 
     2434         zsc_vw_1(:,:) = ff_t * pustke(A2D(0)) * phbl(A2D(0)) 
     2435         zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) *   & 
     2436            &            zsc_vw_1(:,:) 
     2437      ENDWHERE 
     2438      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2439         IF ( ldconv(ji,jj) ) THEN 
     2440            IF ( jk <= kmld(ji,jj) ) THEN 
     2441               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 
     2442               zznd_d  = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2443               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +   & 
     2444                  &              0.3_wp * ( -2.0_wp + 2.5_wp * ( 1.0_wp + 0.1_wp * zznd_ml**4 ) - EXP( -8.0_wp * zznd_ml ) ) *   & 
     2445                  &              zsc_uw_1(ji,jj) 
     2446               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) +   & 
     2447                  &              0.3_wp * 0.1_wp * ( EXP( -1.0_wp * zznd_d ) + EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) *   & 
     2448                  &              zsc_vw_1(ji,jj) 
     2449            END IF 
     2450         ELSE 
     2451            IF ( jk <= kbld(ji,jj) ) THEN 
     2452               znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
     2453               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     2454               IF ( zznd_d <= 2.0 ) THEN 
     2455                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *                                              & 
     2456                     &                                ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) *   & 
     2457                     &                                  ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ) * zsc_uw_1(ji,jj) 
     2458               ELSE 
     2459                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *   & 
     2460                     &                                ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - znd ) ) ) * zsc_uw_2(ji,jj) 
     2461               ENDIF 
     2462               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) *   & 
     2463                  &                                EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj) 
     2464               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 
     2465            END IF 
     2466         END IF 
     2467      END_3D 
     2468      ! 
     2469      IF(ln_dia_osm) THEN 
     2470         IF ( iom_use("ghamu_f") )    CALL iom_put( "ghamu_f",    wmask       *ghamu    ) 
     2471         IF ( iom_use("ghamv_f") )    CALL iom_put( "ghamv_f",    wmask       *ghamv    ) 
     2472         IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 
     2473         IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 
     2474         IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 
     2475         IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 
     2476      END IF 
     2477      ! 
     2478      ! Make surface forced velocity non-gradient terms go to zero at the base 
     2479      ! of the mixed layer. 
     2480      ! 
     2481      ! Make surface forced velocity non-gradient terms go to zero at the base 
     2482      ! of the boundary layer. 
     2483      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2484         IF ( ( .NOT. ldconv(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 
     2485            znd = ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about 
     2486            IF ( znd >= 0.0_wp ) THEN 
     2487               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) 
     2488               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) 
     2489            ELSE 
     2490               ghamu(ji,jj,jk) = 0.0_wp 
     2491               ghamv(ji,jj,jk) = 0.0_wp 
     2492            ENDIF 
     2493         END IF 
     2494      END_3D 
     2495      ! 
     2496      ! Pynocline contributions 
     2497      ! 
     2498      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays for output of pycnocline gradient/shear profiles 
     2499         ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) 
     2500         IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 
     2501         z3ddz_pyc_1(:,:,:) = 0.0_wp 
     2502         z3ddz_pyc_2(:,:,:) = 0.0_wp 
     2503      END IF 
     2504      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2505         IF ( ldconv (ji,jj) ) THEN 
     2506            ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
     2507            !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
     2508            !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
     2509            !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
     2510            !Alan is this right? 
     2511            !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
     2512            !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
     2513            !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
     2514            !                       &      )/ (zdh(ji,jj)  + epsln ) 
     2515            !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
     2516            !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
     2517            !                     IF ( znd <= 0.0 ) THEN 
     2518            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
     2519            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
     2520            !                     ELSE 
     2521            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
     2522            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
     2523            !                     ENDIF 
     2524            !                  END DO 
     2525         ELSE   ! Stable conditions 
     2526            IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     2527               ! Pycnocline profile only defined when depth steady of increasing. 
     2528               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
     2529                  IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 
     2530                     IF ( phol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline 
     2531                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 
     2532                        ztgrad = pdt_bl(ji,jj) * ztmp 
     2533                        zsgrad = pds_bl(ji,jj) * ztmp 
     2534                        zbgrad = pdb_bl(ji,jj) * ztmp 
     2535                        IF ( jk <= kbld(ji,jj) ) THEN 
     2536                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
     2537                           zdtdz_pyc =  ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
     2538                           zdsdz_pyc =  zsgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
     2539                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc 
     2540                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc 
     2541                           IF ( ln_dia_pyc_scl ) THEN 
     2542                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 
     2543                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 
     2544                           END IF 
     2545                        END IF 
     2546                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     2547                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 
     2548                        ztgrad = pdt_bl(ji,jj) * ztmp 
     2549                        zsgrad = pds_bl(ji,jj) * ztmp 
     2550                        zbgrad = pdb_bl(ji,jj) * ztmp 
     2551                        IF ( jk <= kbld(ji,jj) ) THEN 
     2552                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 
     2553                           zdtdz_pyc =  ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
     2554                           zdsdz_pyc =  zsgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
     2555                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc 
     2556                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc 
     2557                           IF ( ln_dia_pyc_scl ) THEN 
     2558                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 
     2559                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 
     2560                           END IF 
     2561                        END IF 
     2562                     ENDIF   ! IF (zhol >=0.5) 
     2563                  ENDIF      ! IF (zdb_bl> 0.) 
     2564               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     2565            END IF 
     2566         END IF 
     2567      END_3D 
     2568      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles 
     2569         IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 
     2570         IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
     2571      END IF 
     2572      DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2573         IF ( .NOT. ldconv (ji,jj) ) THEN 
     2574            IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     2575               zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj) 
     2576               zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj) 
     2577               IF ( jk <= kbld(ji,jj) ) THEN 
     2578                  znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 
     2579                  IF ( znd < 1.0 ) THEN 
     2580                     zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) 
     2581                  ELSE 
     2582                     zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) 
     2583                  ENDIF 
     2584                  zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) 
     2585                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + pviscos(ji,jj,jk) * zdudz_pyc 
     2586                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + pviscos(ji,jj,jk) * zdvdz_pyc 
     2587                  IF ( ln_dia_pyc_shr ) THEN 
     2588                     z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc 
     2589                     z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc 
     2590                  END IF 
     2591               END IF 
     2592            END IF 
     2593         END IF 
     2594      END_3D 
     2595      IF ( ln_dia_pyc_shr ) THEN   ! Output of pycnocline shear profiles 
     2596         IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 
     2597         IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
     2598      END IF 
     2599      IF(ln_dia_osm) THEN 
     2600         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
     2601         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
     2602      END IF 
     2603      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles 
     2604         DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) 
     2605      END IF 
     2606      ! 
     2607      DO_2D( 0, 0, 0, 0 ) 
     2608         ghamt(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
     2609         ghams(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
     2610         ghamu(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
     2611         ghamv(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 
     2612      END_2D 
     2613      ! 
     2614      IF(ln_dia_osm) THEN 
     2615         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu   ) 
     2616         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv   ) 
     2617         IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*pviscos ) 
     2618      END IF 
     2619      ! 
     2620      IF( ln_timing ) CALL timing_stop('zdf_osm_ft') 
     2621      ! 
     2622   END SUBROUTINE zdf_osm_fgr_terms 
     2623 
    25442624   SUBROUTINE zdf_osm_init( Kmm ) 
    25452625     !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.