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 wn 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 zdfmxl ! mixed layer depth 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) #ifdef key_osm_debug INTEGER :: nn_idb = 297, nn_jdb = 193, nn_kdb = 35, nn_narea_db = 109 INTEGER :: iloc_db, jloc_db #endif ! ! 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 /rau0 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 !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $ !! 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 ) ! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ? ! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & ! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') ! IF ( ln_osm_mle ) THEN ! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') ! ENDIF IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) END FUNCTION zdf_osm_alloc SUBROUTINE zdf_osm( kt, 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 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) :: zwb0tot ! Total surface buoyancy flux including insolation 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) :: 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 ! Shear production. 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 ; zwb0tot(:,:) = 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 #ifdef key_osm_debug IF(narea==nn_narea_db)THEN iloc_db=mi0(nn_idb); jloc_db=mj0(nn_jdb) WRITE(narea+100,*) WRITE(narea+100,'(a,i7)')'timestep=',kt WRITE(narea+100,'(3(a,i7))')'narea=',narea,' nn_idb',nn_idb,' nn_jdb=',nn_jdb WRITE(narea+100,'(4(a,i7))')'iloc_db=',iloc_db,' jloc_db',jloc_db,' jpi=',jpi,' jpj=',jpj ji=iloc_db; jj=jloc_db WRITE(narea+100,'(a,i7,5(a,g10.2))')'mbkt=',mbkt(ji,jj),' ht_n',ht_n(ji,jj),& &' hu_n-',hu_n(ji-1,jj),' hu_n+',hu_n(ji,jj), ' hv_n-',hv_n(ji,jj-1),' hv_n+',hv_n(ji,jj) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif ! 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 jj = 2, jpjm1 DO ji = 2, jpim1 ! Surface downward irradiance (so always +ve) zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_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 DO END DO ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL DO jj = 2, jpjm1 DO ji = 2, jpim1 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_rau0_rcp * tmask(ji,jj,1) ! Upwards surface salinity flux for non-local term zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) ! Non radiative upwards surface buoyancy flux zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) ! Total upwards surface buoyancy flux zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * ( zrad0(ji,jj) - zradh(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_rau0 * tmask(ji,jj,1) zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rau0 * 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) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(4(3(a,g11.3),/), 2(a,g11.3),/)') & & 'after calculating fluxes: hbl=', hbl(ji,jj),' zthermal=',zthermal, ' zbeta=', zbeta,& & ' zrad0=', zrad0(ji,jj),' zradh=', zradh(ji,jj), ' zradav=', zradav(ji,jj), & & ' zwth0=', zwth0(ji,jj), ' zwthav=', zwthav(ji,jj), ' zws0=', zws0(ji,jj), & & ' zwb0=', zwb0(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwb0tot_in hbl=', zwb0tot(ji,jj) + grav * zthermal * zradh(ji,jj),& & ' zwbav=', zwbav(ji,jj) FLUSH(narea+100) END IF #endif END DO END DO ! 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO ! Assume Pierson-Moskovitz wind-wave spectrum CASE(1) DO jj = 2, jpjm1 DO ji = 2, jpim1 ! 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 DO END DO ! Use ECMWF wave fields as output from SBCWAVE CASE(2) zfac = 2.0_wp * rpi / 16.0_wp DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END SELECT #ifdef key_osm_debug IF(narea==nn_narea_db)THEN WRITE(narea+100,'(2(a,g11.3))') & & 'Before reduction: zustke=', zustke(iloc_db,jloc_db),' dstokes =',dstokes(iloc_db,jloc_db) FLUSH(narea+100) END IF #endif IF (ln_zdfosm_ice_shelter) THEN ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END SELECT ! Langmuir velocity scale (zwstrl), La # (zla) ! mixed scale (zvstr), convective velocity scale (zwstrc) DO jj = 2, jpjm1 DO ji = 2, jpim1 ! 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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,2(a,g11.3),/)') & & 'After reduction: zustke=', zustke(ji,jj), ' dstokes=', dstokes(ji,jj), & & ' zustar =', zustar(ji,jj), ' zwstrl=', zwstrl(ji,jj), ' zwstrc=', zwstrc(ji,jj),& & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj) FLUSH(narea+100) END IF #endif END DO END DO !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! 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_n(:,:,4) ) ibld(:,:) = 4 DO jk = 5, jpkm1 DO jj = 1, jpj DO ji = 1, jpi IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN ibld(ji,jj) = MIN(mbkt(ji,jj), jk) ENDIF END DO END DO END DO ! ########################################################################## DO jj = 2, jpjm1 DO ji = 2, jpim1 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) - 1 )) , 1 )) zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) END DO END DO #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,2(a,i7),/)') & & 'Before updating hbl: hbl=', hbl(ji,jj), ' dh=', dh(ji,jj), & &' zhbl =',zhbl(ji,jj) , ' zhml=', zhml(ji,jj), ' zdh=', zdh(ji,jj),& &' imld=', imld(ji,jj), ' ibld=', ibld(ji,jj) WRITE(narea+100,'(a,g11.3,a,2g11.3)') 'Physics: ssh ',sshn(ji,jj),' T S surface=',tsn(ji,jj,1,jp_tem),tsn(ji,jj,1,jp_sal) jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a,*(g11.3))') ' T[imld-1..ibld+2] =', ( tsn(ji,jj,jk,jp_tem), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' S[imld-1..ibld+2] =', ( tsn(ji,jj,jk,jp_sal), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' U+[imld-1..ibld+2] =', ( un(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' U-[imld-1..ibld+2] =', ( un(ji-1,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' V+[imld-1..ibld+2] =', ( vn(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' V-[imld-1..ibld+2] =', ( vn(ji,jj-1,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' W[imld-1..ibld+2] =', ( wn(ji,jj-1,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 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(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) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & & 'After averaging, with old hbl (& jp_ext==2), hml: zt_bl=', zt_bl(ji,jj),& & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),& & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),& & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) FLUSH(narea+100) END IF #endif ! 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 ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & & 'After rotation, with old hbl (& jp_ext==2), hml:', & & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) FLUSH(narea+100) END IF #endif ! 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 ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(2(a,l7),a, i7,/,3(a,g11.3),/)') & & 'After zdf_osm_osbl_state: lconv=', lconv(ji,jj), ' lshear=', lshear(ji,jj), ' j_ddh=', j_ddh(ji,jj),& & 'zwb_ent=', zwb_ent(ji,jj), ' zwb_min=', zwb_min(ji,jj), ' zshear=', zshear(ji,jj) FLUSH(narea+100) END IF #endif IF ( ln_osm_mle ) THEN ! Fox-Kemper Scheme mld_prof = 4 DO jk = 5, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END DO END DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) END DO END DO #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(2(a,g11.3), a, i7,/,2(3(a,g11.3),/),4(a,g11.3),/)') & & 'Before updating hmle: hmle =',hmle(ji,jj) , ' zhmle=', zhmle(ji,jj), ' mld_prof=', mld_prof(ji,jj), & & 'averaging over hmle: zt_mle=', zt_mle(ji,jj), ' zs_mle=', zs_mle(ji,jj), ' zb_mle=', zb_mle(ji,jj),& & 'zdt_mle=', zdt_mle(ji,jj), ' zds_mle=', zds_mle(ji,jj), ' zdb_mle=', zdb_mle(ji,jj),& & 'zu_mle =', zu_mle(ji,jj), ' zv_mle=', zv_mle(ji,jj), ' zdu_mle=', zdu_mle(ji,jj), ' zdv_mle=', zdv_mle(ji,jj) FLUSH(narea+100) END IF #endif !! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) !! Calculate vertical gradients immediately below zmld CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) !! calculate max vertical FK flux zwb_fk & set logical descriptors CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) !! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(a,g11.3,a,i7,/, 2(4(a,g11.3),/),2(a,g11.3),/,2(3(a,g11.3),/),a,i7,2(a,g11.3),/,3(a,g11.3),/,/)') & & 'Before updating hmle: zmld =',zmld(ji,jj),' mld_prof=', mld_prof(ji,jj), & & 'zdtdx+=', zdtdx(ji,jj),' zdtdx-=', zdtdx(ji-1,jj),' zdsdx+=', zdsdx(ji,jj),' zdsdx-=',zdsdx(ji-1,jj), & & 'zdtdy+=', zdtdy(ji,jj),' zdtdy-=', zdtdy(ji,jj-1),' zdsdy+=', zdsdy(ji,jj),' zdsdy-=',zdsdy(ji,jj-1), & & 'dbdx_mle+=', dbdx_mle(ji,jj),' dbdx_mle-=', dbdx_mle(ji-1,jj),& & 'dbdy_mle+=', dbdy_mle(ji,jj),' dbdy_mle-=',dbdy_mle(ji,jj-1),' zdbds_mle=',zdbds_mle(ji,jj), & & 'zdtdz_mle_ext=', zdtdz_mle_ext(ji,jj), ' zdsdz_mle_ext=', zdsdz_mle_ext(ji,jj), & & ' zdbdz_mle_ext=', zdbdz_mle_ext(ji,jj), & & 'After updating hmle: mld_prof=', mld_prof(ji,jj),' hmle=', hmle(ji,jj), ' zhmle=', zhmle(ji,jj),& & 'zvel_mle =', zvel_mle(ji,jj), ' zdiff_mle=', zdiff_mle(ji,jj), ' zwb_fk=', zwb_fk(ji,jj) FLUSH(narea+100) END IF #endif ELSE ! ln_osm_mle ! FK not selected, Boundary Layer only. lpyc(:,:) = .TRUE. lflux(:,:) = .FALSE. lmle(:,:) = .FALSE. DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. END DO END DO ENDIF ! ln_osm_mle !! External gradient below BL needed both with and w/o FK CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! Test if pycnocline well resolved DO jj = 2, jpjm1 DO ji = 2,jpim1 IF (lconv(ji,jj) ) THEN ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 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 DO END DO #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(4(a,l7),a,i7,/, 3(a,g11.3),/)') & & 'BL logical descriptors: lconv =',lconv(ji,jj),' lpyc=', lpyc(ji,jj),' lflux=', lflux(ji,jj),' lmle=', lmle(ji,jj),& & ' jp_ext=', jp_ext(ji,jj), & & 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj) FLUSH(narea+100) END IF #endif ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 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) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & & 'After averaging, with old hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),& & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),& & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),& & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) FLUSH(narea+100) END IF #endif ! Rate of change of hbl CALL zdf_osm_calculate_dhdt( zdhdt ) DO jj = 2, jpjm1 DO ji = 2, jpim1 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it ! adjustment to represent limiting by ocean bottom IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) lpyc(ji,jj) = .FALSE. ENDIF #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3))')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),& & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_rdt,' delta hbl from w ', wn(ji,jj,ibld(ji,jj))*rn_rdt FLUSH(narea+100) END IF #endif END DO END DO imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index ibld(:,:) = 4 DO jk = 4, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN ibld(ji,jj) = jk ENDIF END DO END DO END DO ! ! 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? ! Recalculate BL averages and differences using new BL depth 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO 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 ! Recalculate ML averages and differences using new ML depth 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 ! #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & & 'After averaging, with new hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),& & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),& & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),& & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) FLUSH(narea+100) END IF #endif 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 ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & & 'After rotation, with new hbl (& correct jp_ext), hml:', & & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) FLUSH(narea+100) END IF #endif !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! 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 ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a,l7,/,3(a,g11.3),/)') & & 'After pycnocline profiles BL lpyc=', lpyc(ji,jj),& & 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj), & & 'Pycnocline: zalpha_pyc=', zalpha_pyc(ji,jj) WRITE(narea+100,'(a,*(g11.3))') ' zdtdz_pyc[imld-1..ibld+2] =', ( zdtdz_pyc(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zdsdz_pyc[imld-1..ibld+2] =', ( zdsdz_pyc(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zdbdz_pyc[imld-1..ibld+2] =', ( zdbdz_pyc(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zdudz_pyc[imld-1..ibld+2] =', ( zdudz_pyc(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zdvdz_pyc[imld-1..ibld+2] =', ( zdvdz_pyc(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zviscos[imld-1..ibld+2] =', ( zviscos(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif ! ! 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 jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_d = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) / dstokes(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) END DO ENDIF ! endif for check on lconv END DO ! end of ji loop END DO ! end of jj loop ! 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 jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_d = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) / 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 DO ! ji loop END DO ! jj loo #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a,g11.3)')'Stokes contrib to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj) WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) IF( lconv(ji,jj) ) THEN WRITE(narea+100,'(3(a,g11.3))')'Stokes contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), & &' zsc_uw_2=',zsc_uw_2(ji,jj) ELSE WRITE(narea+100,'(2(a,g11.3))')'Stokes contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj) END IF WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif ! 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 ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) ELSEWHERE zsc_wth_1 = 0._wp zsc_ws_1 = 0._wp ENDWHERE DO jj = 2, jpjm1 DO ji = 2, jpim1 IF (lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) ! calculate turbulent time scale zl_c = 0.9 * ( 1.0 - EXP ( - 5.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & & * ( 1.0 - EXP ( -15.0 * ( 1.2 - zznd_ml ) ) ) zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & & * ( 1.0 - EXP ( - 8.0 * ( 1.15 - 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.4 * zsc_wth_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 * zsc_ws_1(ji,jj) * zl_eps / ( 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 DO jk = 2, ibld(ji,jj) zznd_pyc = -( gdepw_n(ji,jj,jk) - 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 ! IF ( dh(ji,jj) < 0.2*hbl(ji,jj) ) THEN 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_n(ji,jj,jk) - 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 END IF 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 DO ! ji loop END DO ! jj loop 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 jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2 , imld(ji,jj) zznd_d = gdepw_n(ji,jj,jk) / 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 DO ! ji loop END DO ! jj loop DO jj = 2, jpjm1 DO ji = 2, jpim1 IF( lconv(ji,jj) ) THEN 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 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 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_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse * & & ( 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_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse * & & ( 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 ENDIF ! lconv END DO ! ji loop END DO ! jj loop #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc contribs to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj) WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) IF( lconv(ji,jj) ) THEN WRITE(narea+100,'(3(a,g11.3))')'Stokes + buoy + pyc contribs to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), & &' zsc_uw_2=',zsc_uw_2(ji,jj) ELSE WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc contribs to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj) END IF WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif 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 jj = 2, jpjm1 DO ji = 2, jpim1 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.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) zsc_ws_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) ENDIF ELSE zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) zsc_ws_1(ji,jj) = zws0(ji,jj) ENDIF END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml=gdepw_n(ji,jj,jk) / 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 ! ! may need to comment out lpyc block IF ( lpyc(ji,jj) ) THEN ! pycnocline DO jk = imld(ji,jj), ibld(ji,jj) zznd_pyc = - ( gdepw_n(ji,jj,jk) - 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_n(ji,jj,jk) / dstokes(ji,jj) znd = gdepw_n(ji,jj,jk) / 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 ENDDO ! ji loop END DO ! jj loop 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 jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) zznd_d = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) / zhbl(ji,jj) zznd_d = gdepw_n(ji,jj,jk) / 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 DO END DO #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/contrib to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj) IF (lpyc(ji,jj)) WRITE(narea+100,'(2(a,g11.3))') 'zsc_wth_pyc=', zsc_wth_pyc(ji,jj), ' zsc_wth_pyc=',zsc_wth_pyc(ji,jj) WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) IF( lconv(ji,jj) ) THEN WRITE(narea+100,'(2(a,g11.3))')'Unstable; transport contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj) ELSE WRITE(narea+100,'(3(a,g11.3))')'Stable; transport contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), & &' zsc_uw_2=',zsc_uw_2(ji,jj) END IF WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif 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 jj = 2, jpjm1 DO ji = 2, jpim1 IF ( .not. lconv(ji,jj) ) THEN DO jk = 2, ibld(ji,jj) znd = -( gdepw_n(ji,jj,jk) - 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 DO END DO ! pynocline contributions DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( .not. lconv(ji,jj) ) THEN IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN DO jk= 2, ibld(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 DO END DO 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 jj=2, jpjm1 DO ji = 2, jpim1 ghamt(ji,jj,ibld(ji,jj)) = 0._wp ghams(ji,jj,ibld(ji,jj)) = 0._wp ghamu(ji,jj,ibld(ji,jj)) = 0._wp ghamv(ji,jj,ibld(ji,jj)) = 0._wp END DO ! ji loop END DO ! jj loop #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a)')'Tweak gham[uv] to go to zero near surface, add pycnocline viscosity/diffusivity & set=0 at ibld' WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO 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 jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) DO jj = 1, jpjm1 DO ji = 1, jpim1 ! vector opt. z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) END DO END DO END DO ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 ! vector opt. ! ! 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 DO END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END IF ! ln_kpprimix = .true. ! KPP-style set diffusivity large if unstable below BL IF( ln_convmix) THEN DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END IF ! ln_convmix = .true. #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a)') ' After including KPP Ri# diffusivity & viscosity' WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' zviscos[imld-1..ibld+2] =', ( zviscos(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing DO jj = 2 , jpjm1 DO ji = 2, jpim1 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_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(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_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(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_n(ji,jj,jk) / 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_n(ji,jj,jk) / 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 DO END DO #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a)') ' After including FK diffusivity & non-local terms' WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif 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( zviscos(:,:,:), 'W', 1. ) ! GN 25/8: need to change tmask --> wmask DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END DO ! 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. , p_avm, 'W', 1., & & ghamu, 'W', 1. , ghamv, 'W', 1. ) DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END DO ! 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 unchanged) CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & & ghamu, 'U', -1. , ghamv, 'V', -1. ) #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(a)') ' Final diffusivity & viscosity, & non-local terms' WRITE(narea+100,'(a,*(g11.3))') ' p_avt[imld-1..ibld+2] =', ( p_avt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' p_avm[imld-1..ibld+2] =', ( p_avm(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm ) WRITE(narea+100,*) FLUSH(narea+100) END IF #endif 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.*rau0*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.*rau0*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("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*zwb0 ) ! IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb buoyancy flux 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("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base 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.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*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("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolved internal to zdf_osm routine IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*j_ddh ) ! index forpyc thicknessh internal to zdf_osm routine IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKE 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("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 jj = 2, jpjm1 DO ji = 2, jpim1 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) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac FLUSH(narea+100) END IF #endif IF ( lpyc(ji,jj) ) THEN zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) FLUSH(narea+100) END IF #endif IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) ENDIF #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) FLUSH(narea+100) END IF #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) ) 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) ) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) FLUSH(narea+100) END IF #endif zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) FLUSH(narea+100) END IF #endif zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_n_sc(ji,jj) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) FLUSH(narea+100) END IF #endif 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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj) FLUSH(narea+100) END IF #endif ELSE ! conv, stable 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) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) FLUSH(narea+100) END IF #endif END IF END DO END DO ! DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) ! mixed layer diffusivity zznd_ml = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) - 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 - zb_cubic ) zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic DO jk = imld(ji,jj) , ibld(ji,jj) zznd_pyc = -( gdepw_n(ji,jj,jk) - 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_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 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_n(ji,jj,jk) / 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_n(ji, jj, ibld(ji,jj)) zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) ENDIF ENDIF ! end if ( lconv ) ! END DO ! end of ji loop END DO ! end of jj loop END SUBROUTINE zdf_osm_diffusivity_viscosity SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) !!--------------------------------------------------------------------- !! *** 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 ! Local Variables INTEGER :: jj, ji REAL(wp), DIMENSION(jpi,jpj) :: zekman REAL(wp), DIMENSION(jpi,jpj) :: 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.8 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 REAL, PARAMETER :: zri_c = 0.25 REAL, PARAMETER :: zek = 4.0 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. ! Determins stability and set flag lconv DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( zhol(ji,jj) < 0._wp ) THEN lconv(ji,jj) = .TRUE. ELSE lconv(ji,jj) = .FALSE. ENDIF END DO END DO zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) zshear(:,:) = 0._wp #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db WRITE(narea+100,'(a,g11.3)') & & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj) FLUSH(narea+100) END IF #endif j_ddh(:,:) = 1 DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN IF ( zdb_bl(ji,jj) > 0._wp ) THEN zri_p(ji,jj) = 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 ) IF ( ff_t(ji,jj) >= 0._wp ) THEN ! Northern Hemisphere zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )**2 + MAX( -zdv_ml(ji,jj), 1.e-5)**2 ) ELSE ! Southern Hemisphere zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )**2 + MAX( zdv_ml(ji,jj), 1.e-5)**2 ) ENDIF 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 ) ) #ifdef key_osm_debug ! IF(narea==nn_narea_db)THEN ! WRITE(narea+100,'(2(a,i10.4))')'ji',ji,'jj',jj ! WRITE(narea+100,'(2(a,i10.4))')'iloc_db',iloc_db,'jloc_db',jloc_db ! WRITE(narea+100,'(2(a,i10.4))')'iloc_db+',mi0(nn_idb),'jloc_db+',mj0(nn_jdb) ! FLUSH(narea+100) ! END IF IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',zshear(ji,jj) WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj) FLUSH(narea+100) END IF #endif ! Stability Dependence zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',zshear(ji,jj) FLUSH(narea+100) END IF #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! ! full code available ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF ( zshear(ji,jj) > 1.e-10 ) THEN IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN(hu_n(ji,jj), hu_n(ji-1,jj), hv_n(ji,jj), hv_n(ji,jj-1))>100._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 ENDIF ELSE j_ddh(ji,jj) = 2 lshear(ji,jj) = .FALSE. ENDIF ! 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 ELSE ! zdb_bl test, note zshear set to zero j_ddh(ji,jj) = 2 lshear(ji,jj) = .FALSE. ENDIF ENDIF END DO END DO ! Calculate entrainment buoyancy flux due to surface fluxes. DO jj = 2, jpjm1 DO ji = 2, jpim1 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 * zalpha_c * zrf_conv * zwbav(ji,jj) & & - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & & + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & & - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj) ! #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',zwb_ent(ji,jj) FLUSH(narea+100) END IF #endif ENDIF END DO ! ji loop END DO ! jj loop zwb_min(:,:) = 0._wp DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr FLUSH(narea+100) END IF #endif IF ( j_ddh(ji,jj) == 0 ) THEN ! ! Developing shear layer, additional shear production possible. ! zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp ) ! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) ! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) ! zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) ! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, & & ' zshear=',zshear(ji,jj),' zshear_u=', zshear_u FLUSH(narea+100) END IF #endif 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 ENDIF ! lshear IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2._wp * zwbav(ji,jj) ENDIF ! lconv #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), & & ' zwb_min=',zwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwbav= ', zwbav(ji,jj) FLUSH(narea+100) END IF #endif END DO ! ji END DO ! jj 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 jj = 2, jpjm1 ! Vertical slab DO ji = 2, jpim1 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_n(ji,jj,jk) zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & & / 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) - tsn(ji,jj,ibld_ext,jp_tem) zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal) zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) & & / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) & & / 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 DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 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) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) END DO END DO 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, zdbdz_mle_int znd_param(:,:) = 0._wp DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 ! 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 zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) DO jk = ibld(ji,jj), mld_prof(ji,jj) zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) ) zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(4(a,g11.3))')'start of zdf_osm_osbl_state_fk: zwb_fk=',zwb_fk(ji,jj), & & ' znd_param=',znd_param(ji,jj), ' zpe_mle_ref=', zpe_mle_ref, ' zpe_mle_layer=', zpe_mle_layer FLUSH(narea+100) END IF #endif END DO END DO ! Diagnosis DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lconv(ji,jj) ) THEN IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 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(ji,jj) = .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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',zwb_ent(ji,jj), & & ' zhmle=',zhmle(ji,jj), ' zhbl=', zhbl(ji,jj), & & ' lpyc= ', lpyc(ji,jj), ' lflux= ', lflux(ji,jj), ' lmle= ', lmle(ji,jj), ' lconv= ', lconv(ji,jj) FLUSH(narea+100) END IF #endif END DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 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) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & & / e3w_n(ji,jj,jkb1) zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & & / e3w_n(ji,jj,jkb1) 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 DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 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) znd = -( gdepw_n(ji,jj,jk) - 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 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,/,3(a,g11.3),/,2(a,g11.3),/)')'end of zdf_osm_pycnocline_scalar_profiles:lconv=lpyc=T',& & 'zzeta_m=', zzeta_m, ' zalpha=', zalpha(ji,jj), ' ztmp=', ztmp,& & ' zbgrad=', zbgrad, ' zgamma_b_nd=', zgamma_b_nd FLUSH(narea+100) END IF #endif 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_n(ji,jj,jk) * 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_n(ji,jj,jk) - 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) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_pycnocline_scalar_profiles:lconv=F ztgrad=',& & ztgrad, ' zsgrad=', zsgrad, ' zbgrad=', zbgrad FLUSH(narea+100) END IF #endif 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 DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 ! 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_n(ji,jj,jk) - 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_n(ji,jj,jk) / 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 DO END DO END SUBROUTINE zdf_osm_pycnocline_shear_profiles SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_calculate_dhdt *** !! !! ** Purpose : Calculates the rate at which hbl changes. !! !! ** Method : !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl INTEGER :: jj, ji REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi REAL(wp) :: zvel_max,zddhdt REAL(wp), PARAMETER :: zzeta_m = 0.3 REAL(wp), PARAMETER :: zgamma_c = 2.0 REAL(wp), PARAMETER :: zdhoh = 0.1 REAL(wp), PARAMETER :: zalpha_b = 0.3 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth DO jj = 2, jpjm1 DO ji = 2, jpim1 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 = ( zvstr(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) > 1.0e-15 ) THEN zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0._wp ) * zdh(ji,jj) / zdb_ml(ji,jj) zpsi = ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * ( zwb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp ) ) * zdh(ji,jj) / zhbl(ji,jj) zpsi = zpsi + 1.75 * ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) )*( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) * MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp ) zpsi = zalpha_b * MAX ( 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) ) + zpsi / ( zvel_max + MAX( zdb_bl(ji,jj), 1.e-15 ) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt, OSBL is deepening, entrainment > restratification: zdhdt=',zdhdt(ji,jj) WRITE(narea+100,'(3(a,g11.3))') ' zpsi=',zpsi, ' zgamma_b_nd=', zgamma_b_nd, ' zdh=', zdh(ji,jj) FLUSH(narea+100) END IF #endif 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 = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt,j_ddh(ji,jj) == 1: zari=',zari FLUSH(narea+100) END IF #endif ELSE IF ( j_ddh(ji,jj) == 0 ) THEN ! Growing shear layer zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt ELSE zddhdt = 0._wp ENDIF ! j_ddh zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( 1.0 -0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * & & zdb_ml(ji,jj) * zddhdt / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 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) = - MIN(zvel_mle(ji,jj), hbl(ji,jj)/10800.) ENDIF ELSE ! Fox-Kemper not used. zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / 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_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) ELSE zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) ENDIF zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) zdhdt(ji,jj) = MAX(zdhdt(ji,jj), -hbl(ji,jj)/5400.) 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) = - MIN(zvel_mle(ji,jj), hbl(ji,jj)/10800.) 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( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) ENDIF zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) zdhdt(ji,jj) = MAX(zdhdt(ji,jj), -hbl(ji,jj)/5400.) ENDIF ! lconv ENDIF ! lshear #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(4(a,g11.3))')'end of 1st major loop of zdf_osm_calculate_dhdt: zdhdt=',zdhdt(ji,jj), & & ' zpert=', zpert, ' zddhdt=', zddhdt, ' zvel_max=', zvel_max IF ( ln_osm_mle ) THEN WRITE(narea+100,'(3(a,g11.3),/)') 'zvel_mle=',zvel_mle(ji,jj), ' zwb_fk_b=', zwb_fk_b(ji,jj), & & ' zwb_ent + 2*zwb_fk_b =', zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) FLUSH(narea+100) END IF END IF #endif END DO END DO 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 jj = 2, jpjm1 DO ji = 2, jpim1 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',imld(ji,jj),' trial ibld=', ibld(ji,jj) FLUSH(narea+100) END IF #endif 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_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird ENDIF #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a,g11.3)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=T: zvel_max=',zvel_max FLUSH(narea+100) END IF #endif DO jk = imld(ji,jj), ibld(ji,jj) zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & & 0.0 ) + zvel_max IF ( ln_osm_mle ) THEN zhbl_s = zhbl_s + MIN( & & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w_n(ji,jj,jm) ) ELSE zhbl_s = zhbl_s + MIN( & & rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w_n(ji,jj,jm) ) ENDIF ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj) FLUSH(narea+100) END IF #endif END DO hbl(ji,jj) = zhbl_s ibld(ji,jj) = jm ELSE ! stable #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=F' FLUSH(narea+100) END IF #endif DO jk = imld(ji,jj), ibld(ji,jj) zdb = MAX( & & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& & 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_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' zhol',zhol(ji,jj),' zdhdt',zdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj) FLUSH(narea+100) END IF #endif END DO ENDIF ! IF ( lconv ) hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) ibld(ji,jj) = MAX(jm, 4 ) ELSE ! change zero or one model level. hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) ENDIF zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', zhbl(ji,jj),' ibld=', ibld(ji,jj) FLUSH(narea+100) END IF #endif END DO END DO 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, zddhdt, zvel_max REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN IF ( zdb_bl(ji,jj) > 1.0e-15) THEN IF ( j_ddh(ji,jj) == 0 ) THEN zvel_max = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) ! ddhdt for pycnocline determined in osm_calculate_dhdt zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) ) ELSE ! Need to recalculate because hbl has been updated. 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_rdt ) dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) ENDIF ELSE ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0 * rn_rdt ) dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(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_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / 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 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_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / 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_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / 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_n(ji,jj,ibld(ji,jj)-1), 1.e-3) ) , 1 ) imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(4(a,g11.3),2(a,i7),/,5(a,g11.3),/)') 'end of zdf_osm_pycnocline_thickness:hml=',hml(ji,jj), & & ' zhml=',zhml(ji,jj),' zdh=', zdh(ji,jj), ' dh=', dh(ji,jj), ' imld=', imld(ji,jj), ' inhml=', inhml, & & 'zvel_max=', zvel_max, ' ztau=', ztau,' zdh_ref=', zdh_ref, ' zar=', zari, ' zddhdt=', zddhdt FLUSH(narea+100) END IF #endif END DO END DO 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_rau0 ! convert density criteria into N^2 criteria DO jk = nlb10, jpkm1 DO jj = 1, jpj ! Mixed layer level: w-level DO ji = 1, jpi ikt = mbkt(ji,jj) zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) END DO END DO ! ensure mld_prof .ge. ibld ! ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation ! ztm(:,:) = 0._wp zsm(:,:) = 0._wp DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer DO jj = 1, jpj DO ji = 1, jpi zc = e3t_n(ji,jj,jk) * 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 * tsn(ji,jj,jk,jp_tem) zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) END DO END DO END DO ! average temperature and salinity. ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) ! calculate horizontal gradients at u & v points DO jj = 2, jpjm1 DO ji = 1, jpim1 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 DO END DO DO jj = 1, jpjm1 DO ji = 2, jpim1 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 DO END DO ! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on ztsm_midu(:,jpj,jp_sal) = 10. ztsm_midv(jpi,:,jp_sal) = 10. CALL eos_rab(ztsm_midu, zmld_midu, zabu) CALL eos_rab(ztsm_midv, zmld_midv, zabv) DO jj = 2, jpjm1 DO ji = 1, jpim1 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) END DO END DO DO jj = 1, jpjm1 DO ji = 2, jpim1 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO END SUBROUTINE zdf_osm_zmld_horizontal_gradients SUBROUTINE zdf_osm_mle_parameters( zmld, 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 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 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 jj = 2, jpjm1 DO ji = 2, jpim1 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 DO END DO ! Timestep mixed layer eddy depth. DO jj = 2, jpjm1 DO ji = 2, jpim1 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 * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) ) zdb_mle = zb_bl(ji,jj) - zbuoy ! Timestep hmle. hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_rdt / zdb_mle ELSE IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ELSE hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau ENDIF ENDIF hmle(ji,jj) = MAX(MIN(hmle(ji,jj), ht_n(ji,jj)), gdepw_n(ji,jj,4)) IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) ) ! For now try just set hmle to zmld hmle(ji,jj) = zmld(ji,jj) END DO END DO mld_prof = 4 DO jk = 5, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END DO END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) END DO END DO END SUBROUTINE zdf_osm_mle_parameters END SUBROUTINE zdf_osm SUBROUTINE zdf_osm_init !!---------------------------------------------------------------------- !! *** 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 :: ios ! local integer INTEGER :: ji, jj, jk ! dummy loop indices REAL z1_t2 !! #ifdef key_osm_debug 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 & & ,nn_idb, nn_jdb, nn_kdb, nn_narea_db #else 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 #endif ! 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 !!---------------------------------------------------------------------- ! REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 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 #ifdef key_osm_debug WRITE(numout,*) 'nn_idb', nn_idb, 'nn_jdb', nn_jdb, 'nn_kdb', nn_kdb, 'nn_narea_db', nn_narea_db iloc_db = mi0(nn_idb) jloc_db = mj0(nn_jdb) WRITE(numout,*) 'iloc_db ', iloc_db , 'jloc_db', jloc_db #endif 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 REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 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 /rau0 ! 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 jj=1,jpj do ji = 1,jpi r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) end do end do ! 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, '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 jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 ! vector opt. 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 DO END DO END DO 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 jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 ! vector opt. 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 DO END DO END DO 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. ! IF( lwxios ) THEN CALL iom_set_rstw_var_active('wn') CALL iom_set_rstw_var_active('hbl') CALL iom_set_rstw_var_active('dh') IF( ln_osm_mle ) THEN CALL iom_set_rstw_var_active('hmle') END IF ENDIF END SUBROUTINE zdf_osm_init SUBROUTINE osm_rst( kt, 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 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_autoglo, 'wn', wn, ldxios = lrxios ) WRITE(numout,*) ' ===>>>> : wn read from restart file' ELSE wn(:,:,:) = 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_autoglo, 'hbl' , hbl , ldxios = lrxios ) CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios ) 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_autoglo, 'hmle' , hmle , ldxios = lrxios ) 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 hbli into the restart file, then return IF(lwp) WRITE(numout,*) '---- osm-rst ----' CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn, ldxios = lwxios ) CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios ) CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) IF( ln_osm_mle ) THEN CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios ) 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( tsn, rab_n ) CALL bn2(tsn, rab_n, rn2) 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_rau0 ! convert density criteria into N^2 criteria ! hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 DO jk = 1, jpkm1 DO jj = 1, jpj ! Mixed layer level: w-level DO ji = 1, jpi ikt = mbkt(ji,jj) hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END DO END DO END DO ! DO jj = 1, jpj DO ji = 1, jpi iiki = MAX(4,imld_rst(ji,jj)) hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth END DO END DO WRITE(numout,*) ' ===>>>> : hbl computed from stratification' IF( ln_osm_mle ) THEN hmle(:,:) = hbl(:,:) ! Initialise MLE depth. WRITE(numout,*) ' ===>>>> : hmle set = to hbl' END IF wn(:,:,:) = 0._wp WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' END SUBROUTINE osm_rst SUBROUTINE tra_osm( kt ) !!---------------------------------------------------------------------- !! *** 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 INTEGER :: ji, jj, jk ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF IF( l_trdtra ) THEN !* Save ta and sa trends ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem) ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal) ENDIF DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & & - ( ghamt(ji,jj,jk ) & & - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk) tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & & - ( ghams(ji,jj,jk ) & & - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) END DO END DO END DO ! save the non-local tracer flux trends for diagnostics IF( l_trdtra ) THEN ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) ENDIF IF(ln_ctl) THEN CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, & & tab3d_2=tsa(:,:,:,jp_sal), 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 ) !!---------------------------------------------------------------------- !! *** 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 ! ! 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 jk = 1, jpkm1 ! add non-local u and v fluxes DO jj = 2, jpjm1 DO ji = 2, jpim1 ua(ji,jj,jk) = ua(ji,jj,jk) & & - ( ghamu(ji,jj,jk ) & & - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) va(ji,jj,jk) = va(ji,jj,jk) & & - ( ghamv(ji,jj,jk ) & & - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) END DO END DO END DO ! ! code for saving tracer trends removed ! END SUBROUTINE dyn_osm !!====================================================================== END MODULE zdfosm