MODULE zdfosm !!====================================================================== !! *** MODULE zdfosm *** !! Ocean physics: vertical mixing coefficient compute from the OSMOSIS !! turbulent closure parameterization !!===================================================================== !! History : NEMO 4.0 ! A. Grant, G. Nurser !! 15/03/2017 Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG !! 15/03/2017 Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G !! 06/06/2017 (1) Checks on sign of buoyancy jump in calculation of OSBL depth. A.G. !! (2) Removed variable zbrad0, zbradh and zbradav since they are not used. !! (3) Approximate treatment for shear turbulence. !! Minimum values for zustar and zustke. !! Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers. !! Limit maximum value for Langmuir number. !! Use zvstr in definition of stability parameter zhol. !! (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla !! (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative. !! (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop. !! (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems. !! (8) Change upper limits from ibld-1 to ibld. !! (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri. !! (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL. !! (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles. !! (12) Replace zwstrl with zvstr in calculation of eddy viscosity. !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information !! (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. !! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) !! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, !! (a) Pycnocline temperature and salinity profies changed for unstable layers !! (b) The stable OSBL depth parametrization changed. !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. !! 23/05/19 (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 'ln_zdfosm' OSMOSIS scheme !!---------------------------------------------------------------------- !! zdf_osm : update momentum and tracer Kz from osm scheme !! zdf_osm_init : initialization, namelist read, and parameters control !! osm_rst : read (or initialize) and write osmosis restart fields !! tra_osm : compute and add to the T & S trend the non-local flux !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) !! dyn_osm : compute and add to u & v trensd the non-local flux !! !! Subroutines in revised code. !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers ! uses ww from previous time step (which is now wb) to calculate hbl USE dom_oce ! ocean space and time domain USE zdf_oce ! ocean vertical physics USE sbc_oce ! surface boundary condition: ocean USE sbcwave ! surface wave parameters USE phycst ! physical constants USE eosbn2 ! equation of state USE traqsr ! details of solar radiation absorption USE zdfddm ! double diffusion mixing (avs array) USE iom ! I/O library USE lib_mpp ! MPP library USE trd_oce ! ocean trends definition USE trdtra ! tracers trends ! USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) IMPLICIT NONE PRIVATE PUBLIC zdf_osm ! routine called by step.F90 PUBLIC zdf_osm_init ! routine called by nemogcm.F90 PUBLIC osm_rst ! routine called by step.F90 PUBLIC tra_osm ! routine called by step.F90 PUBLIC trc_osm ! routine called by trcstp.F90 PUBLIC dyn_osm ! routine called by step.F90 PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamt !: non-local temperature flux (gamma/o) REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghams !: non-local salinity flux (gamma/o) REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift. REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer. ! !!** Namelist namzdf_osm ** LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation REAL(wp) :: rn_osm_la ! Turbulent Langmuir number REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift REAL(wp) :: rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by REAL(wp) :: rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl LOGICAL :: ln_zdfosm_ice_shelter ! flag to activate ice sheltering REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave INTEGER :: nn_osm_SD_reduce ! = 0/1/2 flag for getting effective stokes drift from surface value LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la LOGICAL :: ln_kpprimix = .true. ! Shear instability mixing REAL(wp) :: rn_riinfty = 0.7 ! local Richardson Number limit for shear instability REAL(wp) :: rn_difri = 0.005 ! maximum shear mixing at Rig = 0 (m2/s) LOGICAL :: ln_convmix = .true. ! Convective instability mixing REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s) ! OSMOSIS mixed layer eddy parametrization constants INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt REAL(wp) :: rn_osm_mle_ce ! MLE coefficient ! ! parameters used in nn_osm_mle = 0 case REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer ! ! parameters used in nn_osm_mle = 1 case REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front LOGICAL :: ln_osm_hmle_limit ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld REAL(wp) :: rn_osm_hmle_limit ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base. REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. ! !!! ** General constants ** REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 INTEGER :: idebug = 236 INTEGER :: jdebug = 228 !! * Substitutions # include "do_loop_substitute.h90" # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION zdf_osm_alloc() !!---------------------------------------------------------------------- !! *** FUNCTION zdf_osm_alloc *** !!---------------------------------------------------------------------- ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') END FUNCTION zdf_osm_alloc SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm *** !! !! ** Purpose : Compute the vertical eddy viscosity and diffusivity !! coefficients and non local mixing using the OSMOSIS scheme !! !! ** Method : The boundary layer depth hosm is diagnosed at tracer points !! from profiles of buoyancy, and shear, and the surface forcing. !! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from !! !! Kx = hosm Wx(sigma) G(sigma) !! !! and the non local term ghamt = Cs / Ws(sigma) / hosm !! Below hosm the coefficients are the sum of mixing due to internal waves !! shear instability and double diffusion. !! !! -1- Compute the now interior vertical mixing coefficients at all depths. !! -2- Diagnose the boundary layer depth. !! -3- Compute the now boundary layer vertical mixing coefficients. !! -4- Compute the now vertical eddy vicosity and diffusivity. !! -5- Smoothing !! !! N.B. The computation is done from jk=2 to jpkm1 !! Surface value of avt are set once a time to zero !! in routine zdf_osm_init. !! !! ** Action : update the non-local terms ghamts !! update avt (before vertical eddy coef.) !! !! References : Large W.G., Mc Williams J.C. and Doney S.C. !! Reviews of Geophysics, 32, 4, November 1994 !! Comments in the code refer to this paper, particularly !! the equation number. (LMD94, here after) !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time step INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) !! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: jl ! dummy loop indices INTEGER :: ikbot, jkmax, jkm1, jkp2 ! REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube ! REAL(wp) :: zbeta, zthermal ! REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density INTEGER :: jm ! dummy loop indices REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43 REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages ! Scales REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s) REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units) REAL(wp), DIMENSION(jpi,jpj) :: zradav ! Radiative flux, bl average (Buoyancy Units) REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number. REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux REAL(wp), DIMENSION(jpi,jpj) :: zvw0 ! Surface v-momentum flux REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic) REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux REAL(wp), DIMENSION(jpi,jpj) :: zwb_min REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer. LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness. ! mixed-layer variables INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency REAL(wp), DIMENSION(jpi,jpj) :: zddhdt ! correction to dhdt due to internal structure. REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer ! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. ! Flux-gradient relationship variables REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale. REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline. REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep ! For calculating Ri#-dependent mixing REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du ! u-shear^2 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv ! v-shear^2 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion ! Temporary variables INTEGER :: inhml REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables REAL(wp) :: zthick, zz0, zz1 ! temporary variables REAL(wp) :: zvel_max, zhbl_s ! temporary variables REAL(wp) :: zfac, ztmp ! temporary variable REAL(wp) :: zus_x, zus_y ! temporary Stokes drift REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau REAL(wp) :: zzeta_s = 0._wp REAL(wp) :: zzeta_v = 0.46 REAL(wp) :: zabsstke REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc ! For debugging INTEGER :: ikt !!-------------------------------------------------------------------- ! ibld(:,:) = 0 ; imld(:,:) = 0 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp zvw0(:,:) = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:) = 0._wp ; zwb0(:,:) = 0._wp zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp zhol(:,:) = 0._wp lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. ! mixed layer ! no initialization of zhbl or zhml (or zdh?) zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp zv_bl(:,:) = 0._wp ; zb_bl(:,:) = 0._wp zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp zt_mle(:,:) = 0._wp ; zs_mle(:,:) = 0._wp ; zu_mle(:,:) = 0._wp zb_mle(:,:) = 0._wp zv_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp zdb_ml(:,:) = 0._wp zdt_mle(:,:) = 0._wp ; zds_mle(:,:) = 0._wp ; zdu_mle(:,:) = 0._wp zdv_mle(:,:) = 0._wp ; zdb_mle(:,:) = 0._wp zwth_ent = 0._wp ; zws_ent = 0._wp ! zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp ! zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp IF ( ln_osm_mle ) THEN ! only initialise arrays if needed zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp ENDIF zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt ! Flux-Gradient arrays. zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp zhbl_t(:,:) = 0._wp ; zdhdt(:,:) = 0._wp zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp zddhdt(:,:) = 0._wp ! hbl = MAX(hbl,epsln) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Calculate boundary layer scales !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! Assume two-band radiation model for depth of OSBL zz0 = rn_abs ! surface equi-partition in 2-bands zz1 = 1. - rn_abs DO_2D( 0, 0, 0, 0 ) ! Surface downward irradiance (so always +ve) zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp ! Downwards irradiance at base of boundary layer zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) ! Downwards irradiance averaged over depth of the OSBL zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) END_2D ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL DO_2D( 0, 0, 0, 0 ) zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) ! Upwards surface Temperature flux for non-local term zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface salinity flux for non-local term 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) ! Non radiative upwards surface buoyancy flux zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) ! turbulent heat flux averaged over depth of OSBL zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) ! turbulent salinity flux averaged over depth of the OBSL zwsav(ji,jj) = 0.5 * zws0(ji,jj) ! turbulent buoyancy flux averaged over the depth of the OBSBL zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) ! Surface upward velocity fluxes zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) ! Friction velocity (zustar), at T-point : LMD94 eq. 2 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) END_2D ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) SELECT CASE (nn_osm_wave) ! Assume constant La#=0.3 CASE(0) DO_2D( 0, 0, 0, 0 ) zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 ! Linearly zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) dstokes(ji,jj) = rn_osm_dstokes END_2D ! Assume Pierson-Moskovitz wind-wave spectrum CASE(1) DO_2D( 0, 0, 0, 0 ) ! Use wind speed wndm included in sbc_oce module zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) END_2D ! Use ECMWF wave fields as output from SBCWAVE CASE(2) zfac = 2.0_wp * rpi / 16.0_wp DO_2D( 0, 0, 0, 0 ) IF (hsw(ji,jj) > 1.e-4) THEN ! Use wave fields zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) ELSE ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) ! .. so default to Pierson-Moskowitz zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) END IF END_2D END SELECT IF (ln_zdfosm_ice_shelter) THEN ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice DO_2D( 0, 0, 0, 0 ) zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) END_2D END IF SELECT CASE (nn_osm_SD_reduce) ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020). CASE(0) ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas. ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation. ! It could represent the effects of the spread of wave directions ! around the mean wind. The effect of this adjustment needs to be tested. IF(nn_osm_wave > 0) THEN zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) END IF CASE(1) ! van Roekel (2012): consider average SD over top 10% of boundary layer ! assumes approximate depth profile of SD from Breivik (2016) zsqrtpi = SQRT(rpi) z_two_thirds = 2.0_wp / 3.0_wp DO_2D( 0, 0, 0, 0 ) zthickness = rn_osm_hblfrac*hbl(ji,jj) z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) zsqrt_depth = SQRT(z2k_times_thickness) zexp_depth = EXP(-z2k_times_thickness) zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth & & - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & & + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness END_2D CASE(2) ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer ! assumes approximate depth profile of SD from Breivik (2016) zsqrtpi = SQRT(rpi) DO_2D( 0, 0, 0, 0 ) zthickness = rn_osm_hblfrac*hbl(ji,jj) z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) IF(z2k_times_thickness < 50._wp) THEN zsqrt_depth = SQRT(z2k_times_thickness) zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) ELSE ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness ! See Abramowitz and Stegun, Eq. 7.1.23 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp END IF zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) END_2D END SELECT ! Langmuir velocity scale (zwstrl), La # (zla) ! mixed scale (zvstr), convective velocity scale (zwstrc) DO_2D( 0, 0, 0, 0 ) ! Langmuir velocity scale (zwstrl), at T-point zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) ! Velocity scale that tends to zustar for large Langmuir numbers zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. ! Note zustke and zwstrl are not amended. ! ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv IF ( zwbav(ji,jj) > 0.0) THEN zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) ELSE zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) ENDIF END_2D !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! BL must be always 4 levels deep. ! For calculation of lateral buoyancy gradients for FK in ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must ! previously exist for hbl also. ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway ! ########################################################################## hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) ibld(:,:) = 4 DO_3D( 1, 1, 1, 1, 5, jpkm1 ) IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN ibld(ji,jj) = MIN(mbkt(ji,jj), jk) ENDIF END_3D ! ########################################################################## DO_2D( 0, 0, 0, 0 ) zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 )) zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) END_2D ! Averages over well-mixed and boundary layer jp_ext(:,:) = 2 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) ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) ! Velocity components in frame aligned with surface stress. CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) ! Determine the state of the OSBL, stable/unstable, shear/no shear CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) IF ( ln_osm_mle ) THEN ! Fox-Kemper Scheme mld_prof = 4 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END_3D jp_ext_mle(:,:) = 2 CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) DO_2D( 0, 0, 0, 0 ) zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) END_2D !! External gradient CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) ELSE ! ln_osm_mle ! FK not selected, Boundary Layer only. lpyc(:,:) = .TRUE. lflux(:,:) = .FALSE. lmle(:,:) = .FALSE. DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. END_2D ENDIF ! ln_osm_mle ! Test if pycnocline well resolved DO_2D( 0, 0, 0, 0 ) IF (lconv(ji,jj) ) THEN ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) IF ( ztmp > 6 ) THEN ! pycnocline well resolved jp_ext(ji,jj) = 1 ELSE ! pycnocline poorly resolved jp_ext(ji,jj) = 0 ENDIF ELSE ! Stable conditions jp_ext(ji,jj) = 0 ENDIF END_2D 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 ) ! jp_ext = ibld-imld+1 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) ! Rate of change of hbl CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) DO_2D( 0, 0, 0, 0 ) zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it ! adjustment to represent limiting by ocean bottom IF ( zhbl_t(ji,jj) >= gdepw(ji, jj, mbkt(ji,jj) + 1, Kmm ) ) THEN zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm) - depth_tol)! ht(:,:)) lpyc(ji,jj) = .FALSE. ENDIF END_2D imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index ibld(:,:) = 4 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN ibld(ji,jj) = jk ENDIF END_3D ! ! Step through model levels taking account of buoyancy change to determine the effect on dhdt ! CALL zdf_osm_timestep_hbl( zdhdt ) ! is external level in bounds? 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 ) ! ! ! Check to see if lpyc needs to be changed CALL zdf_osm_pycnocline_thickness( dh, zdh ) DO_2D( 0, 0, 0, 0 ) 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. END_2D dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. ! ! Average over the depth of the mixed layer in the convective boundary layer ! jp_ext = ibld - imld +1 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 ) ! rotate mean currents and changes onto wind align co-ordinates ! CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Pycnocline gradients for scalars and velocity !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) ! ! calculate non-gradient components of the flux-gradient relationships ! ! Stokes term in scalar flux, flux-gradient relationship WHERE ( lconv ) zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln) ! zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) ELSEWHERE zsc_wth_1 = 2.0 * zwthav ! zsc_ws_1 = 2.0 * zwsav ENDWHERE DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj) END DO ! end jk loop ELSE ! else for if (lconv) ! Stable conditions DO jk = 2, ibld(ji,jj) zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) END DO ENDIF ! endif for check on lconv END_2D ! 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) WHERE ( lconv ) 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 ) zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 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 ) ELSEWHERE zsc_uw_1 = zustar**2 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) ENDWHERE IF(ln_dia_osm) THEN IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) END IF DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) ! ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) & & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) END DO ! end jk loop ELSE ! Stable conditions DO jk = 2, ibld(ji,jj) ! corrected to ibld zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp END DO ! end jk loop ENDIF END_2D ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] WHERE ( lconv ) zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) ELSEWHERE zsc_wth_1 = 0._wp zsc_ws_1 = 0._wp ENDWHERE DO_2D( 0, 0, 0, 0 ) IF (lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) ! calculate turbulent length scale zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) ! non-gradient buoyancy terms ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) END DO IF ( lpyc(ji,jj) ) THEN ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) zwth_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) zws_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) ! Cubic profile used for buoyancy term za_cubic = 0.755 * ztau_sc_u(ji,jj) zb_cubic = 0.25 * ztau_sc_u(ji,jj) DO jk = 2, ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 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 ) END DO ! zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 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 ) ) ! zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) ! zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) ! zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) DO jk = 2, ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird ! 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 END DO ENDIF ! End of pycnocline ELSE ! lconv test - stable conditions DO jk = 2, ibld(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) END DO ENDIF END_2D WHERE ( lconv ) zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) zsc_uw_2 = zwb0 * zustke * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0) zsc_vw_1 = 0._wp ELSEWHERE zsc_uw_1 = 0._wp zsc_vw_1 = 0._wp ENDWHERE DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2 , imld(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & & * zsc_uw_2(ji,jj) ) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) END DO ! jk loop ELSE ! stable conditions DO jk = 2, ibld(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) END DO ENDIF END_2D DO_2D( 0, 0, 0, 0 ) IF ( lpyc(ji,jj) ) THEN IF ( j_ddh(ji,jj) == 0 ) THEN ! Place holding code. Parametrization needs checking for these conditions. zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) ELSE zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) ENDIF zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse zc_cubic = zuw_bse - zd_cubic ! need ztau_sc_u to be available. Change to array. DO jk = imld(ji,jj), ibld(ji,jj) zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) END DO zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse zc_cubic = zvw_bse - zd_cubic DO jk = imld(ji,jj), ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) END DO ENDIF ! lpyc END_2D IF(ln_dia_osm) THEN IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) END IF ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] DO_2D( 1, 0, 1, 0 ) IF ( lconv(ji,jj) ) THEN zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) IF ( lpyc(ji,jj) ) THEN ! Pycnocline scales zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) ENDIF ELSE zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) zsc_ws_1(ji,jj) = zws0(ji,jj) ENDIF END_2D DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & & - EXP( - 6.0 * zznd_ml ) ) ) & & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & & - EXP( - 6.0 * zznd_ml ) ) ) & & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) END DO ! IF ( lpyc(ji,jj) ) THEN ! pycnocline DO jk = imld(ji,jj), ibld(ji,jj) zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) END DO ENDIF ELSE IF( zdhdt(ji,jj) > 0. ) THEN DO jk = 2, ibld(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) END DO ENDIF ENDIF END_2D WHERE ( lconv ) zsc_uw_1 = zustar**2 zsc_vw_1 = ff_t * zustke * zhml ELSEWHERE zsc_uw_1 = zustar**2 zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1 zsc_vw_1 = ff_t * zustke * zhbl zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1 ENDWHERE DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) ! ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) END DO ELSE DO jk = 2, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) IF ( zznd_d <= 2.0 ) THEN ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) ! ELSE ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) ! ENDIF ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) END DO ENDIF END_2D IF(ln_dia_osm) THEN IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) END IF ! ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. DO_2D( 0, 0, 0, 0 ) IF ( .not. lconv(ji,jj) ) THEN DO jk = 2, ibld(ji,jj) znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about IF ( znd >= 0.0 ) THEN ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) ELSE ghamu(ji,jj,jk) = 0._wp ghamv(ji,jj,jk) = 0._wp ENDIF END DO ENDIF END_2D ! pynocline contributions DO_2D( 0, 0, 0, 0 ) IF ( .not. lconv(ji,jj) ) THEN IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN DO jk= 2, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) END DO END IF END IF END_2D IF(ln_dia_osm) THEN IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) END IF DO_2D( 0, 0, 0, 0 ) ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp END_2D IF(ln_dia_osm) THEN IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) END IF !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Need to put in code for contributions that are applied explicitly to ! the prognostic variables ! 1. Entrainment flux ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! rotate non-gradient velocity terms back to model reference frame DO_2D( 0, 0, 0, 0 ) DO jk = 2, ibld(ji,jj) ztemp = ghamu(ji,jj,jk) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) END DO END_2D IF(ln_dia_osm) THEN IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) END IF ! KPP-style Ri# mixing IF( ln_kpprimix) THEN DO_3D( 1, 0, 1, 0, 2, jpkm1 ) !* Shear production at uw- and vw-points (energy conserving form) z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) & & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) & & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) END_3D ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! ! shear prod. at w-point weightened by mask zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) ! ! local Richardson number zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) zfri = MIN( zri / rn_riinfty , 1.0_wp ) zfri = ( 1.0_wp - zfri * zfri ) zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) END_3D DO_2D( 0, 0, 0, 0 ) DO jk = ibld(ji,jj) + 1, jpkm1 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri END DO END_2D END IF ! ln_kpprimix = .true. ! KPP-style set diffusivity large if unstable below BL IF( ln_convmix) THEN DO_2D( 0, 0, 0, 0 ) DO jk = ibld(ji,jj) + 1, jpkm1 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv END DO END_2D END IF ! ln_convmix = .true. IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing DO_2D( 0, 0, 0, 0 ) IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer ! Calculate MLE flux contribution from surface fluxes DO jk = 1, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) END DO DO jk = 1, mld_prof(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) END DO ! Viscosity for MLEs DO jk = 1, mld_prof(ji,jj) znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) END DO ELSE ! Surface transports limited to OSBL. ! Viscosity for MLEs DO jk = 1, mld_prof(ji,jj) znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) END DO ENDIF END_2D ENDIF IF(ln_dia_osm) THEN IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) END IF ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids !CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) ! GN 25/8: need to change tmask --> wmask DO_3D( 0, 0, 0, 0, 2, jpkm1 ) p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) END_3D ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) END_3D ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & & ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) IF(ln_dia_osm) THEN SELECT CASE (nn_osm_wave) ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). CASE(0:1) IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind ) ! x surface Stokes drift IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind ) ! y surface Stokes drift IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) ! Stokes drift read in from sbcwave (=2). CASE(2:3) IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & & SQRT(ut0sd**2 + vt0sd**2 ) ) END SELECT IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt END IF CONTAINS ! subroutine code changed, needs syntax checking. SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_diffusivity_viscosity *** !! !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. !! !! ** Method : !! !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:) :: zdiffut REAL(wp), DIMENSION(:,:,:) :: zviscos ! local ! Scales used to calculate eddy diffusivity and viscosity profiles REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc ! REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) IF ( lpyc(ji,jj) ) THEN zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) ENDIF 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) ) zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 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) zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) ENDIF 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) ) ) zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 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 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) ELSE zbeta_d_sc(ji,jj) = 1.0 zbeta_v_sc(ji,jj) = 1.0 ENDIF ELSE zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) END IF END_2D ! DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) ! mixed layer diffusivity zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) ! zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 ! zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & & * ( 1.0 - 0.5 * zznd_ml**2 ) END DO ! pycnocline IF ( lpyc(ji,jj) ) THEN ! Diffusivity profile in the pycnocline given by cubic polynomial. za_cubic = 0.5 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 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 ) & & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic DO jk = imld(ji,jj) , ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) ! 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 ) 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 ) END DO ! viscosity profiles. za_cubic = 0.5 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 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) zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic DO jk = imld(ji,jj) , ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 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 ) 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 ) END DO IF ( zdhdt(ji,jj) > 0._wp ) THEN zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) ELSE zdiffut(ji,jj,ibld(ji,jj)) = 0._wp zviscos(ji,jj,ibld(ji,jj)) = 0._wp ENDIF ENDIF ELSE ! stable conditions DO jk = 2, ibld(ji,jj) zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) END DO IF ( zdhdt(ji,jj) > 0._wp ) THEN zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) ENDIF ENDIF ! end if ( lconv ) ! END_2D END SUBROUTINE zdf_osm_diffusivity_viscosity SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_osbl_state *** !! !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number !! !! ** Method : !! !! !!---------------------------------------------------------------------- INTEGER, DIMENSION(jpi,jpj) :: j_ddh ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low. LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number ! Local Variables INTEGER :: jj, ji REAL(wp), DIMENSION(jpi,jpj) :: zekman REAL(wp) :: zri_p, zri_b ! Richardson numbers REAL(wp) :: zshear_u, zshear_v, zwb_shr REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. ! Determins stability and set flag lconv DO_2D( 0, 0, 0, 0 ) IF ( zhol(ji,jj) < 0._wp ) THEN lconv(ji,jj) = .TRUE. ELSE lconv(ji,jj) = .FALSE. ENDIF END_2D zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) WHERE ( lconv ) zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) END WHERE zshear(:,:) = 0._wp j_ddh(:,:) = 1 DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN IF ( zdb_bl(ji,jj) > 0._wp ) THEN zri_p = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! ! full code available ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN ! Growing shear layer j_ddh(ji,jj) = 0 lshear(ji,jj) = .TRUE. ELSE j_ddh(ji,jj) = 1 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. lshear(ji,jj) = .TRUE. ELSE ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. zshear(ji,jj) = 0.5 * zshear(ji,jj) lshear(ji,jj) = .FALSE. ENDIF ENDIF ELSE ! zdb_bl test, note zshear set to zero j_ddh(ji,jj) = 2 lshear(ji,jj) = .FALSE. ENDIF ENDIF END_2D ! Calculate entrainment buoyancy flux due to surface fluxes. DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) IF (nn_osm_SD_reduce > 0 ) THEN ! Effective Stokes drift already reduced from surface value zr_stokes = 1.0_wp ELSE ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, ! requires further reduction where BL is deep zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) END IF zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) ! ENDIF END_2D zwb_min(:,:) = 0._wp DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_shr = -za_wb_s * zshear(ji,jj) IF ( j_ddh(ji,jj) == 0 ) THEN ! Developing shear layer, additional shear production possible. zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) zwb_shr = -za_wb_s * zshear(ji,jj) ENDIF zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) ELSE ! IF ( lconv ) THEN - ENDIF ! Stable OSBL - shear production not coded for first attempt. ENDIF ! lconv ELSE ! lshear IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_shr = -za_wb_s * zshear(ji,jj) zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) ENDIF ! lconv ENDIF ! lshear END_2D END SUBROUTINE zdf_osm_osbl_state SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_vertical_average *** !! !! ** Purpose : Determines vertical averages from surface to jnlev. !! !! ** Method : Averages are calculated from the surface to jnlev. !! The external level used to calculate differences is ibld+ibld_ext !! !!---------------------------------------------------------------------- INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. INTEGER, DIMENSION(jpi,jpj) :: jp_ext ! Alan: do we need zb? REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. INTEGER :: jk, ji, jj, ibld_ext REAL(wp) :: zthick, zthermal, zbeta zt = 0._wp zs = 0._wp zu = 0._wp zv = 0._wp DO_2D( 0, 0, 0, 0 ) zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? zbeta = rab_n(ji,jj,1,jp_sal) ! average over depth of boundary layer zthick = epsln DO jk = 2, jnlev_av(ji,jj) zthick = zthick + e3t(ji,jj,jk,Kmm) zt(ji,jj) = zt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) zs(ji,jj) = zs(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) zu(ji,jj) = zu(ji,jj) + e3t(ji,jj,jk,Kmm) & & * ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) & & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) zv(ji,jj) = zv(ji,jj) + e3t(ji,jj,jk,Kmm) & & * ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) & & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) END DO zt(ji,jj) = zt(ji,jj) / zthick zs(ji,jj) = zs(ji,jj) / zthick zu(ji,jj) = zu(ji,jj) / zthick zv(ji,jj) = zv(ji,jj) / zthick zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) IF ( ibld_ext < mbkt(ji,jj) ) THEN zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm) zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm) zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) & & / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) & & / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) ELSE zdt(ji,jj) = 0._wp zds(ji,jj) = 0._wp zdu(ji,jj) = 0._wp zdv(ji,jj) = 0._wp zdb(ji,jj) = 0._wp ENDIF END_2D END SUBROUTINE zdf_osm_vertical_average SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_velocity_rotation *** !! !! ** Purpose : Rotates frame of reference of averaged velocity components. !! !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline INTEGER :: ji, jj REAL(wp) :: ztemp DO_2D( 0, 0, 0, 0 ) ztemp = zu(ji,jj) zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) ztemp = zdu(ji,jj) zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) END_2D END SUBROUTINE zdf_osm_velocity_rotation SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_osbl_state_fk *** !! !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme. !! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined !! lflux :: determines whether effects of surface flux extend below the base of the OSBL !! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl. !! !! ** Method : !! !! !!---------------------------------------------------------------------- ! Outputs LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! REAL(wp), DIMENSION(jpi,jpj) :: znd_param REAL(wp) :: zbuoy, ztmp, zpe_mle_layer REAL(wp) :: zpe_mle_ref, zwb_ent, zdbdz_mle_int znd_param(:,:) = 0._wp DO_2D( 0, 0, 0, 0 ) ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) END_2D DO_2D( 0, 0, 0, 0 ) ! IF ( lconv(ji,jj) ) THEN IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) ! Calculate potential energies of actual profile and reference profile. zpe_mle_layer = 0._wp zpe_mle_ref = 0._wp DO jk = ibld(ji,jj), mld_prof(ji,jj) zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) END DO ! Non-dimensional parameter to diagnose the presence of thermocline znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) ENDIF ENDIF END_2D ! Diagnosis DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & & - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & & + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & & - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN ! MLE layer growing IF ( znd_param (ji,jj) > 100. ) THEN ! Thermocline present lflux(ji,jj) = .FALSE. lmle(ji,jj) =.FALSE. ELSE ! Thermocline not present lflux(ji,jj) = .TRUE. lmle(ji,jj) = .TRUE. ENDIF ! znd_param > 100 ! IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN lpyc(ji,jj) = .FALSE. ELSE lpyc = .TRUE. ENDIF ELSE ! MLE layer restricted to OSBL or just below. IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN ! Weak stratification MLE layer can grow. lpyc(ji,jj) = .FALSE. lflux(ji,jj) = .TRUE. lmle(ji,jj) = .TRUE. ELSE ! Strong stratification lpyc(ji,jj) = .TRUE. lflux(ji,jj) = .FALSE. lmle(ji,jj) = .FALSE. ENDIF ! zdb_bl < rn_mle_thresh_bl and ENDIF ! zhmle > 1.2 zhbl ELSE lpyc(ji,jj) = .TRUE. lflux(ji,jj) = .FALSE. lmle(ji,jj) = .FALSE. IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ELSE ! Stable Boundary Layer lpyc(ji,jj) = .FALSE. lflux(ji,jj) = .FALSE. lmle(ji,jj) = .FALSE. ENDIF ! lconv END_2D END SUBROUTINE zdf_osm_osbl_state_fk SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_external_gradients *** !! !! ** Purpose : Calculates the gradients below the OSBL !! !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. !! !!---------------------------------------------------------------------- INTEGER, DIMENSION(jpi,jpj) :: jbase REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. INTEGER :: jj, ji, jkb, jkb1 REAL(wp) :: zthermal, zbeta DO_2D( 0, 0, 0, 0 ) IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? zbeta = rab_n(ji,jj,1,jp_sal) jkb = jbase(ji,jj) jkb1 = MIN(jkb + 1, mbkt(ji,jj)) zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & & / e3t(ji,jj,ibld(ji,jj),Kmm) zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & & / e3t(ji,jj,ibld(ji,jj),Kmm) zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) ELSE zdtdz(ji,jj) = 0._wp zdsdz(ji,jj) = 0._wp zdbdz(ji,jj) = 0._wp END IF END_2D END SUBROUTINE zdf_osm_external_gradients SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline REAL(wp), DIMENSION(jpi,jpj) :: zalpha INTEGER :: jk, jj, ji REAL(wp) :: ztgrad, zsgrad, zbgrad REAL(wp) :: zgamma_b_nd, znd REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 DO_2D( 0, 0, 0, 0 ) IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! convective conditions IF ( lpyc(ji,jj) ) THEN zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 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 ) ) zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) ztmp = 1._wp/MAX(zdh(ji,jj), epsln) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Commented lines in this section are not needed in new code, once tested ! ! can be removed ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) DO jk = 2, ibld(ji,jj)+ibld_ext znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp IF ( znd <= zzeta_m ) THEN ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & & EXP( -6.0 * ( znd -zzeta_m )**2 ) ELSE ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) ENDIF END DO ENDIF ! if no pycnocline pycnocline gradients set to zero ELSE ! stable conditions ! if pycnocline profile only defined when depth steady of increasing. IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady. IF ( zdb_bl(ji,jj) > 0._wp ) THEN IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) ztgrad = zdt_bl(ji,jj) * ztmp zsgrad = zds_bl(ji,jj) * ztmp zbgrad = zdb_bl(ji,jj) * ztmp DO jk = 2, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) * ztmp zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) END DO ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. ztmp = 1._wp/MAX(zdh(ji,jj), epsln) ztgrad = zdt_bl(ji,jj) * ztmp zsgrad = zds_bl(ji,jj) * ztmp zbgrad = zdb_bl(ji,jj) * ztmp DO jk = 2, ibld(ji,jj) znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) END DO ENDIF ! IF (zhol >=0.5) ENDIF ! IF (zdb_bl> 0.) ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero ENDIF ! IF (lconv) ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) END_2D END SUBROUTINE zdf_osm_pycnocline_scalar_profiles SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** !! !! ** Purpose : Calculates velocity shear in the pycnocline !! !! ** Method : !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz INTEGER :: jk, jj, ji REAL(wp) :: zugrad, zvgrad, znd REAL(wp) :: zzeta_v = 0.45 ! DO_2D( 0, 0, 0, 0 ) ! IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN IF ( lconv (ji,jj) ) THEN ! Unstable conditions. Shouldn;t be needed with no pycnocline code. ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) !Alan is this right? ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & ! & )/ (zdh(ji,jj) + epsln ) ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v ! IF ( znd <= 0.0 ) THEN ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) ! ELSE ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) ! ENDIF ! END DO ELSE ! stable conditions zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) DO jk = 2, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) IF ( znd < 1.0 ) THEN zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) ELSE zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) ENDIF zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) END DO ENDIF ! END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) END_2D END SUBROUTINE zdf_osm_pycnocline_shear_profiles SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_calculate_dhdt *** !! !! ** Purpose : Calculates the rate at which hbl changes. !! !! ** Method : !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt ! Rate of change of hbl INTEGER :: jj, ji REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi REAL(wp) :: zvel_max!, zwb_min REAL(wp) :: zzeta_m = 0.3 REAL(wp) :: zgamma_c = 2.0 REAL(wp) :: zdhoh = 0.1 REAL(wp) :: alpha_bc = 0.5 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! Convective IF ( ln_osm_mle ) THEN IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN ! Fox-Kemper buoyancy flux average over OSBL zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) ELSE zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) ENDIF zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN ! OSBL is deepening, entrainment > restratification IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN ! *** Used for shear Needs to be changed to work stabily ! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml ! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) ! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) ! za_1 = 1.0 / zgamma_b**2 - 0.017 ! za_2 = 1.0 / zgamma_b**3 - 0.0025 ! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) zpsi = 0._wp zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) IF ( j_ddh(ji,jj) == 1 ) THEN IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ! Relaxation to dh_ref = zari * hbl zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) ELSE ! j_ddh == 0 ! Growing shear layer zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) ENDIF ! j_ddh zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) ELSE ! zdb_bl >0 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) ENDIF ELSE ! zwb_min + 2*zwb_fk_b < 0 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) zdhdt(ji,jj) = - zvel_mle(ji,jj) ENDIF ELSE ! Fox-Kemper not used. zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) ! added ajgn 23 July as temporay fix ENDIF ! ln_osm_mle ELSE ! lconv - Stable zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) IF ( zdhdt(ji,jj) < 0._wp ) THEN ! For long timsteps factor in brackets slows the rapid collapse of the OSBL zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) ELSE zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) ENDIF zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) ENDIF ! lconv ELSE ! lshear IF ( lconv(ji,jj) ) THEN ! Convective IF ( ln_osm_mle ) THEN IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN ! Fox-Kemper buoyancy flux average over OSBL zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) ELSE zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) ENDIF zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN ! OSBL is deepening, entrainment > restratification IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) ELSE zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) ENDIF ELSE ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) zdhdt(ji,jj) = - zvel_mle(ji,jj) ENDIF ELSE ! Fox-Kemper not used. zvel_max = -zwb_ent(ji,jj) / & & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) ! added ajgn 23 July as temporay fix ENDIF ! ln_osm_mle ELSE ! Stable zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) IF ( zdhdt(ji,jj) < 0._wp ) THEN ! For long timsteps factor in brackets slows the rapid collapse of the OSBL zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) ELSE zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) ENDIF zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) ENDIF ! lconv ENDIF ! lshear END_2D END SUBROUTINE zdf_osm_calculate_dhdt SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_timestep_hbl *** !! !! ** Purpose : Increments hbl. !! !! ** Method : If thechange in hbl exceeds one model level the change is !! is calculated by moving down the grid, changing the buoyancy !! jump. This is to ensure that the change in hbl does not !! overshoot a stable layer. !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl. INTEGER :: jk, jj, ji, jm REAL(wp) :: zhbl_s, zvel_max, zdb REAL(wp) :: zthermal, zbeta DO_2D( 0, 0, 0, 0 ) IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN ! ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. ! zhbl_s = hbl(ji,jj) jm = imld(ji,jj) zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) IF ( lconv(ji,jj) ) THEN !unstable IF( ln_osm_mle ) THEN zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) ELSE zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird ENDIF DO jk = imld(ji,jj), ibld(ji,jj) zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & & 0.0 ) + zvel_max IF ( ln_osm_mle ) THEN zhbl_s = zhbl_s + MIN( & & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w(ji,jj,jm,Kmm) ) ELSE zhbl_s = zhbl_s + MIN( & & rn_Dt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w(ji,jj,jm,Kmm) ) ENDIF ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 END DO hbl(ji,jj) = zhbl_s ibld(ji,jj) = jm ELSE ! stable DO jk = imld(ji,jj), ibld(ji,jj) zdb = MAX( & & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )& & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& & 0.0 ) + & & 2.0 * zvstr(ji,jj)**2 / zhbl_s ! Alan is thuis right? I have simply changed hbli to hbl zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 END DO ENDIF ! IF ( lconv ) hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4,Kmm) ) ibld(ji,jj) = MAX(jm, 4 ) ELSE ! change zero or one model level. hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) ENDIF zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) END_2D END SUBROUTINE zdf_osm_timestep_hbl SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_pycnocline_thickness *** !! !! ** Purpose : Calculates thickness of the pycnocline !! !! ** Method : The thickness is calculated from a prognostic equation !! that relaxes the pycnocine thickness to a diagnostic !! value. The time change is calculated assuming the !! thickness relaxes exponentially. This is done to deal !! with large timesteps. !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. ! INTEGER :: jj, ji INTEGER :: inhml REAL(wp) :: zari, ztau, zdh_ref REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN IF ( j_ddh(ji,jj) == 0 ) THEN ! ddhdt for pycnocline determined in osm_calculate_dhdt dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt ELSE ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt ) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) ENDIF ELSE ! lconv ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here ! boundary layer deepening IF ( zdb_bl(ji,jj) > 0._wp ) THEN ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) ELSE zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! IF(dhdt < 0) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ! IF (dhdt >= 0) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse ! Alan: this hml is never defined or used -- do we need it? ENDIF ELSE ! lshear ! for lshear = .FALSE. calculate ddhdt here IF ( lconv(ji,jj) ) THEN IF( ln_osm_mle ) THEN IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE ! unstable zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = zari * hbl(ji,jj) ELSE ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! ln_osm_mle IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE ! unstable zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = zari * hbl(ji,jj) ELSE ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF END IF ! ln_osm_mle dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! Alan: this hml is never defined or used ELSE ! IF (lconv) ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here ! boundary layer deepening IF ( zdb_bl(ji,jj) > 0._wp ) THEN ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) ELSE zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! IF(dhdt < 0) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ! IF (dhdt >= 0) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse ENDIF ! IF (lconv) ENDIF ! lshear hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj),Kmm), 1.e-3) ) , 1 ) imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) END_2D END SUBROUTINE zdf_osm_pycnocline_thickness SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm_horizontal_gradients *** !! !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. !! !! ** Method : !! !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ii, ij, ik, ikmax ! local integers REAL(wp) :: zc REAL(wp) :: zN2_c ! local buoyancy difference from 10m value REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv !!---------------------------------------------------------------------- ! ! !== MLD used for MLE ==! mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) ikt = mbkt(ji,jj) zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END_3D DO_2D( 1, 1, 1, 1 ) mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) END_2D ! ensure mld_prof .ge. ibld ! ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation ! ztm(:,:) = 0._wp zsm(:,:) = 0._wp DO_3D( 1, 1, 1, 1, 1, ikmax ) zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) END_3D ! average temperature and salinity. ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) ! calculate horizontal gradients at u & v points DO_2D( 0, 0, 1, 0 ) zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) END_2D DO_2D( 1, 0, 0, 0 ) zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) END_2D CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm) CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm) DO_2D( 0, 0, 1, 0 ) dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) END_2D DO_2D( 1, 0, 0, 0 ) dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) END_2D DO_2D( 0, 0, 0, 0 ) ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) END_2D END SUBROUTINE zdf_osm_zmld_horizontal_gradients SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm_mle_parameters *** !! !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. !! !! ** Method : !! !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 INTEGER, DIMENSION(jpi,jpj) :: mld_prof REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers INTEGER , DIMENSION(jpi,jpj) :: inml_mle REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 ENDIF END_2D ! Timestep mixed layer eddy depth. DO_2D( 0, 0, 0, 0 ) IF ( lmle(ji,jj) ) THEN ! MLE layer growing. ! Buoyancy gradient at base of MLE layer. zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) jkb = mld_prof(ji,jj) jkb1 = MIN(jkb + 1, mbkt(ji,jj)) ! zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) ) zdb_mle = zb_bl(ji,jj) - zbuoy ! Timestep hmle. hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / zdb_mle ELSE IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau ELSE hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau ENDIF ENDIF hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) END_2D mld_prof = 4 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END_3D DO_2D( 0, 0, 0, 0 ) zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm) END_2D END SUBROUTINE zdf_osm_mle_parameters END SUBROUTINE zdf_osm SUBROUTINE zdf_osm_init( Kmm ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm_init *** !! !! ** Purpose : Initialization of the vertical eddy diffivity and !! viscosity when using a osm turbulent closure scheme !! !! ** Method : Read the namosm namelist and check the parameters !! called at the first timestep (nit000) !! !! ** input : Namlist namosm !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: Kmm ! time level INTEGER :: ios ! local integer INTEGER :: ji, jj, jk ! dummy loop indices REAL z1_t2 !! NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter ! Namelist for Fox-Kemper parametrization. NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit !!---------------------------------------------------------------------- ! READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) IF(lwm) WRITE ( numond, namzdf_osm ) IF(lwp) THEN ! Control print WRITE(numout,*) WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave SELECT CASE (nn_osm_wave) CASE(0) WRITE(numout,*) ' calculated assuming constant La#=0.3' CASE(1) WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves' CASE(2) WRITE(numout,*) ' calculated from ECMWF wave fields' END SELECT WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce WRITE(numout,*) ' fraction of hbl to average SD over/fit' WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac SELECT CASE (nn_osm_SD_reduce) CASE(0) WRITE(numout,*) ' No reduction' CASE(1) WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL' CASE(2) WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL' END SELECT WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv ENDIF ! ! Check wave coupling settings ! ! ! Further work needed - see ticket #2447 ! IF( nn_osm_wave == 2 ) THEN IF (.NOT. ( ln_wave .AND. ln_sdw )) & & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) END IF ! ! allocate zdfosm arrays IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) IF( ln_osm_mle ) THEN ! Initialise Fox-Kemper parametrization READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') IF(lwm) WRITE ( numond, namosm_mle ) IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' WRITE(numout,*) '~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namosm_mle : ' WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit ENDIF ! ENDIF ! IF(lwp) THEN WRITE(numout,*) IF( ln_osm_mle ) THEN WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' ELSE WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' ENDIF ENDIF ! IF( ln_osm_mle ) THEN ! MLE initialisation ! rb_c = grav * rn_osm_mle_rho_c /rho0 ! Mixed Layer buoyancy criteria IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' ! IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! ! ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) ! ENDIF ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) z1_t2 = 2.e-5 DO_2D( 1, 1, 1, 1 ) r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) END_2D ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) ! ENDIF call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle IF( ln_zdfddm) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' Double diffusion mixing on temperature and salinity ' WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm ' ENDIF ENDIF !set constants not in namelist !----------------------------- IF(lwp) THEN WRITE(numout,*) ENDIF IF (nn_osm_wave == 0) THEN dstokes(:,:) = rn_osm_dstokes END IF ! Horizontal average : initialization of weighting arrays ! ------------------- SELECT CASE ( nn_ave ) CASE ( 0 ) ! no horizontal average IF(lwp) WRITE(numout,*) ' no horizontal average on avt' IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !' ! weighting mean arrays etmean ! ( 1 1 ) ! avt = 1/4 ( 1 1 ) ! etmean(:,:,:) = 0.e0 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) etmean(ji,jj,jk) = tmask(ji,jj,jk) & & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) END_3D CASE ( 1 ) ! horizontal average IF(lwp) WRITE(numout,*) ' horizontal average on avt' ! weighting mean arrays etmean ! ( 1/2 1 1/2 ) ! avt = 1/8 ( 1 2 1 ) ! ( 1/2 1 1/2 ) etmean(:,:,:) = 0.e0 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) etmean(ji,jj,jk) = tmask(ji, jj,jk) & & / MAX( 1., 2.* tmask(ji,jj,jk) & & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) END_3D CASE DEFAULT WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave CALL ctl_stop( ctmp1 ) END SELECT ! Initialization of vertical eddy coef. to the background value ! ------------------------------------------------------------- DO jk = 1, jpk avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) END DO ! zero the surface flux for non local term and osm mixed layer depth ! ------------------------------------------------------------------ ghamt(:,:,:) = 0. ghams(:,:,:) = 0. ghamu(:,:,:) = 0. ghamv(:,:,:) = 0. ! END SUBROUTINE zdf_osm_init SUBROUTINE osm_rst( kt, Kmm, cdrw ) !!--------------------------------------------------------------------- !! *** ROUTINE osm_rst *** !! !! ** Purpose : Read or write BL fields in restart file !! !! ** Method : use of IOM library. If the restart does not contain !! required fields, they are recomputed from stratification !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time step index INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag INTEGER :: id1, id2, id3 ! iom enquiry index INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: iiki, ikt ! local integer REAL(wp) :: zhbf ! tempory scalars REAL(wp) :: zN2_c ! local scalar REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) !!---------------------------------------------------------------------- ! !!----------------------------------------------------------------------------- ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return !!----------------------------------------------------------------------------- IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) IF( id1 > 0 ) THEN ! 'wn' exists; read CALL iom_get( numror, jpdom_auto, 'wn', ww ) WRITE(numout,*) ' ===>>>> : wn read from restart file' ELSE ww(:,:,:) = 0._wp WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' END IF id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return CALL iom_get( numror, jpdom_auto, 'hbl' , hbl ) CALL iom_get( numror, jpdom_auto, 'dh', dh ) WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' IF( ln_osm_mle ) THEN id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. ) IF( id3 > 0) THEN CALL iom_get( numror, jpdom_auto, 'hmle' , hmle ) WRITE(numout,*) ' ===>>>> : hmle read from restart file' ELSE WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' hmle(:,:) = hbl(:,:) ! Initialise MLE depth. END IF END IF RETURN ELSE ! 'hbl' & 'dh' not in restart file, recalculate WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' END IF END IF !!----------------------------------------------------------------------------- ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return !!----------------------------------------------------------------------------- IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbl into the restart file, then return IF(lwp) WRITE(numout,*) '---- osm-rst ----' CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww ) CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl ) CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh ) IF( ln_osm_mle ) THEN CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle ) END IF RETURN END IF !!----------------------------------------------------------------------------- ! Getting hbl, no restart file with hbl, so calculate from surface stratification !!----------------------------------------------------------------------------- IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' ! w-level of the mixing and mixed layers CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm) imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria ! hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ikt = mbkt(ji,jj) hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END_3D ! DO_2D( 1, 1, 1, 1 ) iiki = MAX(4,imld_rst(ji,jj)) hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth END_2D WRITE(numout,*) ' ===>>>> : hbl computed from stratification' IF( ln_osm_mle ) THEN hmle(:,:) = hbl(:,:) ! Initialise MLE depth. WRITE(numout,*) ' ===>>>> : hmle set = to hbl' END IF ww(:,:,:) = 0._wp WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' END SUBROUTINE osm_rst SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_osm *** !! !! ** Purpose : compute and add to the tracer trend the non-local tracer flux !! !! ** Method : ??? !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! time step index INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation ! INTEGER :: ji, jj, jk ! IF( kt == nit000 ) THEN IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF ENDIF IF( l_trdtra ) THEN !* Save ta and sa trends ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) ENDIF DO_3D( 0, 0, 0, 0, 1, jpkm1 ) pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & & - ( ghamt(ji,jj,jk ) & & - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & & - ( ghams(ji,jj,jk ) & & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) END_3D ! save the non-local tracer flux trends for diagnostics IF( l_trdtra ) THEN ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) ENDIF IF(sn_cfctl%l_prtctl) THEN CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ENDIF ! END SUBROUTINE tra_osm SUBROUTINE trc_osm( kt ) ! Dummy routine !!---------------------------------------------------------------------- !! *** ROUTINE trc_osm *** !! !! ** Purpose : compute and add to the passive tracer trend the non-local !! passive tracer flux !! !! !! ** Method : ??? !!---------------------------------------------------------------------- ! !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt WRITE(*,*) 'trc_osm: Not written yet', kt END SUBROUTINE trc_osm SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_osm *** !! !! ** Purpose : compute and add to the velocity trend the non-local flux !! copied/modified from tra_osm !! !! ** Method : ??? !!---------------------------------------------------------------------- INTEGER , INTENT( in ) :: kt ! ocean time step index INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation ! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF !code saving tracer trends removed, replace with trdmxl_oce DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & & - ( ghamu(ji,jj,jk ) & & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & & - ( ghamv(ji,jj,jk ) & & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) END_3D ! ! code for saving tracer trends removed ! END SUBROUTINE dyn_osm !!====================================================================== END MODULE zdfosm