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

Changeset 14305


Ignore:
Timestamp:
2021-01-14T19:43:03+01:00 (2 years ago)
Author:
smueller
Message:

Reduction of memory usage by subroutine zdf_osm 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

    r14280 r14305  
    113113   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value 
    114114   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la 
     115   LOGICAL  ::   ln_dia_pyc_scl = .FALSE.   ! Output of pycnocline scalar-gradient profiles 
     116   LOGICAL  ::   ln_dia_pyc_shr = .FALSE.   ! Output of pycnocline velocity-shear  profiles 
    115117 
    116118 
     
    215217      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points) 
    216218      !! 
    217       INTEGER ::   ji, jj, jk                   ! dummy loop indices 
     219      INTEGER ::   ji, jj, jk, jkflt                               ! dummy loop indices 
    218220 
    219221      INTEGER ::   jl                   ! dummy loop indices 
    220222 
    221       INTEGER ::   ikbot, jkmax, jkm1, jkp2     ! 
     223      INTEGER ::   ikbot, jkm1, jkp2     ! 
    222224 
    223225      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      ! 
     
    233235      REAL(wp) :: zt,zs,zu,zv,zrh               ! variables used in constructing averages 
    234236! Scales 
    235       REAL(wp), DIMENSION(jpi,jpj) :: zrad0     ! Surface solar temperature flux (deg m/s) 
    236       REAL(wp), DIMENSION(jpi,jpj) :: zradh     ! Radiative flux at bl base (Buoyancy units) 
    237       REAL(wp), DIMENSION(jpi,jpj) :: zradav    ! Radiative flux, bl average (Buoyancy Units) 
     237      REAL(wp)   :: zrad0                       ! Surface solar temperature flux (deg m/s) 
     238      REAL(wp)   :: zradh                       ! Radiative flux at bl base (Buoyancy units) 
     239      REAL(wp)   :: zradav                      ! Radiative flux, bl average (Buoyancy Units) 
    238240      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity 
    239241      REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale 
     
    241243      REAL(wp), DIMENSION(jpi,jpj) :: zwstrc    ! Convective velocity scale 
    242244      REAL(wp), DIMENSION(jpi,jpj) :: zuw0      ! Surface u-momentum flux 
    243       REAL(wp), DIMENSION(jpi,jpj) :: zvw0      ! Surface v-momentum flux 
     245      REAL(wp)   :: zvw0                        ! Surface v-momentum flux 
    244246      REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic) 
    245247      REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux 
     
    299301      REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 
    300302      REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline 
    301       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline 
    302       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline 
     303      REAL(wp) ::   zdtdz_pyc    ! Parametrized gradient of temperature in pycnocline 
     304      REAL(wp) ::   zdsdz_pyc    ! Parametrised gradient of salinity in pycnocline 
    303305      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline 
    304       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline 
    305       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline 
     306      REAL(wp) ::  zdudz_pyc    ! u-shear across the pycnocline 
     307      REAL(wp) ::  zdvdz_pyc    ! v-shear across the pycnocline 
    306308      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient. 
    307309      ! Flux-gradient relationship variables 
     
    317319 
    318320      ! For calculating Ri#-dependent mixing 
    319       REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du   ! u-shear^2 
    320       REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv   ! v-shear^2 
    321       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion 
     321      REAL(wp), DIMENSION(jpi,jpj) ::   z2du     ! u-shear^2 
     322      REAL(wp), DIMENSION(jpi,jpj) ::   z2dv     ! v-shear^2 
     323      REAL(wp)                     ::   zrimix   ! Spatial form of ri#-induced diffusion 
    322324 
    323325      ! Temporary variables 
     
    342344      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 
    343345      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 
    344348 
    345349      ! For debugging 
     
    349353      IF( ln_timing ) CALL timing_start('zdf_osm') 
    350354      ibld(:,:)   = 0     ; imld(:,:)  = 0 
    351       zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp 
     355      zustar(:,:) = 0.0_wp 
    352356      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp 
    353       zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp 
     357      zwth0(:,:)  = 0.0_wp ; zws0(:,:)  = 0.0_wp ; zwb0(:,:)      = 0.0_wp 
    354358      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp 
    355359      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 
     
    370374      zwth_ent = 0._wp ; zws_ent = 0._wp 
    371375      ! 
    372       zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 
    373       zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 
     376      zdbdz_pyc(:,:,:) = 0.0_wp 
    374377      ! 
    375378      zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 
     
    396399      ! Calculate boundary layer scales 
    397400      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    398  
    399       ! Assume two-band radiation model for depth of OSBL 
    400      zz0 =       rn_abs       ! surface equi-partition in 2-bands 
    401      zz1 =  1. - rn_abs 
     401      ! 
     402      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 
     403      zz0 =           rn_abs   ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands 
     404      zz1 =  1.0_wp - rn_abs 
     405      DO_2D( 0, 0, 0, 0 ) 
     406         zrad0         = qsr(ji,jj) * r1_rho0_rcp                          ! Surface downward irradiance (so always +ve) 
     407         zradh         = zrad0 *                                       &   ! Downwards irradiance at base of boundary layer 
     408            &            ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) ) 
     409         zradav        = zrad0 *                                       &   ! Downwards irradiance averaged over depth of the OSBL 
     410            &     ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 +   & 
     411            &       zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) 
     412         zwth0(ji,jj)  = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)       ! Upwards surface Temperature flux for non-local term 
     413         zwthav(ji,jj) = 0.5_wp * zwth0(ji,jj) -                       &   ! Turbulent heat flux averaged over depth of OSBL 
     414            &            ( 0.5_wp * ( zrad0 + zradh ) - zradav ) 
     415     END_2D 
    402416     DO_2D( 0, 0, 0, 0 ) 
    403         ! Surface downward irradiance (so always +ve) 
    404         zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 
    405         ! Downwards irradiance at base of boundary layer 
    406         zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 
    407         ! Downwards irradiance averaged over depth of the OSBL 
    408         zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 
    409               &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 
     417        zws0(ji,jj)    = -1.0_wp *                                     &   ! Upwards surface salinity flux for non-local term 
     418           &          ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 
     419        zthermal       = rab_n(ji,jj,1,jp_tem) 
     420        zbeta          = rab_n(ji,jj,1,jp_sal) 
     421        zwb0(ji,jj)    = grav * zthermal * zwth0(ji,jj) -              &   ! Non radiative upwards surface buoyancy flux 
     422           &             grav * zbeta * zws0(ji,jj) 
     423        zwsav(ji,jj)   = 0.5 * zws0(ji,jj)                                 ! Turbulent salinity flux averaged over depth of the OBSL 
     424        zwbav(ji,jj)   = grav  * zthermal * zwthav(ji,jj) -            &   ! Turbulent buoyancy flux averaged over the depth of the 
     425           &             grav  * zbeta * zwsav(ji,jj)                      ! OBSBL 
    410426     END_2D 
    411      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 
    412427     DO_2D( 0, 0, 0, 0 ) 
    413         zthermal = rab_n(ji,jj,1,jp_tem) 
    414         zbeta    = rab_n(ji,jj,1,jp_sal) 
    415         ! Upwards surface Temperature flux for non-local term 
    416         zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 
    417         ! Upwards surface salinity flux for non-local term 
    418         zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 
    419         ! Non radiative upwards surface buoyancy flux 
    420         zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
    421         ! turbulent heat flux averaged over depth of OSBL 
    422         zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 
    423         ! turbulent salinity flux averaged over depth of the OBSL 
    424         zwsav(ji,jj) = 0.5 * zws0(ji,jj) 
    425         ! turbulent buoyancy flux averaged over the depth of the OBSBL 
    426         zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    427         ! Surface upward velocity fluxes 
    428         zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
    429         zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
    430         ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    431         zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
     428        zuw0(ji,jj)    = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) *       &   ! Surface upward velocity fluxes 
     429           &             r1_rho0 * tmask(ji,jj,1) 
     430        zvw0           = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
     431        zustar(ji,jj)  = MAX( SQRT( SQRT( zuw0(ji,jj) *                &   ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
     432           &                              zuw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp ) 
    432433        zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
    433         zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
     434        zsin_wind(ji,jj) = -zvw0        / ( zustar(ji,jj) * zustar(ji,jj) ) 
    434435     END_2D 
    435436     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) 
     
    699700 
    700701      CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    701       CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 
    702       CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 
     702      CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_pyc ) 
    703703       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    704704       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
     
    10281028      END_2D 
    10291029 
    1030       ! pynocline contributions 
    1031        DO_2D( 0, 0, 0, 0 ) 
    1032          IF ( .not. lconv(ji,jj) ) THEN 
    1033           IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    1034              DO jk= 2, ibld(ji,jj) 
    1035                 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    1036                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1037                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1038                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1039                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
    1040              END DO 
    1041           END IF 
     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 
    10421101         END IF 
    1043        END_2D 
     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 
    10441135      IF(ln_dia_osm) THEN 
    1045           IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 
    1046           IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 
    1047        END IF 
     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 
    10481142 
    10491143       DO_2D( 0, 0, 0, 0 ) 
     
    10571151          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 
    10581152          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 
    1059           IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 
    1060           IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 
    10611153          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 
    10621154       END IF 
     1155 
    10631156       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    10641157       ! Need to put in code for contributions that are applied explicitly to 
     
    10801173       END_2D 
    10811174 
    1082        IF(ln_dia_osm) THEN 
    1083           IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
    1084           IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
    1085           IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
    1086        END IF 
    1087  
    1088 ! KPP-style Ri# mixing 
    1089        IF( ln_kpprimix) THEN 
    1090           DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form) 
    1091              z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   & 
    1092                   &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 
    1093                   &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 
    1094              z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   & 
    1095                   &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) & 
    1096                   &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 
    1097           END_3D 
    1098       ! 
    1099          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    1100             !                                          ! shear prod. at w-point weightened by mask 
    1101             zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    1102                &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
    1103             !                                          ! local Richardson number 
    1104             zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 
    1105             zfri =  MIN( zri / rn_riinfty , 1.0_wp ) 
    1106             zfri  = ( 1.0_wp - zfri * zfri ) 
    1107             zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk) 
    1108          END_3D 
    1109  
     1175      ! KPP-style Ri# mixing 
     1176      IF ( ln_kpprimix ) THEN 
     1177          jkflt = jpk 
    11101178          DO_2D( 0, 0, 0, 0 ) 
    1111              DO jk = ibld(ji,jj) + 1, jpkm1 
    1112                 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
    1113                 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
    1114              END DO 
     1179             IF ( ibld(ji,jj) < jkflt ) jkflt = ibld(ji,jj) 
    11151180          END_2D 
    1116  
     1181          DO jk = jkflt+1, jpkm1 
     1182             ! Shear production at uw- and vw-points (energy conserving form) 
     1183             DO_2D( 1, 0, 1, 0 ) 
     1184                IF ( jk > MIN( ibld(ji,jj), ibld(ji+1,jj) ) ) THEN 
     1185                   z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) *                      & 
     1186                      &                   ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) /   & 
     1187                      &                   ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 
     1188                END IF 
     1189                IF ( jk > MIN( ibld(ji,jj), ibld(ji,jj+1) ) ) THEN 
     1190                   z2dv(ji,jj) = 0.5_wp * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) *                      & 
     1191                      &                   ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) /   & 
     1192                      &                   ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 
     1193                END IF 
     1194             END_2D 
     1195             DO_2D( 0, 0, 0, 0 ) 
     1196                IF ( jk > ibld(ji,jj) ) THEN 
     1197                   ! Shear prod. at w-point weightened by mask 
     1198                   zesh2  =  ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     1199                      &    + ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
     1200                   ! Local Richardson number 
     1201                   zri   = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX(zesh2, epsln) 
     1202                   zfri =  MIN( zri / rn_riinfty , 1.0_wp ) 
     1203                   zfri  = ( 1.0_wp - zfri * zfri ) 
     1204                   zrimix  =  zfri * zfri  * zfri * wmask(ji, jj, jk) 
     1205                   zdiffut(ji,jj,jk) = zrimix*rn_difri 
     1206                   zviscos(ji,jj,jk) = zrimix*rn_difri 
     1207                END IF 
     1208             END_2D 
     1209          END DO 
    11171210       END IF ! ln_kpprimix = .true. 
    11181211 
     
    11571250          END_2D 
    11581251       ENDIF 
    1159  
    1160        IF(ln_dia_osm) THEN 
    1161           IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 
    1162           IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 
    1163           IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 
    1164        END IF 
    1165  
    11661252 
    11671253       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 
     
    17271813    END SUBROUTINE zdf_osm_external_gradients 
    17281814 
    1729     SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 
    1730  
    1731      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline 
    1732      REAL(wp), DIMENSION(jpi,jpj) :: zalpha 
    1733  
    1734      INTEGER :: jk, jj, ji 
    1735      REAL(wp) :: ztgrad, zsgrad, zbgrad 
    1736      REAL(wp) :: zgamma_b_nd, znd 
    1737      REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 
    1738      REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 
    1739  
    1740      IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 
    1741      DO_2D( 0, 0, 0, 0 ) 
    1742         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    1743            IF ( lconv(ji,jj) ) THEN  ! convective conditions 
    1744              IF ( lpyc(ji,jj) ) THEN 
    1745                 zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
    1746                 zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 
    1747                 zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 
    1748  
    1749                 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
    1750 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    1751 ! Commented lines in this section are not needed in new code, once tested ! 
    1752 ! can be removed                                                          ! 
    1753 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    1754 !                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
    1755 !                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
    1756                 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 
    1757                 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
    1758                 DO jk = 2, ibld(ji,jj)+ibld_ext 
    1759                   znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 
    1760                   IF ( znd <= zzeta_m ) THEN 
    1761 !                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 
    1762 !                &        EXP( -6.0 * ( znd -zzeta_m )**2 ) 
    1763 !                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 
    1764 !                                  & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
    1765                      zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 
    1766                                & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
    1767                   ELSE 
    1768 !                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
    1769 !                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
    1770                       zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
    1771                   ENDIF 
    1772                END DO 
    1773             ENDIF ! if no pycnocline pycnocline gradients set to zero 
    1774            ELSE 
    1775               ! stable conditions 
    1776               ! if pycnocline profile only defined when depth steady of increasing. 
    1777               IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady. 
    1778                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    1779                     IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    1780                        ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 
    1781                        ztgrad = zdt_bl(ji,jj) * ztmp 
    1782                        zsgrad = zds_bl(ji,jj) * ztmp 
    1783                        zbgrad = zdb_bl(ji,jj) * ztmp 
    1784                        DO jk = 2, ibld(ji,jj) 
    1785                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
    1786                           zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    1787                           zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    1788                           zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    1789                        END DO 
    1790                     ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    1791                        ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 
    1792                        ztgrad = zdt_bl(ji,jj) * ztmp 
    1793                        zsgrad = zds_bl(ji,jj) * ztmp 
    1794                        zbgrad = zdb_bl(ji,jj) * ztmp 
    1795                        DO jk = 2, ibld(ji,jj) 
    1796                           znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 
    1797                           zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    1798                           zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    1799                           zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    1800                        END DO 
    1801                     ENDIF ! IF (zhol >=0.5) 
    1802                  ENDIF    ! IF (zdb_bl> 0.) 
    1803               ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
    1804            ENDIF          ! IF (lconv) 
    1805         ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 
    1806      END_2D 
    1807      IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 
    1808  
    1809     END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 
    1810  
    1811     SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 
    1812       !!--------------------------------------------------------------------- 
    1813       !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  *** 
    1814       !! 
    1815       !! ** Purpose : Calculates velocity shear in the pycnocline 
    1816       !! 
    1817       !! ** Method  : 
    1818       !! 
    1819       !!---------------------------------------------------------------------- 
    1820  
    1821       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 
    1822  
    1823       INTEGER :: jk, jj, ji 
    1824       REAL(wp) :: zugrad, zvgrad, znd 
    1825       REAL(wp) :: zzeta_v = 0.45 
     1815   SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha ) 
     1816      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   pdbdz   ! Gradients in the pycnocline 
     1817      REAL(wp), DIMENSION(:,:),   INTENT( inout ) ::   palpha 
     1818      INTEGER                                     ::   jk, jj, ji 
     1819      REAL(wp)                                    ::   zbgrad 
     1820      REAL(wp)                                    ::   zgamma_b_nd, znd 
     1821      REAL(wp)                                    ::   zzeta_m 
     1822      REAL(wp), PARAMETER                         ::   ppgamma_b = 2.25_wp 
    18261823      ! 
    1827       IF( ln_timing ) CALL timing_start('zdf_osm_pshp') 
     1824      IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 
     1825      ! 
    18281826      DO_2D( 0, 0, 0, 0 ) 
    1829          ! 
    18301827         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    1831             IF ( lconv (ji,jj) ) THEN 
    1832                ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
    1833 !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 
    1834 !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 
    1835 !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 
    1836                !Alan is this right? 
    1837 !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 
    1838 !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 
    1839 !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) & 
    1840 !                       &      )/ (zdh(ji,jj)  + epsln ) 
    1841 !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 
    1842 !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 
    1843 !                     IF ( znd <= 0.0 ) THEN 
    1844 !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 
    1845 !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 
    1846 !                     ELSE 
    1847 !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 
    1848 !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 
    1849 !                     ENDIF 
    1850 !                  END DO 
    1851             ELSE 
    1852                ! stable conditions 
    1853                zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    1854                zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    1855                DO jk = 2, ibld(ji,jj) 
    1856                   znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    1857                   IF ( znd < 1.0 ) THEN 
    1858                      zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    1859                   ELSE 
    1860                      zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    1861                   ENDIF 
    1862                   zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
    1863                END DO 
    1864             ENDIF 
    1865             ! 
    1866          END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 
     1828            IF ( lconv(ji,jj) ) THEN   ! convective conditions 
     1829               IF ( lpyc(ji,jj) ) THEN 
     1830                  zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * zhol(ji,jj) ) ) ) 
     1831                  palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) *   & 
     1832                     &                                zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) /                & 
     1833                     &            ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) ) 
     1834                  palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) 
     1835                  ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 
     1836                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1837                  ! Commented lines in this section are not needed in new code, once tested ! 
     1838                  ! can be removed                                                          ! 
     1839                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1840                  ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 
     1841                  ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 
     1842                  zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 
     1843                  zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 
     1844                  DO jk = 2, ibld(ji,jj)+ibld_ext 
     1845                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 
     1846                     IF ( znd <= zzeta_m ) THEN 
     1847                        ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 
     1848                        ! &        EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1849                        ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 
     1850                        ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 
     1851                        pdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + palpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 
     1852                           & EXP( -6.0_wp * ( znd -zzeta_m )**2 ) 
     1853                     ELSE 
     1854                        ! zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1855                        ! zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 
     1856                        pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.0_wp * ppgamma_b * ( znd - zzeta_m )**2 ) 
     1857                     ENDIF 
     1858                  END DO 
     1859               ENDIF   ! If no pycnocline pycnocline gradients set to zero 
     1860            ELSE   ! Stable conditions 
     1861               ! If pycnocline profile only defined when depth steady of increasing. 
     1862               IF ( zdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady. 
     1863                  IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 
     1864                     IF ( zhol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline 
     1865                        ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 
     1866                        zbgrad = zdb_bl(ji,jj) * ztmp 
     1867                        DO jk = 2, ibld(ji,jj) 
     1868                           znd = gdepw(ji,jj,jk,Kmm) * ztmp 
     1869                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 
     1870                        END DO 
     1871                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     1872                        ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 
     1873                        zbgrad = zdb_bl(ji,jj) * ztmp 
     1874                        DO jk = 2, ibld(ji,jj) 
     1875                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 
     1876                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 
     1877                        END DO 
     1878                     ENDIF   ! IF (zhol >=0.5) 
     1879                  ENDIF      ! IF (zdb_bl> 0.) 
     1880               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 
     1881            ENDIF            ! IF (lconv) 
     1882         ENDIF   ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 
    18671883      END_2D 
    1868       IF( ln_timing ) CALL timing_stop('zdf_osm_pshp') 
    1869     END SUBROUTINE zdf_osm_pycnocline_shear_profiles 
     1884      ! 
     1885      IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 
     1886      ! 
     1887   END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 
    18701888 
    18711889   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 
     
    26092627     END IF 
    26102628 
     2629     ! Flags associated with diagnostic output 
     2630     IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) )                            ln_dia_pyc_shr = .TRUE. 
     2631     IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE. 
     2632      
    26112633     !                              ! allocate zdfosm arrays 
    26122634     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
Note: See TracChangeset for help on using the changeset viewer.