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

Changeset 14679


Ignore:
Timestamp:
2021-04-07T17:52:57+02:00 (3 years ago)
Author:
dancopsey
Message:

Merge in revisions 14541 to 14678 of:

​http://forge.ipsl.jussieu.fr/nemo/log/NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0

Location:
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/field_def_nemo-oce.xml

    r14413 r14679  
    240240        <field id="zds_ml"             long_name="salinity jump at base of ML"                 unit="10^-3"      /> 
    241241        <field id="zdb_ml"             long_name="buoyancy jump at base of ML"                 unit="m/s^2"      />  
     242        <field id="pb_coup"             long_name="bottom coupling velocity"                           unit="m/s"      />  
    242243        <!-- extra OSMOSIS diagnostics for debugging --> 
    243244       <field id="zsc_uw_1_0"       long_name="zsc u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 

  • NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90

    r14543 r14679  
    5757  USE eosbn2         ! equation of state 
    5858  USE traqsr         ! details of solar radiation absorption 
     59  USE zdfdrg,  ONLY  : rCdU_bot  ! bottom friction velocity 
    5960  USE zdfddm         ! double diffusion mixing (avs array) 
    6061  !   USE zdfmxl         ! mixed layer depth 
    … …  
    273274    LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl 
    274275    LOGICAL, DIMENSION(jpi,jpj)  :: lshear    ! Shear layers 
     276    LOGICAL, DIMENSION(jpi,jpj)  :: lcoup     ! Coupling to bottom 
    275277    LOGICAL, DIMENSION(jpi,jpj)  :: lpyc      ! OSBL pycnocline present 
    276278    LOGICAL, DIMENSION(jpi,jpj)  :: lflux     ! surface flux extends below OSBL into MLE layer. 
    … …  
    354356    ! For debugging 
    355357    INTEGER :: ikt 
     358    REAL(wp) :: zlarge = -1.e10_wp, zero = 0._wp 
    356359    !!-------------------------------------------------------------------- 
    357360    ! 
    358361    ibld(:,:)   = 0     ; imld(:,:)  = 0 
    359     zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
    360     zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp 
    361     zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp 
    362     zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp 
    363     zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 
    364     zhol(:,:)   = 0._wp ; zwb0tot(:,:) = 0._wp 
     362    zrad0(:,:)  = zlarge ; zradh(:,:) = zlarge ; zradav(:,:)    = zlarge ; zustar(:,:)    = zlarge 
     363    zwstrl(:,:) = zlarge ; zvstr(:,:) = zlarge ; zwstrc(:,:)    = zlarge ; zuw0(:,:)      = zlarge 
     364    zvw0(:,:)   = zlarge ; zwth0(:,:) = zlarge ; zws0(:,:)      = zlarge ; zwb0(:,:)      = zlarge 
     365    zwthav(:,:) = zlarge ; zwsav(:,:) = zlarge ; zwbav(:,:)     = zlarge ; zwb_ent(:,:)   = zlarge 
     366    zustke(:,:) = zlarge ; zla(:,:)   = zlarge ; zcos_wind(:,:) = zlarge ; zsin_wind(:,:) = zlarge 
     367    zhol(:,:)   = zlarge ; zwb0tot(:,:) = zlarge; zalpha_pyc(:,:) = zlarge 
    365368    lconv(:,:)  = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE. 
    366369    ! mixed layer 
    367370    ! no initialization of zhbl or zhml (or zdh?) 
    368     zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp 
    369     zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp 
    370     zv_bl(:,:)   = 0._wp ; zb_bl(:,:)  = 0._wp 
    371     zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp 
    372     zt_mle(:,:)   = 0._wp ; zs_mle(:,:)    = 0._wp ; zu_mle(:,:)   = 0._wp 
    373     zb_mle(:,:) = 0._wp 
    374     zv_ml(:,:)   = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp 
    375     zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp 
    376     zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp 
    377     zdb_ml(:,:)  = 0._wp 
    378     zdt_mle(:,:)  = 0._wp ; zds_mle(:,:)  = 0._wp ; zdu_mle(:,:)   = 0._wp 
    379     zdv_mle(:,:)  = 0._wp ; zdb_mle(:,:)  = 0._wp 
    380     zwth_ent = 0._wp ; zws_ent = 0._wp 
     371    zhbl(:,:)    = zlarge ; zhml(:,:)    = zlarge ; zdh(:,:)      = zlarge ; zdhdt(:,:)   = zlarge 
     372    zt_bl(:,:)   = zlarge ; zs_bl(:,:)   = zlarge ; zu_bl(:,:)    = zlarge 
     373    zv_bl(:,:)   = zlarge ; zb_bl(:,:)  = zlarge 
     374    zt_ml(:,:)   = zlarge ; zs_ml(:,:)    = zlarge ; zu_ml(:,:)   = zlarge 
     375    zt_mle(:,:)   = zlarge ; zs_mle(:,:)    = zlarge ; zu_mle(:,:)   = zlarge 
     376    zb_mle(:,:) = zlarge 
     377    zv_ml(:,:)   = zlarge ; zdt_bl(:,:)   = zlarge ; zds_bl(:,:)  = zlarge 
     378    zdu_bl(:,:)  = zlarge ; zdv_bl(:,:)  = zlarge ; zdb_bl(:,:)  = zlarge 
     379    zdt_ml(:,:)  = zlarge ; zds_ml(:,:)  = zlarge ; zdu_ml(:,:)   = zlarge ; zdv_ml(:,:)  = zlarge 
     380    zdb_ml(:,:)  = zlarge 
     381    zdt_mle(:,:)  = zlarge ; zds_mle(:,:)  = zlarge ; zdu_mle(:,:)   = zlarge 
     382    zdv_mle(:,:)  = zlarge ; zdb_mle(:,:)  = zlarge 
     383    zwth_ent = zlarge ; zws_ent = zlarge 
    381384    ! 
    382     zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 
    383     zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
     385    zdtdz_pyc(:,:,:) = zlarge ; zdsdz_pyc(:,:,:) = zlarge ; zdbdz_pyc(:,:,:) = zlarge 
     386    zdudz_pyc(:,:,:) = zlarge ; zdvdz_pyc(:,:,:) = zlarge 
     387    zdtdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdsdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdbdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp 
     388    zdudz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdvdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp 
    384389    ! 
    385     zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 
     390    zdtdz_bl_ext(:,:) = zlarge ; zdsdz_bl_ext(:,:) = zlarge ; zdbdz_bl_ext(:,:) = zlarge 
    386391 
    387392    IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
    388        zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 
    389        zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 
    390        zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 
    391        zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp 
     393       zdtdx(:,:) = zlarge ; zdtdy(:,:) = zlarge ; zdsdx(:,:) = zlarge 
     394       zdsdy(:,:) = zlarge ; dbdx_mle(:,:) = zlarge ; dbdy_mle(:,:) = zlarge 
     395       zwb_fk(:,:) = zlarge ; zvel_mle(:,:) = zlarge; zdiff_mle(:,:) = zlarge 
     396       zhmle(:,:) = zlarge  ; zmld(:,:) = zlarge 
    392397    ENDIF 
    393     zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
     398    zwb_fk_b(:,:) = zlarge   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
    394399 
    395400    ! Flux-Gradient arrays. 
    396     zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp 
    397     zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp 
    398     zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp 
    399  
    400     zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp 
    401     ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp 
     401    zsc_wth_1(:,:)  = zlarge ; zsc_ws_1(:,:)   = zlarge ; zsc_uw_1(:,:)   = zlarge 
     402    zsc_uw_2(:,:)   = zlarge ; zsc_vw_1(:,:)   = zlarge ; zsc_vw_2(:,:)   = zlarge 
     403    zhbl_t(:,:)     = zlarge ; zdhdt(:,:)      = zlarge 
     404 
     405    zdiffut(:,:,:) = zlarge ; zviscos(:,:,:) = zlarge 
     406    zdiffut(2:jpim1,2:jpjm1,:) = 0._wp ; zviscos(2:jpim1,2:jpjm1,:) = 0._wp 
     407    ghamt(:,:,:) = zlarge; ghams(:,:,:) = zlarge 
     408    ghamt(2:jpim1,2:jpjm1,:) = 0._wp; ghams(2:jpim1,2:jpjm1,:) = 0._wp 
     409    ghamu(:,:,:)   = zlarge ; ghamv(:,:,:) = zlarge 
     410    ghamu(2:jpim1,2:jpjm1,:)   = 0._wp ; ghamv(2:jpim1,2:jpjm1,:) = 0._wp 
     411    zdiff_mle(2:jpim1,2:jpjm1) = 0._wp 
    402412 
    403413 
    … …  
    614624             zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 
    615625          ELSE 
     626             zwstrc(ji,jj) = 0.0_wp 
    616627             zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln ) 
    617628          ENDIF 
    618629#ifdef key_osm_debug 
    619630          IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    620              WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,2(a,g11.3),/)') & 
     631             WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,3(a,g11.3),/)') & 
    621632                  & 'After reduction: zustke=', zustke(ji,jj), ' dstokes=', dstokes(ji,jj), & 
    622633                  & ' zustar =', zustar(ji,jj), ' zwstrl=', zwstrl(ji,jj), ' zwstrc=', zwstrc(ji,jj),& 
    623                   & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj) 
     634                  & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj), ' zvstr=', zvstr(ji,jj) 
    624635             FLUSH(narea+100) 
    625636          END IF 
    … …  
    644655          DO ji = 1, jpi 
    645656             IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    646                 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
     657                ibld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 
    647658             ENDIF 
    648659          END DO 
    … …  
    682693 
    683694    ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 
    684     jp_ext(:,:) = 2 
     695    jp_ext(:,:) = 1   ! ag 19/03 
    685696    CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 
    686     !      jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 
    687     CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     697    jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 
     698    CALL zdf_osm_vertical_average(imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
    688699#ifdef key_osm_debug 
    689700    IF(narea==nn_narea_db) THEN 
    … …  
    793804 
    794805    !! External gradient below BL needed both with and w/o FK 
    795     CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     806    CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )    ! ag 19/03 
    796807 
    797808    ! Test if pycnocline well resolved 
    798     DO jj = 2, jpjm1 
    799        DO ji = 2,jpim1 
    800           IF (lconv(ji,jj) ) THEN 
    801              ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 
    802              IF ( ztmp > 6 ) THEN 
    803                 ! pycnocline well resolved 
    804                 jp_ext(ji,jj) = 1 
    805              ELSE 
    806                 ! pycnocline poorly resolved 
    807                 jp_ext(ji,jj) = 0 
    808              ENDIF 
    809           ELSE 
    810              ! Stable conditions 
    811              jp_ext(ji,jj) = 0 
    812           ENDIF 
    813        END DO 
    814     END DO 
     809!    DO jj = 2, jpjm1                                             Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 
     810!       DO ji = 2,jpim1                                           should account for this. 
     811!          IF (lconv(ji,jj) ) THEN 
     812!             ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 
     813!             IF ( ztmp > 3 ) THEN           ! ag 19/03 
     814!                ! pycnocline well resolved 
     815!               jp_ext(ji,jj) = 1 
     816!             ELSE 
     817!                ! pycnocline poorly resolved 
     818!                jp_ext(ji,jj) = 0 
     819!             ENDIF 
     820!          ELSE 
     821!             ! Stable conditions 
     822!             jp_ext(ji,jj) = 0 
     823!          ENDIF 
     824!       END DO 
     825!    END DO 
    815826#ifdef key_osm_debug 
    816827    IF(narea==nn_narea_db) THEN 
    … …  
    825836 
    826837    ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 
     838    jp_ext(:,:) = 1    ! ag 19/03 
    827839    CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    828     !      jp_ext = ibld-imld+1 
    829     CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 
     840    jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 
     841    CALL zdf_osm_vertical_average(imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)    ! ag 19/03 
    830842#ifdef key_osm_debug 
    831843    IF(narea==nn_narea_db) THEN 
    … …  
    850862          zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 
    851863          ! adjustment to represent limiting by ocean bottom 
    852           IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN 
    853              zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 
    854              lpyc(ji,jj) = .FALSE. 
     864          IF ( mbkt(ji,jj) >2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access 
     865             IF( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) - 2 ) ) THEN 
     866                zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) - 2) - depth_tol)! ht_n(:,:)) 
     867                lpyc(ji,jj) = .FALSE. 
     868             ENDIF 
    855869          ENDIF 
    856870#ifdef key_osm_debug 
    … …  
    877891    END DO 
    878892 
     893! Test if surface boundary layer coupled to bottom. 
     894 
     895    lcoup(:,:) = .FALSE.                                    ! ag 19/03 
     896    DO jj = 2, jpjm1                                        ! ag 19/03 
     897      DO ji = 2, jpim1                                      ! ag 19/03 
     898        IF ( mbkt(ji,jj) >2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access 
     899           IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,mbkt(ji,jj) - 2) ) THEN ! ag 19/03 
     900              zhbl_t(ji,jj) = gdepw_n(ji,jj,mbkt(ji,jj)-2)     ! ag 19/03 
     901              lcoup(ji,jj) = .TRUE.                            ! ag 19/03 
     902           ENDIF                                               ! ag 19/03 
     903        ENDIF                                               ! ag 19/03 
     904      END DO                                                ! ag 19/03 
     905    END DO                                                  ! ag 19/03 
     906 
    879907    ! 
    880908    ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    … …  
    884912 
    885913    !   Recalculate BL averages and differences using new BL depth 
     914    jp_ext(:,:) = 1                  ! ag 19/03 
    886915    CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 
    887916    ! 
    888     ! 
    889     ! Check to see if lpyc needs to be changed  
    890  
     917     
    891918    CALL zdf_osm_pycnocline_thickness( dh, zdh ) 
     919 
     920! reset l_pyc before calculating terms in the flux-gradient relationships 
    892921 
    893922    DO jj = 2, jpjm1 
    894923       DO ji = 2, jpim1 
    895           IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE.  
     924          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) >= mbkt(ji,jj) -2 .or. ibld(ji,jj)-imld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0._wp) THEN                                            ! ag 19/03 
     925             lpyc(ji,jj) = .FALSE.                         ! ag 19/03 
     926             IF ( ibld(ji,jj) >= mbkt(ji,jj) -2 ) THEN 
     927                imld(ji,jj) = ibld(ji,jj) - 1                 ! ag 19/03 
     928                zdh(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) - gdepw_n(ji,jj,imld(ji,jj))   ! ag 19/03 
     929                zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))      ! ag 19/03 
     930                dh(ji,jj) = zdh(ji,jj)                        ! ag 19/03   
     931                hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)           ! ag 19/03 
     932#ifdef key_osm_debug 
     933                IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     934                   WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= Fl: ibld(ji,jj) >= mbkt(ji,jj) -2' 
     935                   WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj) 
     936                   WRITE(narea+100,'(2(a,g11.3))')'dh=',dh,' hml=',hml(ji,jj) 
     937                   FLUSH(narea+100) 
     938                END IF 
     939#endif 
     940             ENDIF 
     941          ENDIF                                            ! ag 19/03 
    896942       END DO 
    897943    END DO 
    … …  
    902948    !      jp_ext = ibld - imld +1 
    903949    !     Recalculate ML averages and differences using new ML depth 
    904     CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    905     ! rotate mean currents and changes onto wind align co-ordinates 
    906     ! 
     950    jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1         ! ag 19/03 
     951    CALL zdf_osm_vertical_average( imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
     952 
     953    CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    907954#ifdef key_osm_debug 
    908955    IF(narea==nn_narea_db) THEN 
    … …  
    919966    END IF 
    920967#endif 
     968 
     969! rotate mean currents and changes onto wind align co-ordinates 
     970  
    921971    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 
    922972    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 
    … …  
    935985    !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    936986 
    937     CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     987    jp_ext(:,:) = 1                  ! ag 19/03 
    938988    CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 
    939989    CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
    … …  
    9741024    ! 
    9751025    ! Stokes term in scalar flux, flux-gradient relationship 
    976     WHERE ( lconv ) 
    977        zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln) 
     1026    WHERE ( lconv(2:jpim1,2:jpjm1) ) 
     1027       zsc_wth_1(2:jpim1,2:jpjm1) = zwstrl(2:jpim1,2:jpjm1)**3 * zwth0(2:jpim1,2:jpjm1) / & 
     1028            &( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln) 
    9781029       ! 
    979        zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
     1030       zsc_ws_1(2:jpim1,2:jpjm1) = zwstrl(2:jpim1,2:jpjm1)**3 * zws0(2:jpim1,2:jpjm1) / & 
     1031            &( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln ) 
    9801032    ELSEWHERE 
    981        zsc_wth_1 = 2.0 * zwthav 
     1033       zsc_wth_1(2:jpim1,2:jpjm1) = 2.0 * zwthav(2:jpim1,2:jpjm1) 
    9821034       ! 
    983        zsc_ws_1 = 2.0 * zwsav 
     1035       zsc_ws_1(2:jpim1,2:jpjm1) = 2.0 * zwsav(2:jpim1,2:jpjm1) 
    9841036    ENDWHERE 
    9851037 
    … …  
    10091061 
    10101062    ! 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) 
    1011     WHERE ( lconv ) 
    1012        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 ) 
    1013        zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 
    1014        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 ) 
     1063    WHERE ( lconv(2:jpim1,2:jpjm1) ) 
     1064       zsc_uw_1(2:jpim1,2:jpjm1) = ( zwstrl(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 )**pthird * zustke(2:jpim1,2:jpjm1) / & 
     1065            & MAX( ( 1.0 - 1.0 * 6.5 * zla(2:jpim1,2:jpjm1)**(8.0/3.0) ), 0.2 ) 
     1066       zsc_uw_2(2:jpim1,2:jpjm1) = ( zwstrl(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 )**pthird * zustke(2:jpim1,2:jpjm1) / & 
     1067            & MIN( zla(2:jpim1,2:jpjm1)**(8.0/3.0) + epsln, 0.12 ) 
     1068       zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zhml(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1)**3 * MIN( zla(2:jpim1,2:jpjm1)**(8.0/3.0), 0.12 ) / & 
     1069            & ( ( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 )**(2.0/3.0) + epsln ) 
    10151070    ELSEWHERE 
    1016        zsc_uw_1 = zustar**2 
    1017        zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 
     1071       zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2 
     1072       zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zhbl(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1)**3 * & 
     1073            &  MIN( zla(2:jpim1,2:jpjm1)**(8.0/3.0), 0.12 ) / (zvstr(2:jpim1,2:jpjm1)**2 + epsln) 
    10181074    ENDWHERE 
    10191075    IF(ln_dia_osm) THEN 
    … …  
    10661122    ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] 
    10671123 
    1068     WHERE ( lconv ) 
    1069        zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    1070        zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
     1124    WHERE ( lconv(2:jpim1,2:jpjm1) ) 
     1125       zsc_wth_1(2:jpim1,2:jpjm1) = zwbav(2:jpim1,2:jpjm1) * zwth0(2:jpim1,2:jpjm1) * ( 1.0 + EXP ( 0.2 * zhol(2:jpim1,2:jpjm1) ) ) * & 
     1126            & zhml(2:jpim1,2:jpjm1) / ( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln ) 
     1127       zsc_ws_1(2:jpim1,2:jpjm1)  = zwbav(2:jpim1,2:jpjm1) * zws0(2:jpim1,2:jpjm1)  * ( 1.0 + EXP ( 0.2 * zhol(2:jpim1,2:jpjm1) ) ) * & 
     1128            & zhml(2:jpim1,2:jpjm1) / ( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln ) 
    10711129    ELSEWHERE 
    1072        zsc_wth_1 = 0._wp 
    1073        zsc_ws_1 = 0._wp 
     1130       zsc_wth_1(2:jpim1,2:jpjm1) = 0._wp 
     1131       zsc_ws_1(2:jpim1,2:jpjm1) = 0._wp 
    10741132    ENDWHERE 
    10751133 
    … …  
    10891147                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 *  zsc_ws_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) 
    10901148             END DO 
    1091  
    10921149             IF ( lpyc(ji,jj) ) THEN 
    10931150                ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 
    … …  
    11021159                   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 ) 
    11031160                END DO 
    1104                 ! 
    1105                 IF ( dh(ji,jj)  <  0.2*hbl(ji,jj) ) THEN 
     1161#ifdef key_osm_debug 
     1162                jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
     1163                IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     1164                   WRITE(narea+100,'(3(a,g11.3))')'lpyc= lconv=T: ztau_sc_u=',ztau_sc_u(ji,jj),' zwth_ent=',zwth_ent,' zws_ent=',zws_ent 
     1165                   WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
     1166                   WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) 
     1167                   WRITE(narea+100,*) 
     1168                   FLUSH(narea+100) 
     1169                END IF 
     1170#endif 
     1171                   ! 
     1172                IF ( dh(ji,jj)  <  0.2*hbl(ji,jj) .AND. ibld(ji,jj) - imld(ji,jj) > 3 ) THEN 
    11061173                   zbuoy_pyc_sc = 2.0_wp * MAX(zdb_ml(ji,jj), 0._wp) / zdh(ji,jj) 
    11071174                   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 ) ) 
    … …  
    11181185                      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 
    11191186                   END DO 
     1187#ifdef key_osm_debug 
     1188                   IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     1189                      WRITE(narea+100,'(2(a,g11.3))')'lpyc= lconv=T,dh<0.2*hbl: zbuoy_pyc_sc=',zbuoy_pyc_sc,' zdelta_pyc=',zdelta_pyc 
     1190                      WRITE(narea+100,'(3(a,g11.3))')'zwt_pyc_sc_1=',zwt_pyc_sc_1,' zws_pyc_sc_1=',zws_pyc_sc_1,' zzeta_pyc=',zzeta_pyc 
     1191                      FLUSH(narea+100) 
     1192                   END IF 
     1193#endif 
     1194 
    11201195                END IF 
    11211196             ENDIF ! End of pycnocline                   
    … …  
    11291204    END DO      ! jj loop 
    11301205 
    1131     WHERE ( lconv ) 
    1132        zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 
    1133        zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0) 
    1134        zsc_vw_1 = 0._wp 
     1206    WHERE ( lconv(2:jpim1,2:jpjm1) ) 
     1207       zsc_uw_1(2:jpim1,2:jpjm1) = -zwb0(2:jpim1,2:jpjm1) * zustar(2:jpim1,2:jpjm1)**2 * zhml(2:jpim1,2:jpjm1) / & 
     1208            & ( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln ) 
     1209       zsc_uw_2(2:jpim1,2:jpjm1) =  zwb0(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1)    * zhml(2:jpim1,2:jpjm1) / & 
     1210            & ( zvstr(2:jpim1,2:jpjm1)**3 + 0.5 * zwstrc(2:jpim1,2:jpjm1)**3 + epsln )**(2.0/3.0) 
     1211       zsc_vw_1(2:jpim1,2:jpjm1) = 0._wp 
    11351212    ELSEWHERE 
    1136        zsc_uw_1 = 0._wp 
    1137        zsc_vw_1 = 0._wp 
     1213       zsc_uw_1(2:jpim1,2:jpjm1) = 0._wp 
     1214       zsc_vw_1(2:jpim1,2:jpjm1) = 0._wp 
    11381215    ENDWHERE 
    11391216 
    … …  
    12771354    END DO      ! jj loop 
    12781355 
    1279     WHERE ( lconv ) 
    1280        zsc_uw_1 = zustar**2 
    1281        zsc_vw_1 = ff_t * zustke * zhml 
     1356    WHERE ( lconv(2:jpim1,2:jpjm1) ) 
     1357       zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2 
     1358       zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1) * zhml(2:jpim1,2:jpjm1) 
    12821359    ELSEWHERE 
    1283        zsc_uw_1 = zustar**2 
    1284        zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1 
    1285        zsc_vw_1 = ff_t * zustke * zhbl 
    1286        zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1 
     1360       zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2 
     1361       zsc_uw_2(2:jpim1,2:jpjm1) = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1(2:jpim1,2:jpjm1) 
     1362       zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1) * zhbl(2:jpim1,2:jpjm1) 
     1363       zsc_vw_2(2:jpim1,2:jpjm1) = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1(2:jpim1,2:jpjm1) 
    12871364    ENDWHERE 
    12881365 
    … …  
    13261403       ji=iloc_db; jj=jloc_db 
    13271404       jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) 
    1328        WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/contrib to ghamt/s:  zsc_wth_1=',zsc_wth_1(ji,jj), '  zsc_ws_1=',zsc_ws_1(ji,jj) 
     1405       WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/s:  zsc_wth_1=',zsc_wth_1(ji,jj), '  zsc_ws_1=',zsc_ws_1(ji,jj) 
    13291406       IF (lpyc(ji,jj)) WRITE(narea+100,'(2(a,g11.3))') 'zsc_wth_pyc=', zsc_wth_pyc(ji,jj), '  zsc_wth_pyc=',zsc_wth_pyc(ji,jj) 
    13301407       WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) 
    … …  
    17111788 
    17121789      ! Scales used to calculate eddy diffusivity and viscosity profiles 
    1713       REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 
    1714       REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 
    1715       REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 
    1716       REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 
     1790      REAL(wp), DIMENSION(jpi,jpj) :: pdifml_sc, pvisml_sc 
     1791      REAL(wp), DIMENSION(jpi,jpj) :: pdifpyc_n_sc, pdifpyc_s_sc, pdifpyc_shr 
     1792      REAL(wp), DIMENSION(jpi,jpj) :: pvispyc_n_sc, pvispyc_s_sc,pvispyc_shr 
     1793      REAL(wp), DIMENSION(jpi,jpj) :: pbeta_d_sc, pbeta_v_sc 
     1794      REAL(wp), DIMENSION(jpi,jpj) :: pb_coup, pc_coup_vis, pc_coup_dif 
    17171795      ! 
    1718       REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 
     1796      REAL(wp) :: pvel_sc_pyc, pvel_sc_ml, pstab_fac, pz_b 
     1797      REAL(wp) :: pa_cubic, pb_cubic, pc_cubic, pd_cubic 
     1798      REAL(wp) :: pznd_ml, pznd_pyc 
     1799      REAL(wp) :: zmsku, zmskv 
    17191800 
    17201801      REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 
    … …  
    17241805      DO jj = 2, jpjm1 
    17251806         DO ji = 2, jpim1 
     1807            pb_coup(ji,jj) = 0._wp 
    17261808            IF ( lconv(ji,jj) ) THEN 
    17271809 
    1728                zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 
    1729                zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    1730                zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 
    1731  
    1732                zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 
    1733                zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 
     1810               pvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 
     1811               pvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     1812               pstab_fac = ( zhml(ji,jj) / pvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 
     1813 
     1814               pdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * pvel_sc_ml 
     1815               pvisml_sc(ji,jj) = rn_vis_ml * pdifml_sc(ji,jj) 
    17341816#ifdef key_osm_debug 
    17351817               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1736                   WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) 
    1737                   WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac 
     1818                  WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',pdifml_sc(ji,jj),' zvisml_sc=',pvisml_sc(ji,jj) 
     1819                  WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',pvel_sc_pyc,' zvel_sc_ml=',pvel_sc_ml,' pstab_fac=',pstab_fac 
    17381820                  FLUSH(narea+100) 
    17391821               END IF 
    17401822#endif 
    17411823               IF ( lpyc(ji,jj) ) THEN 
    1742                   zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 
    1743                   zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 
    1744                   zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 
     1824                  pdifpyc_n_sc(ji,jj) =  rn_dif_pyc * pvel_sc_ml * zdh(ji,jj) 
     1825                  pvispyc_n_sc(ji,jj) = 0.09 * pvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 
     1826                  pvispyc_n_sc(ji,jj) = rn_vis_pyc * pvel_sc_ml * zdh(ji,jj) + pvispyc_n_sc(ji,jj) * pstab_fac 
    17451827#ifdef key_osm_debug 
    17461828                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1747                      WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 
     1829                     WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj) 
    17481830                     FLUSH(narea+100) 
    17491831                  END IF 
    17501832#endif 
    17511833                  IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 
    1752                      zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 
    1753                      zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 
     1834                     pdifpyc_n_sc(ji,jj) = pdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 
     1835                     pvispyc_n_sc(ji,jj) = pvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 
    17541836                  ENDIF 
    17551837#ifdef key_osm_debug 
    17561838                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1757                      WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 
     1839                     WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj) 
    17581840                     FLUSH(narea+100) 
    17591841                  END IF 
    17601842#endif 
    1761                   zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 
    1762                   zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 
     1843                  pdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * pvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 
     1844                  pvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * pvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 
    17631845#ifdef key_osm_debug 
    17641846                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1765                      WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 
     1847                     WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj) 
    17661848                     FLUSH(narea+100) 
    17671849                  END IF 
    17681850#endif 
    1769                   zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 
    1770                   zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 
     1851                  pdifpyc_s_sc(ji,jj) = 0.09 * pdifpyc_s_sc(ji,jj) * pstab_fac 
     1852                  pvispyc_s_sc(ji,jj) = pvispyc_s_sc(ji,jj) * pstab_fac 
    17711853#ifdef key_osm_debug 
    17721854                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1773                      WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 
     1855                     WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj) 
    17741856                     FLUSH(narea+100) 
    17751857                  END IF 
    17761858#endif 
    17771859 
    1778                   zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 
    1779                   zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_n_sc(ji,jj) ) 
    1780 #ifdef key_osm_debug 
    1781                   IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1782                      WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 
    1783                      FLUSH(narea+100) 
    1784                   END IF 
    1785 #endif 
    1786  
    1787                   zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third 
    1788                   zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 
     1860                  pdifpyc_s_sc(ji,jj) = MAX( pdifpyc_s_sc(ji,jj), -0.5 * pdifpyc_n_sc(ji,jj) ) 
     1861                  pvispyc_s_sc(ji,jj) = MAX( pvispyc_s_sc(ji,jj), -0.5 * pvispyc_n_sc(ji,jj) ) 
     1862 
     1863                  pbeta_d_sc(ji,jj) = 1.0 - ( ( pdifpyc_n_sc(ji,jj) + 1.4 * pdifpyc_s_sc(ji,jj) ) / ( pdifml_sc(ji,jj) + epsln ) )**p2third 
     1864                  pbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( pvispyc_n_sc(ji,jj) + pvispyc_s_sc(ji,jj) ) / ( pvisml_sc(ji,jj) + epsln ) 
    17891865               ELSE 
    1790                   zbeta_d_sc(ji,jj) = 1.0 
    1791                   zbeta_v_sc(ji,jj) = 1.0 
    1792                ENDIF 
     1866                  pdifpyc_n_sc(ji,jj) =  rn_dif_pyc * pvel_sc_ml * zdh(ji,jj)    ! ag 19/03 
     1867                  pdifpyc_s_sc(ji,jj) = 0._wp   ! ag 19/03 
     1868                  pvispyc_n_sc(ji,jj) = rn_vis_pyc * pvel_sc_ml * zdh(ji,jj)  ! ag 19/03 
     1869                  pvispyc_s_sc(ji,jj) = 0._wp  ! ag 19/03 
     1870                  IF(lcoup(ji,jj) ) THEN   ! ag 19/03 
     1871                    ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = Cd|ub| 
     1872                    !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90 
     1873                    !  Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub) 
     1874                    ! wet-cell averaging .. 
     1875                    zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     1876                    zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     1877                    pb_coup(ji,jj) = 0.4 * SQRT(rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
     1878                  &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) ) 
     1879 
     1880                    pz_b = -gdepw_n(ji,jj,mbkt(ji,jj) + 1 )  ! ag 19/03 
     1881                    pc_coup_vis(ji,jj) = -0.5 * ( 0.5 * pvisml_sc(ji,jj) / zhml(ji,jj) - pb_coup(ji,jj) ) / ( zhml(ji,jj) + pz_b )  ! ag 19/03 
     1882#ifdef key_osm_debug 
     1883                    IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     1884                       WRITE(narea+100,'(3(a,g11.3))')' lcoup = T; 1st pz_b= ', pz_b, 'pb_coup', pb_coup(ji,jj), ' pc_coup_vis', pc_coup_vis(ji,jj) 
     1885                       FLUSH(narea+100) 
     1886                    END IF 
     1887#endif 
     1888#ifdef key_osm_debug 
     1889                    WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), 'jjg= ', mjg(jj) 
     1890                    WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', pz_b, 'pb_coup', pb_coup(ji,jj), ' pc_coup_vis', pc_coup_vis(ji,jj) 
     1891                    FLUSH(narea+400) 
     1892#endif 
     1893                    pz_b = -zhml(ji,jj) + gdepw_n(ji,jj,mbkt(ji,jj)+1)  ! ag 19/03  
     1894                    pbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b**2 ) / pvisml_sc(ji,jj)  ! ag 19/03 
     1895                    pbeta_d_sc(ji,jj) = 1.0 - ( ( pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b**2 ) / pdifml_sc(ji,jj) )**p2third 
     1896                    pc_coup_dif(ji,jj) = 0.5 * ( -pdifml_sc(ji,jj) / zhml(ji,jj) * ( 1.0 - pbeta_d_sc(ji,jj) )**1.5 + 1.5 * (pdifml_sc(ji,jj) / zhml(ji,jj) )* pbeta_d_sc(ji,jj) * SQRT( 1.0 - pbeta_d_sc(ji,jj) ) -pb_coup(ji,jj) ) / pz_b  ! ag 19/03 
     1897#ifdef key_osm_debug 
     1898                    IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     1899                       WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', pz_b, ' pc_coup_dif', pc_coup_dif(ji,jj) 
     1900                       FLUSH(narea+100) 
     1901                    END IF 
     1902#endif 
     1903#ifdef key_osm_debug 
     1904                    WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_vis', pc_coup_vis(ji,jj) 
     1905                    FLUSH(narea+400) 
     1906#endif 
     1907                  ELSE     ! ag 19/03 
     1908                    pbeta_d_sc(ji,jj) = 1.0 - ( ( pdifpyc_n_sc(ji,jj) + 1.4 * pdifpyc_s_sc(ji,jj) ) / ( pdifml_sc(ji,jj) + epsln ) )**p2third     ! ag 19/03 
     1909                    pbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( pvispyc_n_sc(ji,jj) + pvispyc_s_sc(ji,jj) ) / ( pvisml_sc(ji,jj) + epsln ) ! ag 19/03 
     1910                  ENDIF ! ag 19/03 
     1911               ENDIF    ! ag 19/03 
    17931912#ifdef key_osm_debug 
    17941913               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1795                   WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj) 
     1914                  WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',pbeta_d_sc(ji,jj) ,' zbeta_v_sc=',pbeta_v_sc(ji,jj) 
     1915                  WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj) 
     1916                  WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj) 
    17961917                  FLUSH(narea+100) 
    17971918               END IF 
    17981919#endif 
    17991920            ELSE ! conv, stable 
    1800                zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
    1801                zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1921               pdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
     1922               pvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 
    18021923#ifdef key_osm_debug 
    18031924               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1804                   WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) 
     1925                  WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',pdifml_sc(ji,jj),' zvisml_sc=',pvisml_sc(ji,jj) 
    18051926                  FLUSH(narea+100) 
    18061927               END IF 
    … …  
    18151936            IF ( lconv(ji,jj) ) THEN 
    18161937               DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
    1817                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
     1938                  pznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    18181939                  ! 
    1819                   zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 
     1940                  zdiffut(ji,jj,jk) = pdifml_sc(ji,jj) * pznd_ml * ( 1.0 - pbeta_d_sc(ji,jj) * pznd_ml )**1.5 
    18201941                  ! 
    1821                   zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 
    1822                        &            *                                      ( 1.0 - 0.5 * zznd_ml**2 ) 
     1942                  zviscos(ji,jj,jk) = pvisml_sc(ji,jj) * pznd_ml * ( 1.0 - pbeta_v_sc(ji,jj) * pznd_ml ) & 
     1943                       &            *                                      ( 1.0 - 0.5 * pznd_ml**2 ) 
    18231944               END DO 
     1945                
     1946! Coupling to bottom 
     1947             
     1948               IF ( lcoup(ji,jj) ) THEN                                                          ! ag 19/03 
     1949                  DO jk = mbkt(ji,jj), imld(ji,jj), -1                                           ! ag 19/03 
     1950                     pz_b = - ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mbkt(ji,jj) + 1 ) )            ! ag 19/03 
     1951                     zviscos(ji,jj,jk) = pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b**2    ! ag 19/03 
     1952                     zdiffut(ji,jj,jk) = pb_coup(ji,jj) * pz_b + pc_coup_dif(ji,jj) * pz_b**2    ! ag 19/03 
     1953                  END DO                                                                         ! ag 19/03 
     1954               ENDIF                                                                             ! ag 19/03 
    18241955               ! pycnocline 
    18251956               IF ( lpyc(ji,jj) ) THEN 
    1826                   ! Diffusivity profile in the pycnocline given by cubic polynomial. 
    1827                   za_cubic = 0.5 
    1828                   zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 
    1829                   zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & 
    1830                        & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 
    1831                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic  - zb_cubic ) 
    1832                   zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1957                  ! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed. 
     1958                  pa_cubic = 0.5 
     1959                  pb_cubic = -1.75 * pdifpyc_s_sc(ji,jj) / pdifpyc_n_sc(ji,jj) 
     1960                  pd_cubic = ( zdh(ji,jj) * pdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - pbeta_d_sc(ji,jj) ) * ( 2.5 * pbeta_d_sc(ji,jj) - 1.0 ) & 
     1961                       & - 0.85 * pdifpyc_s_sc(ji,jj) ) / MAX(pdifpyc_n_sc(ji,jj), 1.e-8) 
     1962                  pd_cubic = pd_cubic - pb_cubic - 2.0 * ( 1.0 - pa_cubic  - pb_cubic ) 
     1963                  pc_cubic = 1.0 - pa_cubic - pb_cubic - pd_cubic 
    18331964                  DO jk = imld(ji,jj) , ibld(ji,jj) 
    1834                      zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1965                     pznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
    18351966                     ! 
    1836                      zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 +   zd_cubic * zznd_pyc**3 ) 
    1837  
    1838                      zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) 
     1967                     zdiffut(ji,jj,jk) = pdifpyc_n_sc(ji,jj) * ( pa_cubic + pb_cubic * pznd_pyc + pc_cubic * pznd_pyc**2 +   pd_cubic * pznd_pyc**3 ) 
     1968 
     1969                     zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + pdifpyc_s_sc(ji,jj) * ( 1.75 * pznd_pyc - 0.15 * pznd_pyc**2 - 0.2 * pznd_pyc**3 ) 
    18391970                  END DO 
    18401971                  ! viscosity profiles. 
    1841                   za_cubic = 0.5 
    1842                   zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
    1843                   zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 
    1844                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 
    1845                   zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
     1972                  pa_cubic = 0.5 
     1973                  pb_cubic = -1.75 * pvispyc_s_sc(ji,jj) / pvispyc_n_sc(ji,jj) 
     1974                  pd_cubic = ( 0.5 * pvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * pvispyc_s_sc(ji,jj)  )  / MAX(pvispyc_n_sc(ji,jj), 1.e-8) 
     1975                  pd_cubic = pd_cubic - pb_cubic - 2.0 * ( 1.0 - pa_cubic - pb_cubic ) 
     1976                  pc_cubic = 1.0 - pa_cubic - pb_cubic - pd_cubic 
    18461977                  DO jk = imld(ji,jj) , ibld(ji,jj) 
    1847                      zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
    1848                      zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 
    1849                      zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) 
     1978                     pznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 
     1979                     zviscos(ji,jj,jk) = pvispyc_n_sc(ji,jj) * ( pa_cubic + pb_cubic * pznd_pyc + pc_cubic * pznd_pyc**2 + pd_cubic * pznd_pyc**3 ) 
     1980                     zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + pvispyc_s_sc(ji,jj) * ( 1.75 * pznd_pyc - 0.15 * pznd_pyc**2 -0.2 * pznd_pyc**3 ) 
    18501981                  END DO 
    1851                   IF ( zdhdt(ji,jj) > 0._wp ) THEN 
    1852                      zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
    1853                      zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
    1854                   ELSE 
    1855                      zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
    1856                      zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     1982!                  IF ( zdhdt(ji,jj) > 0._wp ) THEN 
     1983!                     zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
     1984!                     zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 
     1985!                  ELSE 
     1986!                     zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 
     1987!                     zviscos(ji,jj,ibld(ji,jj)) = 0._wp 
     1988!                  ENDIF 
     1989               ELSE 
     1990! lpyc set false but not coupled to the bottom. 
     1991                  IF ( .not. lcoup(ji,jj) ) THEN 
     1992                     zdiffut(ji,jj,ibld(ji,jj)) = 0.5 * pdifpyc_n_sc(ji,jj) 
     1993                     zviscos(ji,jj,ibld(ji,jj)) = 0.5 * pvispyc_n_sc(ji,jj) 
    18571994                  ENDIF 
    18581995               ENDIF 
    … …  
    18601997               ! stable conditions 
    18611998               DO jk = 2, ibld(ji,jj) 
    1862                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1863                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
    1864                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
     1999                  pznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     2000                  zdiffut(ji,jj,jk) = 0.75 * pdifml_sc(ji,jj) * pznd_ml * ( 1.0 - pznd_ml )**1.5 
     2001                  zviscos(ji,jj,jk) = 0.375 * pvisml_sc(ji,jj) * pznd_ml * (1.0 - pznd_ml) * ( 1.0 - pznd_ml**2 ) 
    18652002               END DO 
    18662003 
    … …  
    18732010         END DO  ! end of ji loop 
    18742011      END DO     ! end of jj loop 
     2012      IF ( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1)*pb_coup(:,:) )         ! BBL-coupling velocity scale 
    18752013 
    18762014    END SUBROUTINE zdf_osm_diffusivity_viscosity 
    … …  
    19222060      END DO 
    19232061 
    1924       zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 
    1925  
    1926       zshear(:,:) = 0._wp 
     2062      zekman(2:jpim1,2:jpjm1) = EXP( - zek * ABS( ff_t(2:jpim1,2:jpjm1) ) * zhbl(2:jpim1,2:jpjm1) / MAX(zustar(2:jpim1,2:jpjm1), 1.e-8 ) ) 
     2063 
     2064      zshear(2:jpim1,2:jpjm1) = 0._wp 
    19272065#ifdef key_osm_debug 
    19282066      IF(narea==nn_narea_db) THEN 
    … …  
    19332071      END IF 
    19342072#endif 
    1935       j_ddh(:,:) = 1      
     2073      j_ddh(2:jpim1,2:jpjm1) = 1      
    19362074 
    19372075      DO jj = 2, jpjm1 
    … …  
    20372175      END DO   ! jj loop  
    20382176 
    2039       zwb_min(:,:) = 0._wp 
     2177      zwb_min(:,:) = zlarge 
    20402178 
    20412179      DO jj = 2, jpjm1 
    … …  
    21162254 
    21172255 
    2118       zt   = 0._wp 
    2119       zs   = 0._wp 
    2120       zu   = 0._wp 
    2121       zv   = 0._wp 
     2256      zt(2:jpim1,2:jpjm1)  = 0._wp 
     2257      zs(2:jpim1,2:jpjm1)  = 0._wp 
     2258      zu(2:jpim1,2:jpjm1)  = 0._wp 
     2259      zv(2:jpim1,2:jpjm1)  = 0._wp 
    21222260      DO jj = 2, jpjm1                                 !  Vertical slab 
    21232261         DO ji = 2, jpim1 
    … …  
    21432281            zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 
    21442282            ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 
    2145             IF ( ibld_ext < mbkt(ji,jj) ) THEN 
     2283            IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN      ! ag 09/03 
     2284 ! Two external levels are available. 
    21462285               zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem) 
    21472286               zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal) 
    … …  
    22132352      REAL(wp)                      :: zpe_mle_ref, zdbdz_mle_int 
    22142353 
    2215       znd_param(:,:) = 0._wp 
     2354      znd_param(2:jpim1,2:jpjm1) = 0._wp 
    22162355 
    22172356      DO jj = 2, jpjm1 
Note: See TracChangeset for help on using the changeset viewer.