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_vertical_average : compute vertical averages over boundary layers !! zdf_osm_init : initialization, namelist read, and parameters control !! osm_rst : read (or initialize) and write osmosis restart fields !! tra_osm : compute and add to the T & S trend the non-local flux !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) !! dyn_osm : compute and add to u & v trensd the non-local flux !! !! Subroutines in revised code. !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers ! uses ww from previous time step (which is now wb) to calculate hbl USE dom_oce ! ocean space and time domain USE zdf_oce ! ocean vertical physics USE sbc_oce ! surface boundary condition: ocean USE sbcwave ! surface wave parameters USE phycst ! physical constants USE eosbn2 ! equation of state USE traqsr ! details of solar radiation absorption USE zdfdrg, ONLY : rCdU_bot ! bottom friction velocity USE zdfddm ! double diffusion mixing (avs array) USE iom ! I/O library USE lib_mpp ! MPP library USE trd_oce ! ocean trends definition USE trdtra ! tracers trends ! USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) USE timing, ONLY : timing_start, timing_stop ! Timer 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_dia_pyc_scl = .FALSE. ! Output of pycnocline scalar-gradient profiles LOGICAL :: ln_dia_pyc_shr = .FALSE. ! Output of pycnocline velocity-shear profiles 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 /rho0 where rho_c is defined in zdfmld REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base. REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. ! !!! ** General constants ** REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 INTEGER :: idebug = 236 INTEGER :: jdebug = 228 !! * Substitutions # include "do_loop_substitute.h90" # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION zdf_osm_alloc() !!---------------------------------------------------------------------- !! *** FUNCTION zdf_osm_alloc *** !!---------------------------------------------------------------------- ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') END FUNCTION zdf_osm_alloc SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm *** !! !! ** Purpose : Compute the vertical eddy viscosity and diffusivity !! coefficients and non local mixing using the OSMOSIS scheme !! !! ** Method : The boundary layer depth hosm is diagnosed at tracer points !! from profiles of buoyancy, and shear, and the surface forcing. !! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from !! !! Kx = hosm Wx(sigma) G(sigma) !! !! and the non local term ghamt = Cs / Ws(sigma) / hosm !! Below hosm the coefficients are the sum of mixing due to internal waves !! shear instability and double diffusion. !! !! -1- Compute the now interior vertical mixing coefficients at all depths. !! -2- Diagnose the boundary layer depth. !! -3- Compute the now boundary layer vertical mixing coefficients. !! -4- Compute the now vertical eddy vicosity and diffusivity. !! -5- Smoothing !! !! N.B. The computation is done from jk=2 to jpkm1 !! Surface value of avt are set once a time to zero !! in routine zdf_osm_init. !! !! ** Action : update the non-local terms ghamts !! update avt (before vertical eddy coef.) !! !! References : Large W.G., Mc Williams J.C. and Doney S.C. !! Reviews of Geophysics, 32, 4, November 1994 !! Comments in the code refer to this paper, particularly !! the equation number. (LMD94, here after) !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time step INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) !! INTEGER :: ji, jj, jk, jkflt ! dummy loop indices INTEGER :: jl ! dummy loop indices INTEGER :: ikbot, 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) :: 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) :: 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) :: lcoup ! Coupling to bottom 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 ! offset for external level INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 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) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in 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), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep ! For calculating Ri#-dependent mixing REAL(wp), DIMENSION(jpi,jpj) :: z2du ! u-shear^2 REAL(wp), DIMENSION(jpi,jpj) :: z2dv ! v-shear^2 REAL(wp) :: 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 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 REAL(wp) :: zlarge = -1.0e10_wp, zero = 0.0_wp !!-------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('zdf_osm') ibld(:,:) = 0 ; imld(:,:) = 0 zrad0(:,:) = zlarge ; zradh(:,:) = zlarge ; zustar(:,:) = zlarge zwstrl(:,:) = zlarge ; zvstr(:,:) = zlarge ; zwstrc(:,:) = zlarge ; zuw0(:,:) = zlarge zwth0(:,:) = zlarge ; zws0(:,:) = zlarge ; zwb0(:,:) = zlarge zwthav(:,:) = zlarge ; zwsav(:,:) = zlarge ; zwbav(:,:) = zlarge ; zwb_ent(:,:) = zlarge zustke(:,:) = zlarge ; zla(:,:) = zlarge ; zcos_wind(:,:) = zlarge ; zsin_wind(:,:) = zlarge zhol(:,:) = zlarge ; zwb0tot(:,:) = zlarge ; zalpha_pyc(:,:) = zlarge lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. ! mixed layer ! no initialization of zhbl or zhml (or zdh?) zhbl(:,:) = zlarge ; zhml(:,:) = zlarge ; zdh(:,:) = zlarge ; zdhdt(:,:) = zlarge zt_bl(:,:) = zlarge ; zs_bl(:,:) = zlarge ; zu_bl(:,:) = zlarge zv_bl(:,:) = zlarge ; zb_bl(:,:) = zlarge zt_ml(:,:) = zlarge ; zs_ml(:,:) = zlarge ; zu_ml(:,:) = zlarge zt_mle(:,:) = zlarge ; zs_mle(:,:) = zlarge ; zu_mle(:,:) = zlarge zb_mle(:,:) = zlarge zv_ml(:,:) = zlarge ; zdt_bl(:,:) = zlarge ; zds_bl(:,:) = zlarge zdu_bl(:,:) = zlarge ; zdv_bl(:,:) = zlarge ; zdb_bl(:,:) = zlarge zdt_ml(:,:) = zlarge ; zds_ml(:,:) = zlarge ; zdu_ml(:,:) = zlarge ; zdv_ml(:,:) = zlarge zdb_ml(:,:) = zlarge ! zdbdz_pyc(:,:,:) = zlarge zdbdz_pyc(A2D(0),:) = 0.0_wp ! zdtdz_bl_ext(:,:) = zlarge ; zdsdz_bl_ext(:,:) = zlarge ; zdbdz_bl_ext(:,:) = zlarge IF ( ln_osm_mle ) THEN ! only initialise arrays if needed zdtdx(:,:) = zlarge ; zdtdy(:,:) = zlarge ; zdsdx(:,:) = zlarge zdsdy(:,:) = zlarge ; dbdx_mle(:,:) = zlarge ; dbdy_mle(:,:) = zlarge zwb_fk(:,:) = zlarge ; zvel_mle(:,:) = zlarge ; zdiff_mle(:,:) = zlarge zhmle(:,:) = zlarge ; zmld(:,:) = zlarge ENDIF zwb_fk_b(:,:) = zlarge ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt zhbl_t(:,:) = zlarge zdiffut(:,:,:) = zlarge ; zviscos(:,:,:) = zlarge zdiffut(A2D(0),:) = 0.0_wp ; zviscos(A2D(0),:) = 0.0_wp ghamt(:,:,:) = zlarge ; ghams(:,:,:) = zlarge ghamt(A2D(0),:) = 0.0_wp ; ghams(A2D(0),:) = 0.0_wp ghamu(:,:,:) = zlarge ; ghamv(:,:,:) = zlarge ghamu(A2D(0),:) = 0.0_wp ; ghamv(A2D(0),:) = 0.0_wp zdiff_mle(A2D(0)) = 0.0_wp #ifdef key_osm_debug IF(mi0(nn_idb)==mi1(nn_idb) .AND. mj0(nn_jdb)==mj1(nn_jdb) .AND. & & mi0(nn_idb) > 1 .AND. mi0(nn_idb) < jpi .AND. mj0(nn_jdb) > 1 .AND. mj0(nn_jdb) < jpj) THEN nn_narea_db = narea 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(ji,jj),& &' hu_n-',hu(ji-1,jj,Kmm),' hu_n+',hu(ji,jj,Kmm), ' hv_n-',hv(ji,jj-1,Kmm),' hv_n+',hv(ji,jj,Kmm) WRITE(narea+100,*) FLUSH(narea+100) ELSE nn_narea_db = -1000 END IF #endif ! hbl = MAX(hbl,epsln) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Calculate boundary layer scales !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL zz0 = rn_abs ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands zz1 = 1.0_wp - rn_abs DO_2D( 0, 0, 0, 0 ) zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp ! Surface downward irradiance (so always +ve) zradh(ji,jj) = zrad0(ji,jj) * & ! Downwards irradiance at base of boundary layer & ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) ) zradav = zrad0(ji,jj) * & ! Downwards irradiance averaged over depth of the OSBL & ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + & & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface Temperature flux for non-local term zwthav(ji,jj) = 0.5_wp * zwth0(ji,jj) - & ! Turbulent heat flux averaged over depth of OSBL & ( 0.5_wp * ( zrad0(ji,jj) + zradh(ji,jj) ) - zradav ) END_2D DO_2D( 0, 0, 0, 0 ) zws0(ji,jj) = -1.0_wp * & ! Upwards surface salinity flux for non-local term & ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - & ! Non radiative upwards surface buoyancy flux & grav * zbeta * zws0(ji,jj) zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * & ! Total upwards surface buoyancy flux & ( zrad0(ji,jj) - zradh(ji,jj) ) zwsav(ji,jj) = 0.5 * zws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - & ! Turbulent buoyancy flux averaged over the depth of the & grav * zbeta * zwsav(ji,jj) ! OBSBL END_2D DO_2D( 0, 0, 0, 0 ) zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * & ! Surface upward velocity fluxes & r1_rho0 * tmask(ji,jj,1) zvw0 = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * & ! Friction velocity (zustar), at T-point : LMD94 eq. 2 & zuw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp ) zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) zsin_wind(ji,jj) = -zvw0 / ( zustar(ji,jj) * zustar(ji,jj) ) #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) zradav = zrad0(ji,jj) * ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + & & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) 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, & & ' 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_2D ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) SELECT CASE (nn_osm_wave) ! Assume constant La#=0.3 CASE(0) DO_2D( 0, 0, 0, 0 ) zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 ! Linearly zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) dstokes(ji,jj) = rn_osm_dstokes END_2D ! Assume Pierson-Moskovitz wind-wave spectrum CASE(1) DO_2D( 0, 0, 0, 0 ) ! Use wind speed wndm included in sbc_oce module zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) END_2D ! Use ECMWF wave fields as output from SBCWAVE CASE(2) zfac = 2.0_wp * rpi / 16.0_wp DO_2D( 0, 0, 0, 0 ) IF (hsw(ji,jj) > 1.e-4) THEN ! Use wave fields zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) ELSE ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) ! .. so default to Pierson-Moskowitz zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) END IF END_2D END SELECT #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_2D( 0, 0, 0, 0 ) zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) END_2D END IF SELECT CASE (nn_osm_SD_reduce) ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020). CASE(0) ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas. ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation. ! It could represent the effects of the spread of wave directions ! around the mean wind. The effect of this adjustment needs to be tested. IF(nn_osm_wave > 0) THEN zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) END IF CASE(1) ! van Roekel (2012): consider average SD over top 10% of boundary layer ! assumes approximate depth profile of SD from Breivik (2016) zsqrtpi = SQRT(rpi) z_two_thirds = 2.0_wp / 3.0_wp DO_2D( 0, 0, 0, 0 ) zthickness = rn_osm_hblfrac*hbl(ji,jj) z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) zsqrt_depth = SQRT(z2k_times_thickness) zexp_depth = EXP(-z2k_times_thickness) zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth & & - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & & + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness END_2D CASE(2) ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer ! assumes approximate depth profile of SD from Breivik (2016) zsqrtpi = SQRT(rpi) DO_2D( 0, 0, 0, 0 ) zthickness = rn_osm_hblfrac*hbl(ji,jj) z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) IF(z2k_times_thickness < 50._wp) THEN zsqrt_depth = SQRT(z2k_times_thickness) zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) ELSE ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness ! See Abramowitz and Stegun, Eq. 7.1.23 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp END IF zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) END_2D END SELECT ! Langmuir velocity scale (zwstrl), La # (zla) ! mixed scale (zvstr), convective velocity scale (zwstrc) DO_2D( 0, 0, 0, 0 ) ! Langmuir velocity scale (zwstrl), at T-point zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) ! Velocity scale that tends to zustar for large Langmuir numbers zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. ! Note zustke and zwstrl are not amended. ! ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv IF ( zwbav(ji,jj) > 0.0) THEN zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) ELSE zwstrc(ji,jj) = 0.0_wp 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),/,3(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), ' zvstr=', zvstr(ji,jj) FLUSH(narea+100) END IF #endif END_2D !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! BL must be always 4 levels deep. ! For calculation of lateral buoyancy gradients for FK in ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must ! previously exist for hbl also. ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway ! ########################################################################## hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) ibld(:,:) = 4 DO_3D( 1, 1, 1, 1, 5, jpkm1 ) IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN ibld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) ENDIF END_3D ! ########################################################################## DO_2D( 0, 0, 0, 0 ) zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj) - 1, Kmm )) , 1 )) zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) END_2D #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 ',ssh(ji,jj,Kmm),' T S surface=',ts(ji,jj,1,jp_tem,Kmm),ts(ji,jj,1,jp_sal,Kmm) 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] =', ( ts(ji,jj,jk,jp_tem,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' S[imld-1..ibld+2] =', ( ts(ji,jj,jk,jp_sal,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' U+[imld-1..ibld+2] =', ( uu(ji,jj,jk,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' U-[imld-1..ibld+2] =', ( uu(ji-1,jj,jk,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' V+[imld-1..ibld+2] =', ( vv(ji,jj,jk,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' V-[imld-1..ibld+2] =', ( vv(ji,jj-1,jk,Kmm), jk=jl,jm ) WRITE(narea+100,'(a,*(g11.3))') ' W[imld-1..ibld+2] =', ( ww(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(:,:) = 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, jp_ext, & & 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_3D( 0, 0, 0, 0, 5, jpkm1 ) IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END_3D CALL zdf_osm_vertical_average( Kbb, Kmm, & & mld_prof, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle ) DO_2D( 0, 0, 0, 0 ) zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) END_2D #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,/,(3(a,g11.3),/),2(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),& & 'zu_mle =', zu_mle(ji,jj), ' zv_mle=', zv_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_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. END_2D ENDIF ! ln_osm_mle !! External gradient below BL needed both with and w/o FK CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03 ! Test if pycnocline well resolved ! DO_2D( 0, 0, 0, 0 ) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity ! IF (lconv(ji,jj) ) THEN should account for this. ! ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) ! IF ( ztmp > 6 ) THEN ! ! pycnocline well resolved ! jp_ext(ji,jj) = 1 ! ELSE ! ! pycnocline poorly resolved ! jp_ext(ji,jj) = 0 ! ENDIF ! ELSE ! ! Stable conditions ! jp_ext(ji,jj) = 0 ! ENDIF ! END_2D #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.. jp_ext(:,:) = 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, & & jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) ! ag 19/03 #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 ) ! Test if surface boundary layer coupled to bottom lcoup(:,:) = .FALSE. ! ag 19/03 DO_2D( 0, 0, 0, 0 ) zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it ! adjustment to represent limiting by ocean bottom IF ( mbkt(ji,jj) > 2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access IF ( zhbl_t(ji,jj) > gdepw(ji, jj,mbkt(ji,jj)-2,Kmm) ) THEN zhbl_t(ji,jj) = MIN( zhbl_t(ji,jj), gdepw(ji,jj,mbkt(ji,jj)-2,Kmm) ) ! ht(:,:)) lpyc(ji,jj) = .FALSE. lcoup(ji,jj) = .TRUE. ! ag 19/03 END IF END IF #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)),2(a,l7)')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),& & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_Dt,' delta hbl from w ', ww(ji,jj,ibld(ji,jj))*rn_Dt, & & ' lcoup= ', lcoup(ji,jj), ' lpyc= ', lpyc(ji,jj) FLUSH(narea+100) END IF #endif END_2D imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index ibld(:,:) = 4 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN ibld(ji,jj) = jk ENDIF END_3D ! ! Step through model levels taking account of buoyancy change to determine the effect on dhdt ! CALL zdf_osm_timestep_hbl( zdhdt ) ! is external level in bounds? ! Recalculate BL averages and differences using new BL depth jp_ext(:,:) = 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) CALL zdf_osm_pycnocline_thickness( dh, zdh ) ! Reset l_pyc before calculating terms in the flux-gradient relationship DO_2D( 0, 0, 0, 0 ) IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) >= mbkt(ji,jj) - 2 .or. & & ibld(ji,jj)-imld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0.0_wp ) THEN ! ag 19/03 lpyc(ji,jj) = .FALSE. ! ag 19/03 IF ( ibld(ji,jj) >= mbkt(ji,jj) -2 ) THEN imld(ji,jj) = ibld(ji,jj) - 1 ! ag 19/03 zdh(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) - gdepw(ji,jj,imld(ji,jj),Kmm) ! ag 19/03 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) ! ag 19/03 dh(ji,jj) = zdh(ji,jj) ! ag 19/03 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) ! ag 19/03 #ifdef key_osm_debug IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= F5: ibld(ji,jj) >= mbkt(ji,jj) -2' WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj) WRITE(narea+100,'(2(a,g11.3))')'dh=',dh(ji,jj),' hml=',hml(ji,jj) FLUSH(narea+100) END IF #endif ENDIF ENDIF ! ag 19/03 END_2D dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. ! ! Average over the depth of the mixed layer in the convective boundary layer ! jp_ext = ibld - imld +1 ! Recalculate ML averages and differences using new ML depth jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 CALL zdf_osm_vertical_average( Kbb, Kmm, & & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, & & jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) #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 ! rotate mean currents and changes onto wind align co-ordinates CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) #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 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< jp_ext(:,:) = 1 ! ag 19/03 CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_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 ! -------------------------------------------------------------------- CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear, & & zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla, & & zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, & & zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Need to put in code for contributions that are applied explicitly to ! the prognostic variables ! 1. Entrainment flux ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! rotate non-gradient velocity terms back to model reference frame DO_2D( 0, 0, 0, 0 ) DO jk = 2, ibld(ji,jj) ztemp = ghamu(ji,jj,jk) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) END DO END_2D ! KPP-style Ri# mixing IF ( ln_kpprimix ) THEN jkflt = jpk DO_2D( 0, 0, 0, 0 ) IF ( ibld(ji,jj) < jkflt ) jkflt = ibld(ji,jj) END_2D DO jk = jkflt+1, jpkm1 ! Shear production at uw- and vw-points (energy conserving form) DO_2D( 1, 0, 1, 0 ) IF ( jk > MIN( ibld(ji,jj), ibld(ji+1,jj) ) ) THEN z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * & & ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) / & & ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) END IF IF ( jk > MIN( ibld(ji,jj), ibld(ji,jj+1) ) ) THEN z2dv(ji,jj) = 0.5_wp * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) * & & ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) / & & ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) END IF END_2D DO_2D( 0, 0, 0, 0 ) IF ( jk > ibld(ji,jj) ) THEN ! Shear prod. at w-point weightened by mask zesh2 = ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & & + ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) ! Local Richardson number zri = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX(zesh2, epsln) zfri = MIN( zri / rn_riinfty , 1.0_wp ) zfri = ( 1.0_wp - zfri * zfri ) zrimix = zfri * zfri * zfri * wmask(ji, jj, jk) zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zrimix*rn_difri ) zviscos(ji,jj,jk) = MAX( zviscos(ji,jj,jk), zrimix*rn_difri ) END IF END_2D END DO END IF ! ln_kpprimix = .true. ! KPP-style set diffusivity large if unstable below BL IF( ln_convmix) THEN DO_2D( 0, 0, 0, 0 ) DO jk = ibld(ji,jj) + 1, jpkm1 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) END DO END_2D 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_2D( 0, 0, 0, 0 ) IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer ! Calculate MLE flux contribution from surface fluxes DO jk = 1, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - 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(ji,jj,jk,Kmm) / 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(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) END DO ELSE ! Surface transports limited to OSBL. ! Viscosity for MLEs DO jk = 1, mld_prof(ji,jj) znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) END DO ENDIF END_2D #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 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids !CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) ! GN 25/8: need to change tmask --> wmask DO_3D( 0, 0, 0, 0, 2, jpkm1 ) p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) END_3D ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) END_3D ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & & ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) #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.*rho0*tmask(:,:,1)*zustar**2*zustke ) ! Stokes drift read in from sbcwave (=2). CASE(2:3) IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & & SQRT(ut0sd**2 + vt0sd**2 ) ) END SELECT IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*zwb0 ) ! IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*zwth0 ) ! 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.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine IF ( iom_use("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("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy 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 IF( ln_timing ) CALL timing_stop('zdf_osm') 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), DIMENSION(jpi,jpj) :: zb_coup, zc_coup_vis, zc_coup_dif ! REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! Coefficients in cubic polynomial specifying diffusivity in pycnocline REAL(wp) :: zznd_ml, zznd_pyc REAL(wp) :: zmsku, zmskv 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 IF( ln_timing ) CALL timing_start('zdf_osm_dv') zb_coup(:,:) = 0.0_wp DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) #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_wp * zvispyc_n_sc(ji,jj) ) zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) ELSE zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/03 zdifpyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/03 zvispyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 IF(lcoup(ji,jj) ) THEN ! ag 19/03 ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub| ! already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90 ! Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub) ! wet-cell averaging .. zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) * & & SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) ) zz_b = -gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / zhml(ji,jj) - zb_coup(ji,jj) ) / & & ( zhml(ji,jj) + zz_b ) ! ag 19/03 #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))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj), & & ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj) WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv FLUSH(narea+100) END IF #endif !#ifdef key_osm_debug ! WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj) ! WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', zz_b, 'pb_coup', zb_coup(ji,jj), & ! & ' pc_coup_vis', zc_coup_vis(ji,jj) ! FLUSH(narea+400) !#endif zz_b = -zhml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / & & zvisml_sc(ji,jj) ! ag 19/03 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / & & zdifml_sc(ji,jj) )**p2third zc_coup_dif(ji,jj) = 0.5_wp * ( -zdifml_sc(ji,jj) / zhml(ji,jj) * ( 1.0_wp - zbeta_d_sc(ji,jj) )**1.5_wp + & & 1.5_wp * ( zdifml_sc(ji,jj) / zhml(ji,jj) ) * zbeta_d_sc(ji,jj) * & & SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b ! ag 19/03 #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 pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj) FLUSH(narea+100) END IF #endif !#ifdef key_osm_debug ! WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', zc_coup_dif(ji,jj) ! FLUSH(narea+400) !#endif ELSE ! ag 19/03 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) / & & ( zdifml_sc(ji,jj) + epsln ) )**p2third ! ag 19/03 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / & & ( zvisml_sc(ji,jj) + epsln ) ! ag 19/03 ENDIF ! ag 19/03 ENDIF ! ag 19/03 #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) WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 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 ELSE zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) #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_2D ! DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN DO jk = 2, imld(ji,jj) ! mixed layer diffusivity zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) ! zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 ! zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & & * ( 1.0 - 0.5 * zznd_ml**2 ) END DO ! Coupling to bottom IF ( lcoup(ji,jj) ) THEN ! ag 19/03 DO jk = mbkt(ji,jj), imld(ji,jj), -1 ! ag 19/03 zz_b = - ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) ! ag 19/03 zviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ! ag 19/03 zdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2 ! ag 19/03 END DO ! ag 19/03 ENDIF ! ag 19/03 ! pycnocline IF ( lpyc(ji,jj) ) THEN ! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed. za_cubic = 0.5 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic DO jk = imld(ji,jj) , ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) ! zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) END DO ! viscosity profiles. za_cubic = 0.5 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj) ) / MAX(zvispyc_n_sc(ji,jj), 1.e-8) zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic DO jk = imld(ji,jj) , ibld(ji,jj) zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) END DO ! IF ( zdhdt(ji,jj) > 0._wp ) THEN ! zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) ! zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) ! ELSE ! zdiffut(ji,jj,ibld(ji,jj)) = 0._wp ! zviscos(ji,jj,ibld(ji,jj)) = 0._wp ! ENDIF ENDIF ELSE ! stable conditions DO jk = 2, ibld(ji,jj) zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) END DO IF ( zdhdt(ji,jj) > 0._wp ) THEN zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) ENDIF ENDIF ! end if ( lconv ) ! END_2D IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) ) ! BBL-coupling velocity scale IF( ln_timing ) CALL timing_stop('zdf_osm_dv') 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. IF( ln_timing ) CALL timing_start('zdf_osm_os') ! Determins stability and set flag lconv DO_2D( 0, 0, 0, 0 ) IF ( zhol(ji,jj) < 0._wp ) THEN lconv(ji,jj) = .TRUE. ELSE lconv(ji,jj) = .FALSE. ENDIF END_2D zekman(A2D(0)) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * zhbl(A2D(0)) / MAX(zustar(A2D(0)), 1.e-8 ) ) zshear(A2D(0)) = 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(A2D(0)) = 1 DO_2D( 0, 0, 0, 0 ) 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.0_wp ) THEN ! Northern hemisphere zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( -1.0_wp * zdv_ml(ji,jj), 1e-5_wp)**2 ) ELSE ! Southern hemisphere zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( zdv_ml(ji,jj), 1e-5_wp)**2 ) END IF 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_wp * MAX( 0.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) > 1e-10 ) THEN IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN ! Growing shear layer j_ddh(ji,jj) = 0 lshear(ji,jj) = .TRUE. ELSE j_ddh(ji,jj) = 1 ! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. lshear(ji,jj) = .TRUE. ! ELSE END IF ELSE j_ddh(ji,jj) = 2 lshear(ji,jj) = .FALSE. END IF ! 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_2D ! Calculate entrainment buoyancy flux due to surface fluxes. DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) IF (nn_osm_SD_reduce > 0 ) THEN ! Effective Stokes drift already reduced from surface value zr_stokes = 1.0_wp ELSE ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, ! requires further reduction where BL is deep zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) END IF zwb_ent(ji,jj) = -2.0_wp * 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_wp * 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_2D zwb_min(:,:) = zlarge DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_shr = -1.0_wp * 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 END IF ! lshear IF ( lconv(ji,jj) ) THEN ! Unstable OSBL zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2.0_wp * zwbav(ji,jj) END IF ! 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_2D IF( ln_timing ) CALL timing_stop('zdf_osm_os') END SUBROUTINE zdf_osm_osbl_state 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 IF( ln_timing ) CALL timing_start('zdf_osm_vr') DO_2D( 0, 0, 0, 0 ) ztemp = zu(ji,jj) zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) ztemp = zdu(ji,jj) zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) END_2D IF( ln_timing ) CALL timing_stop('zdf_osm_vr') 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 IF( ln_timing ) CALL timing_start('zdf_osm_osf') znd_param(A2D(0)) = 0.0_wp DO_2D( 0, 0, 0, 0 ) ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) END_2D DO_2D( 0, 0, 0, 0 ) ! IF ( lconv(ji,jj) ) THEN IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) ! Calculate potential energies of actual profile and reference profile. zpe_mle_layer = 0._wp zpe_mle_ref = 0._wp 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 * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) END DO ! Non-dimensional parameter to diagnose the presence of thermocline znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) ENDIF ENDIF #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_2D ! Diagnosis DO_2D( 0, 0, 0, 0 ) 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_2D IF( ln_timing ) CALL timing_stop('zdf_osm_osf') 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 IF( ln_timing ) CALL timing_start('zdf_osm_eg') DO_2D( 0, 0, 0, 0 ) IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? zbeta = rab_n(ji,jj,1,jp_sal) jkb = jbase(ji,jj) jkb1 = MIN(jkb + 1, mbkt(ji,jj)) zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & & / e3w(ji,jj,jkb1,Kmm) zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & & / e3w(ji,jj,jkb1,Kmm) zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) ELSE zdtdz(ji,jj) = 0._wp zdsdz(ji,jj) = 0._wp zdbdz(ji,jj) = 0._wp END IF END_2D IF( ln_timing ) CALL timing_stop('zdf_osm_eg') END SUBROUTINE zdf_osm_external_gradients SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha ) REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: pdbdz ! Gradients in the pycnocline REAL(wp), DIMENSION(:,:), INTENT( inout ) :: palpha INTEGER :: jk, jj, ji REAL(wp) :: zbgrad REAL(wp) :: zgamma_b_nd, znd REAL(wp) :: zzeta_m REAL(wp), PARAMETER :: ppgamma_b = 2.25_wp ! IF( ln_timing ) CALL timing_start('zdf_osm_pscp') ! DO_2D( 0, 0, 0, 0 ) IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! convective conditions IF ( lpyc(ji,jj) ) THEN zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * zhol(ji,jj) ) ) ) palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) * & & zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / & & ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) ) palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) ztmp = 1.0_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 = palpha(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 = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp IF ( znd <= zzeta_m ) THEN ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) pdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + palpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & & EXP( -6.0_wp * ( 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 ) pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * ppgamma_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_buoyancy_profiles:lconv=lpyc=T',& & 'zzeta_m=', zzeta_m, ' zalpha=', palpha(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_wp ) THEN ! Depth increasing, or steady. IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN IF ( zhol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) zbgrad = zdb_bl(ji,jj) * ztmp DO jk = 2, ibld(ji,jj) znd = gdepw(ji,jj,jk,Kmm) * ztmp pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) END DO ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) zbgrad = zdb_bl(ji,jj) * ztmp DO jk = 2, ibld(ji,jj) znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**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,'(1(a,g11.3))')'end of zdf_osm_pycnocline_buoyancy_profiles:lconv=F zbgrad=', zbgrad ! WRITE(narea+100,'(1(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_2D ! IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) ) END IF ! IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') ! END SUBROUTINE zdf_osm_pycnocline_buoyancy_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_wp REAL(wp), PARAMETER :: zgamma_c = 2.0_wp REAL(wp), PARAMETER :: zdhoh = 0.1_wp REAL(wp), PARAMETER :: zalpha_b = 0.3_wp REAL(wp), PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5 ! also in pycnocline_depth IF( ln_timing ) CALL timing_start('zdf_osm_cd') DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN ! Convective IF ( ln_osm_mle ) THEN IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN ! Fox-Kemper buoyancy flux average over OSBL zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) ELSE zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) ENDIF zvel_max = ( 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) > 1e-15 ) THEN zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0.0_wp ) * zdh(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) zpsi = ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) * & & ( zwb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) ) * zdh(ji,jj) / zhbl(ji,jj) zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) * & & ( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) * MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) zpsi = zalpha_b * MAX( zpsi, 0.0_wp ) zdhdt(ji,jj) = -1.0_wp * ( zwb_ent(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) / & & ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15_wp ) ) + & & zpsi / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-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 = -1.0_wp * a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / & & ( zvel_max + MAX( zdb_bl(ji,jj), 1e-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,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 = -1.0_wp * a_ddh * ( 1.0 - 1.6_wp * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & & ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1e-8_wp ) ) * zddhdt ELSE zddhdt = 0.0_wp ENDIF ! j_ddh zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) * & & zdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-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) = -1.0_wp * MIN( zvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) ENDIF ELSE ! Fox-Kemper not used. zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) ! added ajgn 23 July as temporay fix ENDIF ! ln_osm_mle ELSE ! lconv - Stable zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) IF ( zdhdt(ji,jj) < 0._wp ) THEN ! For long timsteps factor in brackets slows the rapid collapse of the OSBL zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) ELSE zpert = MAX( 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), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 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) = -1.0_wp * MIN( zvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) 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), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 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_2D IF( ln_timing ) CALL timing_stop('zdf_osm_cd') 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 IF( ln_timing ) CALL timing_start('zdf_osm_th') DO_2D( 0, 0, 0, 0 ) #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_Dt / 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) - ts(ji,jj,jm,jp_tem,Kmm) ) & & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & & 0.0 ) + zvel_max IF ( ln_osm_mle ) THEN zhbl_s = zhbl_s + MIN( & & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w(ji,jj,jm,Kmm) ) ELSE zhbl_s = zhbl_s + MIN( & & rn_Dt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & & e3w(ji,jj,jm,Kmm) ) ENDIF ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 #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) - ts(ji,jj,jm,jp_tem,Kmm) )& & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& & 0.0 ) + & & 2.0 * zvstr(ji,jj)**2 / zhbl_s ! Alan is thuis right? I have simply changed hbli to hbl zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) lpyc(ji,jj) = .FALSE. ENDIF IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 #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(ji,jj,4,Kmm) ) ibld(ji,jj) = MAX(jm, 4 ) ELSE ! change zero or one model level. hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) ENDIF zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) #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_2D IF( ln_timing ) CALL timing_stop('zdf_osm_th') 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 IF( ln_timing ) CALL timing_start('zdf_osm_pt') DO_2D( 0, 0, 0, 0 ) IF ( lshear(ji,jj) ) THEN IF ( lconv(ji,jj) ) THEN IF ( zdb_bl(ji,jj) > 1e-15_wp ) THEN IF ( j_ddh(ji,jj) == 0 ) THEN zvel_max = ( zvstr(ji,jj)**3 + 0.5_wp * 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), 1e-15 ) ) zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX( zustar(ji,jj), 1e-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_Dt, 0.625_wp * 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_Dt ) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) ENDIF ELSE ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5_wp * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt ) dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.2_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj) END IF ELSE ! lconv ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here ! boundary layer deepening IF ( zdb_bl(ji,jj) > 0._wp ) THEN ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) ELSE zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! IF(dhdt < 0) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ! IF (dhdt >= 0) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse ENDIF ELSE ! lshear ! for lshear = .FALSE. calculate ddhdt here IF ( lconv(ji,jj) ) THEN IF( ln_osm_mle ) THEN IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE ! unstable zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = zari * hbl(ji,jj) ELSE ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! ln_osm_mle IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ELSE ! unstable zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) ENDIF ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = zari * hbl(ji,jj) ELSE ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) zdh_ref = 0.2 * hbl(ji,jj) ENDIF END IF ! ln_osm_mle dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! Alan: this hml is never defined or used ELSE ! IF (lconv) ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here ! boundary layer deepening IF ( zdb_bl(ji,jj) > 0._wp ) THEN ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) ELSE zdh_ref = 0.2 * hbl(ji,jj) ENDIF ELSE ! IF(dhdt < 0) zdh_ref = 0.2 * hbl(ji,jj) ENDIF ! IF (dhdt >= 0) dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse ENDIF ! IF (lconv) ENDIF ! lshear hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj) - 1,Kmm), 1.e-3) ) , 1 ) imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) #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_2D IF( ln_timing ) CALL timing_stop('zdf_osm_pt') 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 !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('zdf_osm_zhg') ! !== MLD used for MLE ==! mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) ikt = mbkt(ji,jj) zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END_3D DO_2D( 1, 1, 1, 1 ) mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) END_2D ! ensure mld_prof .ge. ibld ! ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation ! ztm(:,:) = 0._wp zsm(:,:) = 0._wp DO_3D( 1, 1, 1, 1, 1, ikmax ) zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) END_3D ! average temperature and salinity. ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) ! calculate horizontal gradients at u & v points zmld_midu(:,:) = 0.0_wp ztsm_midu(:,:,:) = 10.0_wp DO_2D( 0, 0, 1, 0 ) zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) END_2D zmld_midv(:,:) = 0.0_wp ztsm_midv(:,:,:) = 10.0_wp DO_2D( 1, 0, 0, 0 ) zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) END_2D CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm) CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm) DO_2D( 0, 0, 1, 0 ) dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) END_2D DO_2D( 1, 0, 0, 0 ) dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) END_2D DO_2D( 0, 0, 0, 0 ) ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) END_2D IF( ln_timing ) CALL timing_stop('zdf_osm_zhg') END SUBROUTINE zdf_osm_zmld_horizontal_gradients SUBROUTINE zdf_osm_mle_parameters( pmld, 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) :: pmld ! == 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 IF( ln_timing ) CALL timing_start('zdf_osm_mp') ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. DO_2D( 0, 0, 0, 0 ) IF ( lconv(ji,jj) ) THEN ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 ENDIF END_2D ! Timestep mixed layer eddy depth. DO_2D( 0, 0, 0, 0 ) IF ( lmle(ji,jj) ) THEN ! MLE layer growing. ! Buoyancy gradient at base of MLE layer. zthermal = rab_n(ji,jj,1,jp_tem) zbeta = rab_n(ji,jj,1,jp_sal) jkb = mld_prof(ji,jj) jkb1 = MIN(jkb + 1, mbkt(ji,jj)) ! zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) ) zdb_mle = zb_bl(ji,jj) - zbuoy ! Timestep hmle. hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_Dt / zdb_mle ELSE IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau ELSE hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau ENDIF ENDIF hmle(ji,jj) = MAX( MIN( hmle(ji,jj), ht(ji,jj) ), gdepw(ji,jj,4,Kmm) ) 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) = pmld(ji,jj) END_2D mld_prof = 4 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) END_3D DO_2D( 0, 0, 0, 0 ) zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm) END_2D IF( ln_timing ) CALL timing_stop('zdf_osm_mp') END SUBROUTINE zdf_osm_mle_parameters END SUBROUTINE zdf_osm SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, & & knlev, pt, ps, pb, pu, pv, & & kp_ext, pdt, pds, pdb, pdu, pdv ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_vertical_average *** !! !! ** Purpose : Determines vertical averages from surface to knlev, !! and optionally the differences between these vertical !! averages and values at an external level !! !! ** Method : Averages are calculated from the surface to knlev. !! The external level used to calculate differences is !! knlev+kp_ext !!---------------------------------------------------------------------- INTEGER, INTENT(in ) :: Kbb, Kmm ! Ocean time-level indices INTEGER, DIMENSION(jpi,jpj), INTENT(in ) :: knlev ! Number of levels to average over. REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt, ps ! Average temperature and salinity REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pb ! Average buoyancy REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pu, pv ! Average current components INTEGER, DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: kp_ext ! External-level offsets REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pdt, pds ! Difference between average temperature, salinity, REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pdb ! buoyancy, REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pdu, pdv ! velocity components and the OSBL ! INTEGER :: jk, jkflt, jkmax, ji, jj ! Loop indices INTEGER :: ibld_ext ! External-layer index REAL(wp), DIMENSION(jpi,jpj) :: zthick ! Layer thickness REAL(wp) :: zthermal, zbeta ! Thermal/haline expansion/contraction coefficients !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('zdf_osm_va') ! ! Averages over depth of boundary layer pt(A2D(0)) = 0.0_wp ps(A2D(0)) = 0.0_wp pu(A2D(0)) = 0.0_wp pv(A2D(0)) = 0.0_wp zthick(:,:) = epsln jkflt = jpk jkmax = 0 DO_2D( 0, 0, 0, 0 ) IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj) IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) END_2D DO_3D( 0, 0, 0, 0, 2, jkflt ) ! Upper, flat part of layer zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) pt(ji,jj) = pt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) ps(ji,jj) = ps(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) pu(ji,jj) = pu(ji,jj) + e3t(ji,jj,jk,Kmm) * & & ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) / & & MAX( 1.0_wp , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) pv(ji,jj) = pv(ji,jj) + e3t(ji,jj,jk,Kmm) * & & ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) / & & MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) END_3D DO_3D( 0, 0, 0, 0, jkflt+1, jkmax ) ! Lower, non-flat part of layer IF ( knlev(ji,jj) >= jk ) THEN zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) pt(ji,jj) = pt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) ps(ji,jj) = ps(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) pu(ji,jj) = pu(ji,jj) + e3t(ji,jj,jk,Kmm) * & & ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) / & & MAX( 1.0_wp , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) pv(ji,jj) = pv(ji,jj) + e3t(ji,jj,jk,Kmm) * & & ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) / & & MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) END IF END_3D DO_2D( 0, 0, 0, 0 ) pt(ji,jj) = pt(ji,jj) / zthick(ji,jj) ps(ji,jj) = ps(ji,jj) / zthick(ji,jj) pu(ji,jj) = pu(ji,jj) / zthick(ji,jj) pv(ji,jj) = pv(ji,jj) / zthick(ji,jj) zthermal = rab_n(ji,jj,1,jp_tem) ! ideally use ibld not 1?? zbeta = rab_n(ji,jj,1,jp_sal) pb(ji,jj) = grav * zthermal * pt(ji,jj) - grav * zbeta * ps(ji,jj) END_2D ! ! Differences between vertical averages and values at an external layer IF ( PRESENT( kp_ext ) ) THEN DO_2D( 0, 0, 0, 0 ) ibld_ext = knlev(ji,jj) + kp_ext(ji,jj) IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN ! ag 09/03 ! Two external levels are available pdt(ji,jj) = pt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm) pds(ji,jj) = ps(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm) pdu(ji,jj) = pu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) / & & MAX(1.0_wp , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) pdv(ji,jj) = pv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) / & & MAX(1.0_wp , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) zthermal = rab_n(ji,jj,1,jp_tem) ! ideally use ibld not 1?? zbeta = rab_n(ji,jj,1,jp_sal) pdb(ji,jj) = grav * zthermal * pdt(ji,jj) - grav * zbeta * pds(ji,jj) ELSE pdt(ji,jj) = 0.0_wp pds(ji,jj) = 0.0_wp pdu(ji,jj) = 0.0_wp pdv(ji,jj) = 0.0_wp pdb(ji,jj) = 0.0_wp ENDIF END_2D END IF ! IF( ln_timing ) CALL timing_stop('zdf_osm_va') ! END SUBROUTINE zdf_osm_vertical_average SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear, & & pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla, & & pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml, & & pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) !!--------------------------------------------------------------------- !! *** ROUTINE zdf_osm_fgr_terms *** !! !! ** Purpose : Compute non-gradient terms in flux-gradient relationship !! !! ** Method : !! !!---------------------------------------------------------------------- INTEGER, INTENT(in ) :: Kmm ! Time-level index INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer INTEGER, DIMENSION(:,:), INTENT(in ) :: kmld ! ML base layer INTEGER, DIMENSION(:,:), INTENT(in ) :: kp_ext ! Offset for external level LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldpyc ! Pycnocline flags INTEGER, DIMENSION(:,:), INTENT(in ) :: k_ddh ! Type of shear layer REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth REAL(wp), DIMENSION(:,:), INTENT(in ) :: phml ! ML depth REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline depth REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdhdt ! BL depth tendency REAL(wp), DIMENSION(:,:), INTENT(in ) :: phol ! Stability parameter for boundary layer REAL(wp), DIMENSION(:,:), INTENT(in ) :: pshear ! Shear production REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustar ! Friction velocity REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrl ! Langmuir velocity scale REAL(wp), DIMENSION(:,:), INTENT(in ) :: pvstr ! Velocity scale (approaches zustar for large Langmuir number) REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrc ! Convective velocity scale REAL(wp), DIMENSION(:,:), INTENT(in ) :: puw0 ! Surface u-momentum flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwth0 ! Surface heat flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pws0 ! Surface freshwater flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb0 ! Surface buoyancy flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwthav ! BL average heat flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwsav ! BL average freshwater flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwbav ! BL average buoyancy flux REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustke ! Surface Stokes drift REAL(wp), DIMENSION(:,:), INTENT(in ) :: pla ! Langmuir number REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdt_bl ! Temperature diff. between BL average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pds_bl ! Salinity diff. between BL average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_bl ! Velocity diff. (u) between BL average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_bl ! Velocity diff. (u) between BL average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdt_ml ! Temperature diff. between mixed-layer average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pds_ml ! Salinity diff. between mixed-layer average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_ml ! Buoyancy diff. between mixed-layer average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_ml ! Velocity diff. (u) between mixed-layer average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_ml ! Velocity diff. (v) between mixed-layer average and basal value REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdtdz_bl_ext ! External temperature gradients REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdsdz_bl_ext ! External salinity gradients REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdbdz_pyc ! Pycnocline buoyancy gradients REAL(wp), DIMENSION(:,:), INTENT(in ) :: palpha_pyc ! REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdiffut ! t-diffusivity REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pviscos ! Viscosity ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z3ddz_pyc_1, z3ddz_pyc_2 ! Pycnocline gradient/shear profiles ! INTEGER :: ji, jj, jk, jkm_bld, jkf_mld, jkm_mld ! Loop indices #ifdef key_osm_debug INTEGER :: jl, jm #endif INTEGER :: istat ! Memory allocation status REAL(wp) :: zznd_d, zznd_ml, zznd_pyc, znd ! Temporary non-dimensional depths REAL(wp), DIMENSION(A2D(0)) :: zsc_wth_1,zsc_ws_1 ! Temporary scales REAL(wp), DIMENSION(A2D(0)) :: zsc_uw_1, zsc_uw_2 ! Temporary scales REAL(wp), DIMENSION(A2D(0)) :: zsc_vw_1, zsc_vw_2 ! Temporary scales REAL(wp), DIMENSION(A2D(0)) :: ztau_sc_u ! Dissipation timescale at base of WML REAL(wp) :: zbuoy_pyc_sc, zdelta_pyc ! REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale REAL(wp), DIMENSION(A2D(0)) :: za_cubic, zb_cubic ! Coefficients in cubic polynomial specifying REAL(wp), DIMENSION(A2D(0)) :: zc_cubic, zd_cubic ! diffusivity in pycnocline REAL(wp), DIMENSION(A2D(0)) :: zwt_pyc_sc_1, zws_pyc_sc_1 ! REAL(wp), DIMENSION(A2D(0)) :: zzeta_pyc ! REAL(wp) :: zomega, zvw_max ! REAL(wp), DIMENSION(A2D(0)) :: zuw_bse,zvw_bse ! Momentum, heat, and salinity fluxes REAL(wp), DIMENSION(A2D(0)) :: zwth_ent,zws_ent ! at the top of the pycnocline REAL(wp), DIMENSION(A2D(0)) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term REAL(wp) :: ztmp ! REAL(wp) :: ztgrad, zsgrad, zbgrad ! Variables used to calculate pycnocline gradients REAL(wp) :: zugrad, zvgrad ! Variables for calculating pycnocline shear REAL(wp) :: zdtdz_pyc ! Parametrized gradient of temperature in pycnocline REAL(wp) :: zdsdz_pyc ! Parametrised gradient of salinity in pycnocline REAL(wp) :: zdudz_pyc ! u-shear across the pycnocline REAL(wp) :: zdvdz_pyc ! v-shear across the pycnocline !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('zdf_osm_ft') ! ! Auxiliary indices ! ----------------- jkm_bld = 0 jkf_mld = jpk jkm_mld = 0 DO_2D( 0, 0, 0, 0 ) IF ( kbld(ji,jj) > jkm_bld ) jkm_bld = kbld(ji,jj) IF ( kbld(ji,jj) < jkf_mld ) jkf_mld = kbld(ji,jj) IF ( kmld(ji,jj) > jkm_mld ) jkm_mld = kmld(ji,jj) END_2D ! ! Stokes term in scalar flux, flux-gradient relationship ! ------------------------------------------------------ WHERE ( ldconv(A2D(0)) ) zsc_wth_1(:,:) = pwstrl(A2D(0))**3 * pwth0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) zsc_ws_1(:,:) = pwstrl(A2D(0))**3 * pws0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) ELSEWHERE zsc_wth_1(:,:) = 2.0_wp * pwthav(A2D(0)) zsc_ws_1(:,:) = 2.0_wp * pwsav(A2D(0)) ENDWHERE DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( jk <= kmld(ji,jj) ) THEN zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) * & & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) * & & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) END IF ELSE ! Stable conditions IF ( jk <= kbld(ji,jj) ) THEN zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) * & & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) * & & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) END IF END IF ! Check on ldconv END_3D ! 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 ! ! 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 ( ldconv(A2D(0)) ) zsc_uw_1(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / & & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * pla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) zsc_uw_2(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / & & MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & & ( ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) ELSEWHERE zsc_uw_1(:,:) = pustar(A2D(0))**2 zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & & ( pvstr(A2D(0))**2 + epsln ) ENDWHERE DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( jk <= kmld(ji,jj) ) THEN zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) + & & 0.00125_wp * EXP( -1.0_wp * zznd_d ) * zsc_uw_2(ji,jj) ) * & & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65_wp * 0.15_wp * EXP( -1.0_wp * zznd_d ) * & & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_vw_1(ji,jj) END IF ELSE ! Stable conditions IF ( jk <= kbld(ji,jj) ) THEN ! Corrected to ibld zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) * & & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_uw_1(ji,jj) END IF END IF END_3D #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = kmld(ji,jj) - 1; jm = MIN(kbld(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( ldconv(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 ( ldconv(A2D(0)) ) zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / & & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) zsc_ws_1(:,:) = pwbav(A2D(0)) * pws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / & & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) ELSEWHERE zsc_wth_1(:,:) = 0.0_wp zsc_ws_1(:,:) = 0.0_wp ENDWHERE DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( jk <= kmld(ji,jj) ) THEN zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) ! Calculate turbulent time scale zl_c = 0.9_wp * ( 1.0_wp - EXP( -5.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) * & & ( 1.0_wp - EXP( -15.0_wp * ( 1.2_wp - zznd_ml ) ) ) zl_l = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) * & & ( 1.0_wp - EXP( -8.0_wp * ( 1.15_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) ! Non-gradient buoyancy terms ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_wth_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_ws_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) END IF ELSE ! Stable conditions IF ( jk <= kbld(ji,jj) ) THEN 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 IF END IF END_3D DO_2D( 0, 0, 0, 0 ) IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN ztau_sc_u(ji,jj) = phml(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**1.5_wp ) zwth_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj) zws_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN zbuoy_pyc_sc = 2.0_wp * MAX( pdb_ml(ji,jj), 0.0_wp ) / pdh(ji,jj) zdelta_pyc = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird / & & SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) * & & zdelta_pyc**2 / pdh(ji,jj) zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) * & & zdelta_pyc**2 / pdh(ji,jj) zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(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))')'lpyc= lconv=T,dh<0.2*hbl: zbuoy_pyc_sc=',zbuoy_pyc_sc,' zdelta_pyc=',zdelta_pyc WRITE(narea+100,'(3(a,g11.3))')'zwt_pyc_sc_1=',zwt_pyc_sc_1(ji,jj),' zws_pyc_sc_1=',zws_pyc_sc_1(ji,jj), & & ' zzeta_pyc=',zzeta_pyc(ji,jj) FLUSH(narea+100) END IF #endif END IF END IF END_2D DO_3D( 0, 0, 0, 0, 2, jkm_bld ) IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - & & 0.045_wp * ( ( zwth_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * & & MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) ghams(ji,jj,jk) = ghams(ji,jj,jk) - & & 0.045_wp * ( ( zws_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * & & MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) #ifdef key_osm_debug END IF END_3D jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN WRITE(narea+100,'(3(a,g11.3))')'lpyc= lconv=T: ztau_sc_u=',ztau_sc_u(ji,jj),' zwth_ent=',zwth_ent(ji,jj), & & ' zws_ent=',zws_ent(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 ) WRITE(narea+100,*) FLUSH(narea+100) END IF DO_3D( 0, 0, 0, 0, 2, jkm_bld ) IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) #endif IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) .AND. kbld(ji,jj) - kmld(ji,jj) > 3 ) THEN ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp * zws_pyc_sc_1(ji,jj) * & & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird END IF END IF ! End of pycnocline END_3D ! IF(ln_dia_osm) THEN IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! Upward turb. temperature entrainment flux IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! Upward turb. salinity entrainment flux END IF ! zsc_vw_1(:,:) = 0.0_wp WHERE ( ldconv(A2D(0)) ) zsc_uw_1(:,:) = -1.0_wp * pwb0(A2D(0)) * pustar(A2D(0))**2 * phml(A2D(0)) / & & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) zsc_uw_2(:,:) = pwb0(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) / & & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp ) ELSEWHERE zsc_uw_1(:,:) = 0.0_wp ENDWHERE DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( jk <= kmld(ji,jj) ) THEN zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp * & & ( zsc_uw_1(ji,jj) + 0.125_wp * EXP( -0.5_wp * zznd_d ) * & & ( 1.0_wp - EXP( -0.5_wp * zznd_d ) ) * zsc_uw_2(ji,jj) ) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) END IF ELSE ! Stable conditions IF ( jk <= kbld(ji,jj) ) THEN 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 IF ENDIF END_3D ! DO_2D( 0, 0, 0, 0 ) IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN IF ( k_ddh(ji,jj) == 0 ) THEN ! Place holding code. Parametrization needs checking for these conditions. zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) ELSE zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) ENDIF zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj) za_cubic(ji,jj) = zuw_bse(ji,jj) - zb_cubic(ji,jj) zvw_max = 0.7_wp * ff_t(ji,jj) * ( pustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * pustar(ji,jj) * phml(ji,jj) ) zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse(ji,jj) zc_cubic(ji,jj) = zvw_bse(ji,jj) - zd_cubic(ji,jj) END IF END_2D DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld ) ! Need ztau_sc_u to be available. Change to array. IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) * & & ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) * & & ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse(ji,jj) * & & ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) * & & ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) END IF ! ldconv .AND. ldpyc END_3D ! #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = kmld(ji,jj) - 1; jm = MIN(kbld(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( ldconv(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) ] ! ----------------------------------------------------------------------- WHERE ( ldconv(A2D(0)) ) zsc_wth_1(:,:) = pwth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) zsc_ws_1(:,:) = pws0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) WHERE ( ldpyc(A2D(0)) ) ! Pycnocline scales zsc_wth_pyc(:,:) = -0.003_wp * pwstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pdt_ml(A2D(0)) zsc_ws_pyc(:,:) = -0.003_wp * pwstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pds_ml(A2D(0)) END WHERE ELSEWHERE zsc_wth_1(:,:) = 2.0 * pwthav(A2D(0)) zsc_ws_1(:,:) = pws0(A2D(0)) END WHERE DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( ( jk > 1 ) .AND. ( jk <= kmld(ji,jj) ) ) THEN zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) * & & ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) * & & ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) END IF ! ! may need to comment out lpyc block IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN ! Pycnocline zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) END IF ELSE IF( pdhdt(ji,jj) > 0. ) THEN IF ( ( jk > 1 ) .AND. ( jk <= kbld(ji,jj) ) ) THEN zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) + & 7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_wth_1(ji,jj) ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) + & 7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_ws_1(ji,jj) END IF ENDIF ENDIF END_3D ! WHERE ( ldconv(A2D(0)) ) zsc_uw_1(:,:) = pustar(A2D(0))**2 zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) ELSEWHERE zsc_uw_1(:,:) = pustar(A2D(0))**2 zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) * & & zsc_uw_1(:,:) zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phbl(A2D(0)) zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) * & & zsc_vw_1(:,:) ENDWHERE DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) IF ( ldconv(ji,jj) ) THEN IF ( jk <= kmld(ji,jj) ) THEN zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + & & 0.3_wp * ( -2.0_wp + 2.5_wp * ( 1.0_wp + 0.1_wp * zznd_ml**4 ) - EXP( -8.0_wp * zznd_ml ) ) * & & zsc_uw_1(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + & & 0.3_wp * 0.1_wp * ( EXP( -1.0_wp * zznd_d ) + EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) * & & zsc_vw_1(ji,jj) END IF ELSE IF ( jk <= kbld(ji,jj) ) THEN znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) IF ( zznd_d <= 2.0 ) THEN ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp * & & ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) * & & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ) * zsc_uw_1(ji,jj) ELSE ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp * & & ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - znd ) ) ) * zsc_uw_2(ji,jj) ENDIF ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) * & & EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) END IF END IF END_3D #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = kmld(ji,jj) - 1; jm = MIN(kbld(ji,jj) + 2, mbkt(ji,jj) ) WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj) IF (ldpyc(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( ldconv(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_3D( 0, 0, 0, 0, 2, jkm_bld ) IF ( ( .NOT. ldconv(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj) ! ALMG to think about IF ( znd >= 0.0_wp ) THEN ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) ELSE ghamu(ji,jj,jk) = 0.0_wp ghamv(ji,jj,jk) = 0.0_wp ENDIF END IF END_3D ! ! Pynocline contributions ! IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Allocate arrays for output of pycnocline gradient/shear profiles ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) z3ddz_pyc_1(:,:,:) = 0.0_wp z3ddz_pyc_2(:,:,:) = 0.0_wp END IF DO_3D( 0, 0, 0, 0, 2, jkm_bld ) IF ( ldconv (ji,jj) ) THEN ! Unstable conditions. Shouldn;t be needed with no pycnocline code. ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) !Alan is this right? ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & ! & )/ (zdh(ji,jj) + epsln ) ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v ! IF ( znd <= 0.0 ) THEN ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) ! ELSE ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) ! ENDIF ! END DO ELSE ! Stable conditions IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN ! Pycnocline profile only defined when depth steady of increasing. IF ( pdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN IF ( phol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) ztgrad = pdt_bl(ji,jj) * ztmp zsgrad = pds_bl(ji,jj) * ztmp zbgrad = pdb_bl(ji,jj) * ztmp IF ( jk <= kbld(ji,jj) ) THEN znd = gdepw(ji,jj,jk,Kmm) * ztmp zdtdz_pyc = ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) zdsdz_pyc = zsgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc IF ( ln_dia_pyc_scl ) THEN z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc END IF END IF ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) ztgrad = pdt_bl(ji,jj) * ztmp zsgrad = pds_bl(ji,jj) * ztmp zbgrad = pdb_bl(ji,jj) * ztmp IF ( jk <= kbld(ji,jj) ) THEN znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp zdtdz_pyc = ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) zdsdz_pyc = zsgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc IF ( ln_dia_pyc_scl ) THEN z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc END IF END IF ENDIF ! IF (zhol >=0.5) ENDIF ! IF (zdb_bl> 0.) ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero END IF END IF END_3D IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) END IF DO_3D( 0, 0, 0, 0, 2, jkm_bld ) IF ( .NOT. ldconv (ji,jj) ) THEN IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj) zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj) IF ( jk <= kbld(ji,jj) ) THEN znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) IF ( znd < 1.0 ) THEN zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) ELSE zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) ENDIF zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + pviscos(ji,jj,jk) * zdudz_pyc ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + pviscos(ji,jj,jk) * zdvdz_pyc IF ( ln_dia_pyc_shr ) THEN z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc END IF END IF END IF END IF END_3D IF ( ln_dia_pyc_shr ) THEN ! Output of pycnocline shear profiles IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) END IF 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 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Deallocate arrays used for output of pycnocline gradient/shear profiles DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) END IF ! DO_2D( 0, 0, 0, 0 ) ghamt(ji,jj,kbld(ji,jj)) = 0.0_wp ghams(ji,jj,kbld(ji,jj)) = 0.0_wp ghamu(ji,jj,kbld(ji,jj)) = 0.0_wp ghamv(ji,jj,kbld(ji,jj)) = 0.0_wp END_2D #ifdef key_osm_debug IF(narea==nn_narea_db) THEN ji=iloc_db; jj=jloc_db jl = kmld(ji,jj) - 1; jm = MIN(kbld(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("zviscos") ) CALL iom_put( "zviscos", wmask*pviscos ) END IF ! IF( ln_timing ) CALL timing_stop('zdf_osm_ft') ! END SUBROUTINE zdf_osm_fgr_terms SUBROUTINE zdf_osm_init( Kmm ) !!---------------------------------------------------------------------- !! *** ROUTINE zdf_osm_init *** !! !! ** Purpose : Initialization of the vertical eddy diffivity and !! viscosity when using a osm turbulent closure scheme !! !! ** Method : Read the namosm namelist and check the parameters !! called at the first timestep (nit000) !! !! ** input : Namlist namosm !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: Kmm ! time level INTEGER :: ios ! local integer INTEGER :: ji, jj, jk ! dummy loop indices REAL z1_t2 !! NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & #ifdef key_osm_debug & ,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 & ,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 !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('zdf_osm_init') READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) IF(lwm) WRITE ( numond, namzdf_osm ) IF(lwp) THEN ! Control print WRITE(numout,*) WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave SELECT CASE (nn_osm_wave) CASE(0) WRITE(numout,*) ' calculated assuming constant La#=0.3' CASE(1) WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves' CASE(2) WRITE(numout,*) ' calculated from ECMWF wave fields' END SELECT WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce WRITE(numout,*) ' fraction of hbl to average SD over/fit' WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac SELECT CASE (nn_osm_SD_reduce) CASE(0) WRITE(numout,*) ' No reduction' CASE(1) WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL' CASE(2) WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL' END SELECT WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv #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 ! Flags associated with diagnostic output IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) ) ln_dia_pyc_shr = .TRUE. IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE. ! ! allocate zdfosm arrays IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) IF( ln_osm_mle ) THEN ! Initialise Fox-Kemper parametrization READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') IF(lwm) WRITE ( numond, namosm_mle ) IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' WRITE(numout,*) '~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namosm_mle : ' WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit ENDIF ! ENDIF ! IF(lwp) THEN WRITE(numout,*) IF( ln_osm_mle ) THEN WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' ELSE WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' ENDIF ENDIF ! IF( ln_osm_mle ) THEN ! MLE initialisation ! rb_c = grav * rn_osm_mle_rho_c /rho0 ! Mixed Layer buoyancy criteria IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' ! IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! ! ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) ! ENDIF ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) z1_t2 = 2.e-5 DO_2D( 1, 1, 1, 1 ) r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) END_2D ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) ! ENDIF call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle IF( ln_zdfddm) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' Double diffusion mixing on temperature and salinity ' WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm ' ENDIF ENDIF !set constants not in namelist !----------------------------- IF(lwp) THEN WRITE(numout,*) ENDIF IF (nn_osm_wave == 0) THEN dstokes(:,:) = rn_osm_dstokes END IF ! Horizontal average : initialization of weighting arrays ! ------------------- SELECT CASE ( nn_ave ) CASE ( 0 ) ! no horizontal average IF(lwp) WRITE(numout,*) ' no horizontal average on avt' IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !' ! weighting mean arrays etmean ! ( 1 1 ) ! avt = 1/4 ( 1 1 ) ! etmean(:,:,:) = 0.e0 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) etmean(ji,jj,jk) = tmask(ji,jj,jk) & & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) END_3D CASE ( 1 ) ! horizontal average IF(lwp) WRITE(numout,*) ' horizontal average on avt' ! weighting mean arrays etmean ! ( 1/2 1 1/2 ) ! avt = 1/8 ( 1 2 1 ) ! ( 1/2 1 1/2 ) etmean(:,:,:) = 0.e0 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) etmean(ji,jj,jk) = tmask(ji, jj,jk) & & / MAX( 1., 2.* tmask(ji,jj,jk) & & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) END_3D CASE DEFAULT WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave CALL ctl_stop( ctmp1 ) END SELECT ! Initialization of vertical eddy coef. to the background value ! ------------------------------------------------------------- DO jk = 1, jpk avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) END DO ! zero the surface flux for non local term and osm mixed layer depth ! ------------------------------------------------------------------ ghamt(:,:,:) = 0. ghams(:,:,:) = 0. ghamu(:,:,:) = 0. ghamv(:,:,:) = 0. ! IF( ln_timing ) CALL timing_stop('zdf_osm_init') END SUBROUTINE zdf_osm_init SUBROUTINE osm_rst( kt, Kmm, cdrw ) !!--------------------------------------------------------------------- !! *** ROUTINE osm_rst *** !! !! ** Purpose : Read or write BL fields in restart file !! !! ** Method : use of IOM library. If the restart does not contain !! required fields, they are recomputed from stratification !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time step index INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag INTEGER :: id1, id2, id3 ! iom enquiry index INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: iiki, ikt ! local integer REAL(wp) :: zhbf ! tempory scalars REAL(wp) :: zN2_c ! local scalar REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('osm_rst') !!----------------------------------------------------------------------------- ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return !!----------------------------------------------------------------------------- IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) IF( id1 > 0 ) THEN ! 'wn' exists; read CALL iom_get( numror, jpdom_auto, 'wn', ww ) WRITE(numout,*) ' ===>>>> : wn read from restart file' ELSE ww(:,:,:) = 0._wp WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' END IF id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return CALL iom_get( numror, jpdom_auto, 'hbl' , hbl ) CALL iom_get( numror, jpdom_auto, 'dh', dh ) WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' IF( ln_osm_mle ) THEN id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. ) IF( id3 > 0) THEN CALL iom_get( numror, jpdom_auto, 'hmle' , hmle ) WRITE(numout,*) ' ===>>>> : hmle read from restart file' ELSE WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' hmle(:,:) = hbl(:,:) ! Initialise MLE depth. END IF END IF RETURN ELSE ! 'hbl' & 'dh' not in restart file, recalculate WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' END IF END IF !!----------------------------------------------------------------------------- ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return !!----------------------------------------------------------------------------- IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbl into the restart file, then return IF(lwp) WRITE(numout,*) '---- osm-rst ----' CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww ) CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl ) CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh ) IF( ln_osm_mle ) THEN CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle ) END IF RETURN END IF !!----------------------------------------------------------------------------- ! Getting hbl, no restart file with hbl, so calculate from surface stratification !!----------------------------------------------------------------------------- IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' ! w-level of the mixing and mixed layers CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm) imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria ! hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ikt = mbkt(ji,jj) hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level END_3D ! DO_2D( 1, 1, 1, 1 ) iiki = MAX(4,imld_rst(ji,jj)) hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) END_2D WRITE(numout,*) ' ===>>>> : hbl computed from stratification' IF( ln_osm_mle ) THEN hmle(:,:) = hbl(:,:) ! Initialise MLE depth. WRITE(numout,*) ' ===>>>> : hmle set = to hbl' END IF ww(:,:,:) = 0._wp WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' IF( ln_timing ) CALL timing_stop('osm_rst') END SUBROUTINE osm_rst SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_osm *** !! !! ** Purpose : compute and add to the tracer trend the non-local tracer flux !! !! ** Method : ??? !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! time step index INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation ! INTEGER :: ji, jj, jk ! IF( ln_timing ) CALL timing_start('tra_osm') IF( kt == nit000 ) THEN IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF ENDIF IF( l_trdtra ) THEN !* Save ta and sa trends ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) ENDIF DO_3D( 0, 0, 0, 0, 1, jpkm1 ) pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & & - ( ghamt(ji,jj,jk ) & & - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & & - ( ghams(ji,jj,jk ) & & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) END_3D ! save the non-local tracer flux trends for diagnostics IF( l_trdtra ) THEN ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) ENDIF IF(sn_cfctl%l_prtctl) THEN CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ENDIF ! IF( ln_timing ) CALL timing_stop('tra_osm') 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 IF( ln_timing ) CALL timing_start('trc_osm') WRITE(*,*) 'trc_osm: Not written yet', kt IF( ln_timing ) CALL timing_stop('trc_osm') END SUBROUTINE trc_osm SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_osm *** !! !! ** Purpose : compute and add to the velocity trend the non-local flux !! copied/modified from tra_osm !! !! ** Method : ??? !!---------------------------------------------------------------------- INTEGER , INTENT( in ) :: kt ! ocean time step index INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation ! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('dyn_osm') IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF !code saving tracer trends removed, replace with trdmxl_oce DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & & - ( ghamu(ji,jj,jk ) & & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & & - ( ghamv(ji,jj,jk ) & & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) END_3D ! ! code for saving tracer trends removed ! IF( ln_timing ) CALL timing_stop('dyn_osm') END SUBROUTINE dyn_osm !!====================================================================== END MODULE zdfosm