From f0153619842fb10bc2f9899de18d116ed328f362 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Mon, 7 Oct 2024 14:17:36 +0200 Subject: [PATCH 1/4] Deleted redundant line --- hycom/hycom_ALL/hycom_2.2.72_ALL/bin/hycom_profile_hybgen.f | 1 - 1 file changed, 1 deletion(-) diff --git a/hycom/hycom_ALL/hycom_2.2.72_ALL/bin/hycom_profile_hybgen.f b/hycom/hycom_ALL/hycom_2.2.72_ALL/bin/hycom_profile_hybgen.f index 40d900b5..d7d3e4ab 100644 --- a/hycom/hycom_ALL/hycom_2.2.72_ALL/bin/hycom_profile_hybgen.f +++ b/hycom/hycom_ALL/hycom_2.2.72_ALL/bin/hycom_profile_hybgen.f @@ -305,7 +305,6 @@ subroutine hybgen(temp,saln,th3d,dp,theta,kdm, real, parameter :: hybrlx = 1.0/8.0 !1.0/no. baroclinic time steps c logical mxlkta,thermo -cTill integer i,j,kk,n, lp,nstep, i0,j0,itest,jtest integer i,j,kk,n, lp,nstep, i0,j0,itest,jtest integer klist(1,1) real epsil,onem,tencm,onemm,qonem,thbase From a6314b41d69956db7f8c45b3878dcf9d417441fd Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 13 Oct 2024 12:33:24 +0200 Subject: [PATCH 2/4] Updates from upstream+ a few bugfix. change epmass, ocnscl and hybthk --- hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 | 58 +++++-- hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90 | 2 +- hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90 | 133 +++++++++++++++- .../RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90 | 2 +- .../HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 | 14 +- .../RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 | 144 +++++++++++++++--- hycom/RELO/HYCOM_NERSC_src_v2.3/mxkprf.F90 | 37 +++-- hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrt.F90 | 13 +- hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrtm.F90 | 6 +- hycom/RELO/HYCOM_NERSC_src_v2.3/mxpwp.F90 | 13 +- hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 | 134 +++++++++++----- 11 files changed, 461 insertions(+), 95 deletions(-) diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 index 1a4f0a39..362f0be3 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 @@ -885,6 +885,9 @@ subroutine blkdat(linit) ! --- 'hybthn' = HYBGEN: ratio of layer thicknesses to select the thiner ! (1.0 for original behaviour, i.e. key only on thickness; ! 5.0 to favor the larger density anomaly over thickness) +! --- 'hybthk' = HYBGEN: thick-thin-thick ratio to expand the thin layer +! (0.0 for original behaviour, i.e. don't expand, +! 0.9 for maximum allowed expansion) ! --- 'hybmap' = HYBGEN: remapper flag (0=PCM, 1=PLM, 2=PPM, 3=WENO-like) ! --- 'hybflg' = HYBGEN: generator flag (0=T&S, 1=th&S, 2=th&T) ! --- 'advflg' = thermal advection flag (0=T&S, 1=th&S, 2=th&T) @@ -922,6 +925,7 @@ subroutine blkdat(linit) call blkinr(hybrlx,'hybrlx','(a6," =",f10.4," time steps")') call blkinr(hybiso,'hybiso','(a6," =",f10.4," kg/m^3")') call blkinr(hybthn,'hybthn','(a6," =",f10.4," ")') + call blkinr(hybthk,'hybthk','(a6," =",f10.4," ")') call blkini(hybmap,'hybmap') call blkini(hybflg,'hybflg') call blkini(advflg,'advflg') @@ -969,6 +973,16 @@ subroutine blkdat(linit) call xcstop('(blkdat)') stop '(blkdat)' endif +! + if (hybthk.lt.0.0 .or. hybthk.gt.0.9) then + if (mnproc.eq.1) then + write(lp,'(/ a /)') & + &'error - hybthk must be betweeen 0.0 and 0.9' + call flush(lp) + endif !1st tile + call xcstop('(blkdat)') + stop '(blkdat)' + endif ! if (btrmas .and. .not.btrlfr) then if (mnproc.eq.1) then @@ -1833,6 +1847,7 @@ subroutine blkdat(linit) ! --- (=3 wind speed from wind stress; =4,5 wind stress from wind) ! --- (=4 for COARE 3.0; =5 for COREv2 bulk parameterization) ! --- (4,5 use relative wind U10-Uocn; -4,-5 use absolute wind U10) +! --- 'ocnscl' = scale factor for Uocn in relative wind (0.0: absolute wind) ! --- 'ustflg' = ustar forcing flag (3=input,1,2=wndspd,4=stress) ! --- 'flxflg' = thermal forcing flag (0=none,3=net_flux,1-2,4-6=sst-based) ! --- (=1 MICOM bulk parameterization) @@ -1861,8 +1876,27 @@ subroutine blkdat(linit) endif !1st tile call blkini(clmflg,'clmflg') call blkini(wndflg,'wndflg') - if (wndflg.lt.0) then - wndflg = -wndflg + call blkinr(ocnscl,'ocnscl','(a6," =",f10.4," ")') + if (ocnscl.eq.0.0 .and. wndflg.ge.4) then + if (mnproc.eq.1) then + write(lp,'(/ a /)') & + &'error - ocnscl must be >0.0 for relative winds (wndflg=4,5)' + call flush(lp) + endif !1st tile + call xcstop('(blkdat)') + stop '(blkdat)' + endif + if (ocnscl.ne.0.0 .and. wndflg.lt.4) then + if (mnproc.eq.1) then + write(lp,'(/ a /)') & + &'error - ocnscl must be 0.0 for absolute winds (wndflg<4)' + call flush(lp) + endif !1st tile + call xcstop('(blkdat)') + stop '(blkdat)' + endif + wndflg = abs(wndflg) + if (ocnscl.eq.0.0) then amoflg = 0 !U10 else amoflg = 1 !U10-Uocn @@ -2180,7 +2214,6 @@ subroutine blkdat(linit) ! --- pensol use penetrating solar radiation (input above) ! --- pcipf use E-P forcing (may be redefined in forfun) ! --- priver rivers as a precipitation bogas -! --- epmass treat evap-precip as a mass exchange ! ! --- srelax activate surface salinity climatological nudging ! --- trelax activate surface temperature climatological nudging @@ -2202,7 +2235,10 @@ subroutine blkdat(linit) call blkinl(relax, 'relax ') call blkinl(trcrlx,'trcrlx') call blkinl(priver,'priver') - call blkinl(epmass,'epmass') +! +! --- 'epmass' = E-P mass exchange flag (0=no,1=yes,2=river) + call blkini(epmass,'epmass') +! if (mnproc.eq.1) then write(lp,*) endif !1st tile @@ -2217,20 +2253,20 @@ subroutine blkdat(linit) stop '(blkdat)' endif ! - if (epmass .and. .not.thermo) then + if (epmass.ne.0 .and. .not.thermo) then if (mnproc.eq.1) then write(lp,'(/ a /)') & - &'error - epmass must be .false. for flxflg=0' + &'error - epmass must be 0 for flxflg=0' call flush(lp) endif !1st tile call xcstop('(blkdat)') stop '(blkdat)' endif !!Alex add condition epmass.and.btrlfr - if (epmass .and. .not.btrlfr) then + if (epmass.gt.0 .and. .not.btrlfr) then if (mnproc.eq.1) then write(lp,'(/ a /)') & - 'error - btrlfr must be .true. for epmass=1' + 'error - btrlfr must be .true. for epmass>0' call flush(lp) endif !1st tile call xcstop('(blkdat)') @@ -2240,7 +2276,8 @@ subroutine blkdat(linit) ! --- Strongly discouraged in the coupled model. ! --- epmass removes mass from ocean, but ice on top of ocean is treated as ! --- massless, so... - if (iceflg.eq.2 .and. epmass) then +!!! TILL IS THIS WHAT IT SHOULD BE or epmass .ne. 0???? + if (iceflg.eq.2 .and. epmass .eq. 1) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - epmass should be false when iceflg.eq.2' @@ -2920,3 +2957,6 @@ subroutine blkins(svar,cvar,sdefault,append_dot) !> Sep. 2022 - added hybthn !> Apr. 2023 - added optional dx0k !> July 2023 - added mtracr +!> May 2024 - added epmass=2 for river only mass exchange +!> Aug. 2024 - added ocnscl +!> Sep. 2024 - added hybthk diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90 index 07b74d7f..7a420c69 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90 @@ -325,7 +325,7 @@ subroutine forfuna enddo !j endif !stroff:else ! - endif ! windf = .true. + endif ! windf = .true. ! if (thermo) then ! diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90 index 07f24c30..85b7c3f3 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90 @@ -331,7 +331,8 @@ subroutine hybgenaj(m,n,j ) logical ltop !modify top layer real p_hat,p_hat0,p_hat2,p_hat3,hybrlx, & delt,deltm,dels,delsm,q,qdep,qtr,qts,thkbop, & - zthk,dpthin + zthk,dpthin, & + dp_top,dp_cen,dp_bot,dp_far,dp_inc integer i,k,ka,kp,ktr,fixall,fixlay,nums1d character*12 cinfo ! @@ -1040,8 +1041,128 @@ subroutine hybgenaj(m,n,j ) ! ! --- do not maintain constant thickness, k > fixlay ! - if (th3d(i,j,k,n).gt.theta(i,j,k)+epsil .and. & - k.gt.fixlay+1) then +! TILL The lower pressure should always be higher than the layer above. +! A few exceptions has been found. + dp_top = p(i,j,k) - p(i,j,k-1) + epsil + dp_cen = max(0.0,p(i,j,k+1) - p(i,j,k)) + epsil !NEEDED MODIFIED TILL + if (k.eq.kdm) then + dp_bot = dp_cen !disables thick-thin test + dp_far = dp_cen + elseif (k.eq.kdm-1) then + dp_bot = max(0.0,p(i,j,k+2) - p(i,j,k+1)) + epsil ! NEEDED MODIFIED TILL + dp_far = dp_bot !disables thick-2xthin test + else + dp_bot = p(i,j,k+2) - p(i,j,k+1) + epsil + dp_far = p(i,j,k+3) - p(i,j,k+2) + epsil + endif + if ( dp_cen .lt. hybthk*min(dp_top,dp_bot) ) then +! +! --- thick-thin-thick layers with thin < hybthk*thick +! --- expand layer k with water from layers k-1 and k+1 +! + if (min(dp_top,dp_bot) .lt. 0.2*max(dp_top,dp_bot) ) then +! --- weight by layer thicknesses, increment density is wrong + q = dp_top / (dp_top + dp_bot) + elseif (theta(i,j,k).gt.th3d(i,j,k-1,n) .and. & + theta(i,j,k).lt.th3d(i,j,k+1,n) ) then +! --- weight by actual densities, increment at theta(i,j,k) + q = (theta(i,j,k)-th3d(i,j,k+1,n))/ & + (th3d(i,j,k-1,n)-th3d(i,j,k+1,n)) + else +! --- weight by target densities, increment density is approximate + q = (theta(i,j,k)-theta(i,j,k+1))/(theta(i,j,k-1)-theta(i,j,k+1)) + endif +! --- allow for thick layers getting thinner + if (dp_top .le. dp_bot) then + dp_inc = (hybthk*dp_top - dp_cen)/(1.0+q*hybthk) + else + dp_inc = (hybthk*dp_bot - dp_cen)/(1.0+(1.0-q)*hybthk) + endif + p_hat = p(i,j,k) - q*dp_inc + p(i,j,k)=(1.0-qhrlx(k))*p(i,j,k) + & + qhrlx(k) *p_hat + p_hat = p(i,j,k+1) + (1.0-q)*dp_inc + p(i,j,k+1)=(1.0-qhrlx(k))*p(i,j,k+1) + & + qhrlx(k) *p_hat +! +!diag if (i.eq.itest .and. j.eq.jtest) then +!diag write(lp,'(a,3x,i2.2,4f10.3)') & +!diag 'hybgen, thk,thnDP:',k,dp_top*qonem,dp_cen*qonem,dp_bot*qonem, & +!diag hybthk*min(dp_top,dp_bot)*qonem +!diag write(lp,'(a,3x,i2.2,f10.6,2f10.3,f10.5)') & +!diag 'hybgen, thk,thn Q:',k,q,dp_inc*qonem, & +!diag (dp_inc+dp_cen)*qonem, & +!diag q*th3d(i,j,k-1,n)+(1.0-q)*th3d(i,j,k+1,n) +!diag dp_top = p(i,j, k) - p(i,j,k-1) + epsil +!diag dp_cen = p(i,j, k+1) - p(i,j,k) + epsil +!diag dp_bot = p(i,j,min(kdm+1,k+2)) - p(i,j,k+1) + epsil +!diag write(lp,'(a,3x,i2.2,4f10.3)') & +!diag 'hybgen, thk,thnNP:',k,dp_top*qonem,dp_cen*qonem,dp_bot*qonem, & +!diag hybthk*min(dp_top,dp_bot)*qonem +!diag call flush(lp) +!diag endif !debug +! + elseif ( dp_cen+dp_bot .lt. hybthk*min(dp_top,dp_far) ) then +! +! --- thick-thin-thin-thick layers with 2xthin < hybthk*thick +! --- expand layers k,k+1 with water from layers k-1 and k+2 +! --- this should be rare +! + if (min(dp_top,dp_far) .lt. 0.2*max(dp_top,dp_far) ) then +! --- weight by layer thicknesses, increment density is wrong + q = dp_top / (dp_top + dp_far) + elseif (0.5*(theta(i,j,k)+theta(i,j,k+1)).gt.th3d(i,j,k-1,n) .and. & + 0.5*(theta(i,j,k)+theta(i,j,k+1)).lt.th3d(i,j,k+2,n) ) then +! --- weight by actual densities, increment at 0.5*(theta(i,j,k)+theta(i,j,k+1)) + q = (0.5*(theta(i,j,k)+theta(i,j,k+1))-th3d(i,j,k+2,n)) / & + (th3d(i,j,k-1,n)-th3d(i,j,k+2,n)) + else +! --- weight by target densities, increment density is approximate + q = (0.5*(theta(i,j,k)+theta(i,j,k+1))-theta(i,j,k+2)) / & + (theta(i,j,k-1)-theta(i,j,k+2)) + endif +! --- allow for thick layers getting thinner + if (dp_top .le. dp_far) then + dp_inc = (hybthk*dp_top - (dp_cen+dp_bot))/(1.0+q*hybthk) + else + dp_inc = (hybthk*dp_far - (dp_cen+dp_bot))/(1.0+(1.0-q)*hybthk) + endif + p_hat = p(i,j,k) - q*dp_inc + p(i,j,k)=(1.0-qhrlx(k))*p(i,j,k) + & + qhrlx(k) *p_hat + p_hat = p(i,j,k+2) + (1.0-q)*dp_inc + p(i,j,k+2)=(1.0-qhrlx(k))*p(i,j,k+2) + & + qhrlx(k) *p_hat +! +!diag if (i.eq.itest .and. j.eq.jtest) then +!diag if ( min(dp_top,dp_far) .lt. 0.2*max(dp_top,dp_far)) then +!diag write(lp,'(a,3x,i2.2,4f10.3)') & +!diag 'hybgen, thk,thnX2:',k,dp_top*qonem, & +!diag (dp_cen+dp_bot)*qonem,dp_far*qonem, & +!diag hybthk*min(dp_top,dp_far)*qonem +!diag else +!diag write(lp,'(a,3x,i2.2,4f10.3)') & +!diag 'hybgen, thk,thnD2:',k,dp_top*qonem, & +!diag (dp_cen+dp_bot)*qonem,dp_far*qonem, & +!diag hybthk*min(dp_top,dp_far)*qonem +!diag endif +!diag write(lp,'(a,3x,i2.2,f10.6,2f10.3,f10.5)') & +!diag 'hybgen, thk,thn Q:',k,q,dp_inc*qonem, & +!diag (dp_inc+(dp_cen+dp_bot))*qonem, & +!diag q*th3d(i,j,k-1,n)+(1.0-q)*th3d(i,j,k+2,n) +!diag dp_top = p(i,j,k) - p(i,j,k-1) + epsil +!diag dp_cen = p(i,j,k+1) - p(i,j,k) + epsil +!diag dp_bot = p(i,j,k+2) - p(i,j,k+1) + epsil +!diag dp_far = p(i,j,k+3) - p(i,j,k+2) + epsil +!diag write(lp,'(a,3x,i2.2,4f10.3)') & +!diag 'hybgen, thk,thnNP:',k,dp_top*qonem, & +!diag (dp_cen+dp_bot)*qonem,dp_far*qonem, & +!diag hybthk*min(dp_top,dp_far)*qonem +!diag call flush(lp) +!diag endif !debug +! + elseif (th3d(i,j,k,n).gt.theta(i,j,k)+epsil .and. & + k.gt.fixlay+1) then ! ! --- water in layer k is too dense ! --- try to dilute with water from layer k-1 @@ -1299,10 +1420,12 @@ subroutine hybgenaj(m,n,j ) ! enddo !k vertical coordinate relocation - do k=1,kk-1 + do k=2,kk-1 ! If layer is too thick, move interface up if ((p(i,j,k+1)-p(i,j,k)) .gt. dx0k(k)) then p(i,j,k+1)=p(i,j,k)+dx0k(k) + lcm(k) = .false. !dx0k layers are never PCM + lcm(k-1) = .false. ! !diag if (i.eq.itest .and. j.eq.jtest) then !diag write(lp,'(a,i3.2,f8.2)') & @@ -2970,3 +3093,5 @@ end subroutine hybgen_weno_remap !> July 2023 - added trcflg=801: change in layer thickness due to hybgen !> July 2023 - added trcflg=802: change in layer temperature due to hybgen !> July 2023 - added trcflg=803: change in layer salinity due to hybgen +!> Sep. 2024 - dx0k-ed layers are never remapped with PCM +!> Sep. 2024 - added hybthk diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90 index 910dcaf5..d873e4ef 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90 @@ -727,7 +727,7 @@ subroutine barotp(m,n) ! Barotropic Stokes flow included here ! pbudel = (ubavg(i+1,j,ml)+usdbavg(i+1,j))* & - depthu(i+1,j)*scuy(i+1,j)) & + (depthu(i+1,j)*scuy(i+1,j)) & -(ubavg(i, j,ml)+usdbavg(i, j))* & (depthu(i ,j)*scuy(i, j)) pbvdel = (vbavg(i,j+1,ml)+vsdbavg(i,j+1))* & diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 index 3e7c4986..1c146f73 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 @@ -352,7 +352,7 @@ module mod_cb_arrays ! logical, save :: & btrlfr,btrmas,diagno,thermo,windf,mslprf, & - pcipf,epmass,priver,rivera,kparan,lbmont, & + pcipf,priver,rivera,kparan,lbmont, & relax,srelax,trelax,trcrlx,relaxf,relaxs,relaxt, & locsig,vsigma,hybrid,isopyc,icegln,hybraf,isopcm, & mxl_no,mxlkta,mxlktb,mxlkrt,pensol, & @@ -670,6 +670,7 @@ module mod_cb_arrays ! ---'qhybrlx' = HYBGEN: relaxation coefficient (inverse baroclinic time steps) ! --- 'hybiso' = HYBGEN: Use PCM if layer is within hybiso of target density ! --- 'hybthn' = HYBGEN: ratio of layer thicknesses to select the thiner +! --- 'hybthk' = HYBGEN: thick-thin-thick ratio to expand the thin layer ! --- 'visco2' = deformation-dependent Laplacian viscosity factor ! --- 'visco4' = deformation-dependent biharmonic viscosity factor ! --- 'facdf4' = speed-dependent biharmonic viscosity factor @@ -722,6 +723,7 @@ module mod_cb_arrays ! --- 'thkbot' = thickness of bottom boundary layer (m) ! --- 'sigjmp' = minimum density jump across interfaces (theta units) ! --- 'tmljmp' = equivalent temperature jump across the mixed layer (degC) +! --- 'ocnscl' = scale factor for Uocn in relative wind (0.0 for absolute wind) ! --- 'prsbas' = msl pressure is input field + prsbas (Pa) ! --- 'salmin' = minimum salinity allowed in an isopycnic layer (psu) ! --- 'dx00' = maximum layer thickness minimum, optional (m) @@ -780,6 +782,7 @@ module mod_cb_arrays ! --- 'sstflg' = SST relaxation flag (0=none,1=clim,2=atmos,3=obs) ! --- 'icmflg' = ice mask flag (0=none,1=clim,2=atmos,3=obs) ! --- 'difsmo' = KPROF: number of layers with horiz smooth diff coeffs +! --- 'epmass' = E-P mass exchange flag (0=no,1=yes,2=river) ! #if defined(RELO) real, save, allocatable, dimension(:) :: & @@ -791,7 +794,7 @@ module mod_cb_arrays ! real, save :: & thbase,saln0,baclin,batrop, & - qhybrlx,hybiso,hybthn, & + qhybrlx,hybiso,hybthn,hybthk, & visco2,visco4,veldf2,veldf4,facdf4, & temdf2,temdfc,thkdf2,thkdf4,vertmx,diapyc, & tofset,sofset,dtrate,slip,cb,cbar, & @@ -804,7 +807,7 @@ module mod_cb_arrays thkcdw,thkfrz,tfrz_0,tfrz_s,ticegr,hicemn,hicemx, & dx00,dx00f,dx00x, & dp00,dp00f,dp00x,ds00,ds00f,ds00x,dp00i,isotop, & - oneta0,sigjmp,tmljmp,prsbas,emptgt + oneta0,sigjmp,tmljmp,ocnscl,prsbas,emptgt ! integer, save :: & tsofrq,mixfrq,icefrq,icpfrq,nhybrd,nsigma, & @@ -815,7 +818,7 @@ module mod_cb_arrays iversn,iexpt,jerlv0, & iceflg,ishelf,icmflg,wndflg,amoflg,ustflg, & flxflg,empflg,dswflg,albflg,lwflag,sstflg,sssflg, & - empbal,sssbal, & + epmass,empbal,sssbal, & difsmo,disp_count ! real, parameter :: & @@ -1923,3 +1926,6 @@ end module mod_cb_arrays !> July 2023 - added mtracr and itracr !> Sep. 2023 - initialize si_h to zero !> Jan. 2024 - added pbotmin +!> May 2024 - added epmass=2 for river only mass exchange +!> Aug. 2024 - added ocnscl +!> Sep. 2024 - added hybthk diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 index 3ac79c52..fe44b0f6 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 @@ -90,7 +90,13 @@ subroutine momtum_hs(m,n) ! --- ----------------------------------------- ! logical, parameter :: lpipe_momtum=.false. !usually .false. - +! +#if defined (USE_NUOPC_CESMBETA) + logical, parameter :: cesmbeta =.true. +#else + logical, parameter :: cesmbeta =.false. +#endif +! ! Reduced dragw_rho to half in order to align with CICE real, parameter :: dragw_rho=0.00536d0*1026.0d0 !ice-ocean drag from CICE ! real, parameter :: dragw_rho=0.01072d0*1026.0d0 !ice-ocean drag original @@ -388,7 +394,7 @@ subroutine momtum_hs(m,n) endif !usur,vsur ! if (wndflg.eq.2 .or. wndflg.eq.3) then ! tau on p grid - if (cpl_taux) then + if (cesmbeta .and. cpl_taux) then surtx(i,j) = imp_taux(i,j,1) elseif (natm.eq.2) then surtx(i,j) = taux(i,j,l0)*w0+taux(i,j,l1)*w1 @@ -397,7 +403,7 @@ subroutine momtum_hs(m,n) + taux(i,j,l2)*w2+taux(i,j,l3)*w3 endif ! cpl_taux - if (cpl_tauy) then + if (cesmbeta .and. cpl_tauy) then surty(i,j) = imp_tauy(i,j,1) elseif (natm.eq.2) then surty(i,j) = tauy(i,j,l0)*w0+tauy(i,j,l1)*w1 @@ -491,11 +497,13 @@ subroutine momtum_hs(m,n) if (wndflg.eq.4 .and. flxflg.eq.6) then if (amoflg.ne.0) then ! --- use wind-current in place of wind for everything - samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) +! --- set ocnscl to 1.0 for full relative wind + samo = sqrt( (wndx-ocnscl*usur)**2 +& + (wndy-ocnscl*vsur)**2 ) cdw = 1.0e-3*cd_coarep(samo,vpmx,airt,pair, & temp(i,j,1,n)) - surtx( i,j) = rair*cdw*samo*(wndx-usur) - surty( i,j) = rair*cdw*samo*(wndy-vsur) + surtx( i,j) = rair*cdw*samo*(wndx-ocnscl*usur) + surty( i,j) = rair*cdw*samo*(wndy-ocnscl*vsur) wndocn(i,j) = samo !save for thermf else ! --- use wind for everything @@ -509,9 +517,11 @@ subroutine momtum_hs(m,n) temp(i,j,1,n)) if (amoflg.ne.0) then ! --- use wind-current magnitude and direction for stress - samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) - surtx(i,j) = rair*cdw*samo*(wndx-usur) - surty(i,j) = rair*cdw*samo*(wndy-vsur) +! --- set ocnscl to 1.0 for full relative wind + samo = sqrt( (wndx-ocnscl*usur)**2 +& + (wndy-ocnscl*vsur)**2 ) + surtx(i,j) = rair*cdw*samo*(wndx-ocnscl*usur) + surty(i,j) = rair*cdw*samo*(wndy-ocnscl*vsur) else ! --- use wind for everything surtx(i,j) = rair*cdw*wind*wndx @@ -523,9 +533,11 @@ subroutine momtum_hs(m,n) temp(i,j,1,n)) if (amoflg.ne.0) then ! --- use wind-current magnitude and direction for stress - samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) - surtx(i,j) = rair*cdw*samo*(wndx-usur) - surty(i,j) = rair*cdw*samo*(wndy-vsur) +! --- set ocnscl to 1.0 for full relative wind + samo = sqrt( (wndx-ocnscl*usur)**2 +& + (wndy-ocnscl*vsur)**2 ) + surtx(i,j) = rair*cdw*samo*(wndx-ocnscl*usur) + surty(i,j) = rair*cdw*samo*(wndy-ocnscl*vsur) else ! --- use U10 magnitude and direction for stress surtx(i,j) = rair*cdw*wind*wndx @@ -1077,12 +1089,25 @@ subroutine momtum(m,n) ! logical, parameter :: lpipe_momtum=.false. !usually .false. ! +#if defined(MOMTUM_CFL) + logical, parameter :: momtum_cfl=.true. !set by a CPP macro + !include an explicit CFL limiter +#else + logical, parameter :: momtum_cfl=.false. !usually .false. + !include an explicit CFL limiter +#endif + real, parameter :: clip_min=0.5 !minimum clipping to report +! #if defined(RELO) real, save, allocatable, dimension(:,:) :: & vis2u,vis4u,vis2v,vis4v,vort, & wgtia,wgtib,wgtja,wgtjb, & dl2u,dl2uja,dl2ujb,dl2v,dl2via,dl2vib, & pnk0,pnkp,stresl + real, save, allocatable, dimension(:,:) :: & + uvjclp + real, save, allocatable, dimension(:) :: & + uvkmax #else real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & vis2u,vis4u,vis2v,vis4v,vort, & @@ -1095,6 +1120,10 @@ subroutine momtum(m,n) dl2u,dl2uja,dl2ujb,dl2v,dl2via,dl2vib, & pnk0,pnkp,stresl save /momtumr4/ + real, save, dimension(1-nbdy:jdm+nbdy,1:kdm) :: & + uvjclp + real, save, dimension(1:kdm) :: & + uvkmax=clip_min #endif ! #if defined(STOKES) @@ -1112,10 +1141,9 @@ subroutine momtum(m,n) #endif dvdzu,dudzu,dvdzv,dudzv,dzdx,dzdy ! - real uatvk0,uatvkm,uatvkp,usd_v,vatuk0,vatukm,vatukp,vsd_u, & zbot,ztop - +! #endif /* STOKES */ integer, save :: ifirst ! @@ -1126,6 +1154,7 @@ subroutine momtum(m,n) dmontg,dthstr,dragu,dragv,qdpu,qdpv,dpthin, & dpun,uhm,uh0,uhp,dpvn,vhm,vh0,vhp,sum_m,sum_n real dp12,dp23,dp123,dp3m1,ql1,ql2,ql3 + real cfl,uvclpm,uvclpn,uvkclp(kdm) integer i,ia,ib,j,ja,jb,k,ka,l,mbdy,ktop,kmid,kbot,margin ! ! real*8 wtime @@ -1184,6 +1213,16 @@ subroutine momtum(m,n) pnk0 = 0.0 !r_init pnkp = 0.0 !r_init stresl = 0.0 !r_init + if (momtum_cfl) then + allocate( & + uvjclp(1-nbdy:jdm+nbdy,1:kdm) ) + call mem_stat_add( (jdm+2*nbdy)*kdm ) + uvjclp = 0.0 + allocate( & + uvkmax(1:kdm) ) + call mem_stat_add( kdm ) + uvkmax = clip_min + endif !cfl endif !vis2u #if defined(STOKES) if (.not.allocated(usdflux)) then @@ -2955,17 +2994,32 @@ subroutine momtum(m,n) ! --- substitute depth-weighted averages at massless grid points. ! --- extract barotropic velocities generated during most recent ! --- baroclinic time step and use them to force barotropic flow field. +! + if (momtum_cfl) then + uvjclp(:,:) = 0.0 + endif ! margin = 0 ! -!$OMP PARALLEL DO PRIVATE(j,i,k,q,sum_m,sum_n) & +!$OMP PARALLEL DO PRIVATE(j,i,k,q,sum_m,sum_n,cfl,uvclpm,uvclpn) & !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin +! --- limit is conservative, so clipping occurs before CFL is exceedd + cfl = 0.707*0.5*min(scpx(i,j),scpy(i,j))/delt1 if (SEA_U) then k=1 u(i,j,k,m) = u(i,j,k,m)/max(dpu(i,j,k,m),dpthin) u(i,j,k,n) = u(i,j,k,n)/max(dpu(i,j,k,n),dpthin) + if (momtum_cfl) then !this should not be needed + uvclpm = abs(u(i,j,k,m)) + uvclpn = abs(u(i,j,k,n)) + u(i,j,k,m) = max( -cfl, min( cfl, u(i,j,k,m) ) ) + u(i,j,k,n) = max( -cfl, min( cfl, u(i,j,k,n) ) ) + uvclpm = uvclpm - abs(u(i,j,k,m)) + uvclpn = uvclpn - abs(u(i,j,k,n)) + uvjclp(j,k) = max( uvjclp(j,k), uvclpm, uvclpn) + endif !cfl sum_m = u(i,j,k,m)* dpu(i,j,k,m) sum_n = u(i,j,k,n)* dpu(i,j,k,n) do k=2,kk @@ -2977,6 +3031,15 @@ subroutine momtum(m,n) q = min(dpu(i,j,k,n),cutoff) u(i,j,k,n) = (u(i,j,k,n)*q+u(i,j,k-1,n)*(cutoff-q))* & qcutoff + if (momtum_cfl) then !this should not be needed + uvclpm = abs(u(i,j,k,m)) + uvclpn = abs(u(i,j,k,n)) + u(i,j,k,m) = max( -cfl, min( cfl, u(i,j,k,m) ) ) + u(i,j,k,n) = max( -cfl, min( cfl, u(i,j,k,n) ) ) + uvclpm = uvclpm - abs(u(i,j,k,m)) + uvclpn = uvclpn - abs(u(i,j,k,n)) + uvjclp(j,k) = max( uvjclp(j,k), uvclpm, uvclpn) + endif !cfl sum_m = sum_m + u(i,j,k,m)*dpu(i,j,k,m) sum_n = sum_n + u(i,j,k,n)*dpu(i,j,k,n) enddo !k @@ -2995,6 +3058,15 @@ subroutine momtum(m,n) k=1 v(i,j,k,m) = v(i,j,k,m)/max(dpv(i,j,k,m),dpthin) v(i,j,k,n) = v(i,j,k,n)/max(dpv(i,j,k,n),dpthin) + if (momtum_cfl) then !this should not be needed + uvclpm = abs(v(i,j,k,m)) + uvclpn = abs(v(i,j,k,n)) + v(i,j,k,m) = max( -cfl, min( cfl, v(i,j,k,m) ) ) + v(i,j,k,n) = max( -cfl, min( cfl, v(i,j,k,n) ) ) + uvclpm = uvclpm - abs(v(i,j,k,m)) + uvclpn = uvclpn - abs(v(i,j,k,n)) + uvjclp(j,k) = max( uvjclp(j,k), uvclpm, uvclpn) + endif !cfl sum_m = v(i,j,1,m)* dpv(i,j,1,m) sum_n = v(i,j,1,n)* dpv(i,j,1,n) do k=2,kk @@ -3006,6 +3078,15 @@ subroutine momtum(m,n) q = min(dpv(i,j,k,n),cutoff) v(i,j,k,n) = (v(i,j,k,n)*q+v(i,j,k-1,n)*(cutoff-q))* & qcutoff + if (momtum_cfl) then !this should not be needed + uvclpm = abs(v(i,j,k,m)) + uvclpn = abs(v(i,j,k,n)) + v(i,j,k,m) = max( -cfl, min( cfl, v(i,j,k,m) ) ) + v(i,j,k,n) = max( -cfl, min( cfl, v(i,j,k,n) ) ) + uvclpm = uvclpm - abs(v(i,j,k,m)) + uvclpn = uvclpn - abs(v(i,j,k,n)) + uvjclp(j,k) = max( uvjclp(j,k), uvclpm, uvclpn) + endif !cfl sum_m = sum_m + v(i,j,k,m)*dpv(i,j,k,m) sum_n = sum_n + v(i,j,k,n)*dpv(i,j,k,n) enddo !k @@ -3022,6 +3103,30 @@ subroutine momtum(m,n) enddo !i enddo !j - do 30 loop !$OMP END PARALLEL DO +! +! --- check for velocity clipping at cfl limit +! + if (momtum_cfl .and. mod(nstep,3).eq.0) then !skip some time steps for efficiency + do k= 1,kk + uvkclp(k) = 0.0 + do j=1,jj + uvkclp(k) = max( uvkclp(k), uvjclp(j,k) ) + enddo !j + enddo !k + call xcmaxr(uvkclp(1:kk)) + do k= 1,kk + if (uvkclp(k).gt.uvkmax(k)) then + if (mnproc.eq.1) then + write(lp, & + '(i9,a,i3,a,f7.2,a)') & + nstep,' k=',k, & + ' velocty clipped, max=',uvkclp(k),' m/s' + endif !mnproc + uvkmax(k) = uvkclp(k) + endif + enddo !k + call xcsync(flush_lp) + endif !every 3 time steps ! if (lpipe .and. lpipe_momtum) then ! --- compare two model runs. @@ -5401,8 +5506,9 @@ subroutine momtum4(m,n) !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin +! --- limit is conservative, so clipping occurs before CFL is exceedd + cfl = 0.707*0.5*min(scpx(i,j),scpy(i,j))/delt1 if (SEA_U) then - cfl = 0.5*min(scpx(i,j),scpy(i,j))/delt1 k=1 u(i,j,k,m) = u(i,j,k,m)/max(dpu(i,j,k,m),dpthin) u(i,j,k,n) = u(i,j,k,n)/max(dpu(i,j,k,n),dpthin) @@ -5438,11 +5544,8 @@ subroutine momtum4(m,n) ! utotn(i,j) = dt1inv*(sum_n - ubavg(i,j,n)) !for barotp endif !iu - enddo !i ! - do i=1-margin,ii+margin if (SEA_V) then - cfl = 0.5*min(scpx(i,j),scpy(i,j))/delt1 k=1 v(i,j,k,m) = v(i,j,k,m)/max(dpv(i,j,k,m),dpthin) v(i,j,k,n) = v(i,j,k,n)/max(dpv(i,j,k,n),dpthin) @@ -5627,3 +5730,6 @@ end module mod_momtum !> Sep. 2019 - added momtum_init !> Oct. 2019 - placed momtum4_cfl in a CPP macro !> Nov. 2019 - added amoflg for U10 vs U10-Uocn +!> Mar. 2023 - added momtum_cfl in a CPP macro +!> Dec. 2023 - add cesmbeta as a master switch to cpl_ +!> Aug. 2024 - replace U10-Uocn with U10-ocnscl*Uocn diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkprf.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkprf.F90 index fa5ae6df..f81ca44b 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkprf.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkprf.F90 @@ -696,13 +696,17 @@ subroutine mxkppaij(m,n, i,j) sflux1=surflx(i,j)-sswflx(i,j) dtemp=(sflux1+(1.-swfrac(k+1))*sswflx(i,j))* & delt1*g*qspcifh*(qdpmm(k)*qoneta) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln= salflx(i,j)* & - delt1*g* (qdpmm(k)*qoneta) + delt1*g*(qdpmm(k)*qoneta) + elseif (epmass.eq.2) then !river only is mass flux + dsaln=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g*(qdpmm(k)*qoneta) else !water flux treated as a virtual salt flux dsaln=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & - delt1*g* (qdpmm(k)*qoneta) - endif + delt1*g*(qdpmm(k)*qoneta) + endif !epmass !diag if (i.eq.itest.and.j.eq.jtest) then !diag write (lp,101) nstep,i+i0,j+j0,k, & !diag 1.0,swfrac(k+1),dtemp,dsaln @@ -2224,13 +2228,17 @@ subroutine mxmyaij(m,n, i,j) sflux1=surflx(i,j)-sswflx(i,j) dtemp=(sflux1+(1.-swfrac(k+1))*sswflx(i,j))* & delt1*g*qspcifh*(qdpmm(k)*qoneta) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln= salflx(i,j)* & - delt1*g* (qdpmm(k)*qoneta) + delt1*g*(qdpmm(k)*qoneta) + elseif (epmass.eq.2) then !river only is mass flux + dsaln=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g*(qdpmm(k)*qoneta) else !water flux treated as a virtual salt flux dsaln=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & - delt1*g* (qdpmm(k)*qoneta) - endif + delt1*g*(qdpmm(k)*qoneta) + endif !epmass !diag if (i.eq.itest .and. j.eq.jtest) then !diag write (lp,102) nstep,i+i0,j+j0,k, & !diag 0.,1.-swfrac(k+1),dtemp,dsaln @@ -2469,13 +2477,17 @@ subroutine mxgissaij(m,n, i,j) sflux1=surflx(i,j)-sswflx(i,j) dtemp=(sflux1+(1.-swfrac(k+1))*sswflx(i,j))* & delt1*g*qspcifh*(qdpmm(k)*qoneta) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln= salflx(i,j)* & - delt1*g* (qdpmm(k)*qoneta) + delt1*g*(qdpmm(k)*qoneta) + elseif (epmass.eq.2) then !river only is mass flux + dsaln=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g*(qdpmm(k)*qoneta) else !water flux treated as a virtual salt flux dsaln=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & - delt1*g* (qdpmm(k)*qoneta) - endif + delt1*g*(qdpmm(k)*qoneta) + endif !epmass !diag if (i.eq.itest.and.j.eq.jtest) then !diag write (lp,101) nstep,i+i0,j+j0,k, & !diag 0.,1.-swfrac(k+1),dtemp,dsaln @@ -3876,3 +3888,4 @@ subroutine wscale(i,j,zlevel,dnorm,bflux,wm,ws,isb) !> Dec 2018 - added /* USE_NUOPC_CESMBETA */ macro and riv_input !> Mar 2023 - added /* MASSLESS_1MM */ macro !> July 2023 - detrain negative near-surface salinitites +!> May 2024 - added epmass=2 for river only mass exchange diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrt.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrt.F90 index 73309e4b..74543c3b 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrt.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrt.F90 @@ -454,9 +454,13 @@ subroutine mxkrtaaj(m,n, j, depnew) sflux1=surflx(i,j)-sswflx(i,j) dtemp(i)=(sflux1+(1.-swfrac)*sswflx(i,j))* & delt1*g/(spcifh*thk1ta) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln(i)= salflx(i,j)* & delt1*g/thk1ta + elseif (epmass.eq.2) then !river only is mass flux + dsaln(i)=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g/thk1ta else !water flux treated as a virtual salt flux dsaln(i)=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & delt1*g/thk1ta @@ -470,9 +474,13 @@ subroutine mxkrtaaj(m,n, j, depnew) ! dtemp(i)=surflx(i,j)* & delt1*g/(spcifh*thk1ta) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln(i)= salflx(i,j)* & delt1*g/thk1ta + elseif (epmass.eq.2) then !river only is mass flux + dsaln(i)=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g/thk1ta else !water flux treated as a virtual salt flux dsaln(i)=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & delt1*g/thk1ta @@ -1249,3 +1257,4 @@ subroutine mxkrtbbj(m,n, j) !> Nov. 2018 - added wtrflx, salflx now only actual salt flux !> Nov. 2018 - allow for wtrflx in buoyancy flux !> Nov. 2018 - allow for oneta in swfrac and surface fluxes +!> May 2024 - added epmass=2 for river only mass exchange diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrtm.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrtm.F90 index 241fb3f2..658f63fe 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrtm.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxkrtm.F90 @@ -419,8 +419,11 @@ subroutine mxkrtmbj(m,n, sdot, j) ! do i=1,ii if (SEA_P) then - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux vsflx(i)= salflx(i,j) + elseif (epmass.eq.2) then !river only is mass flux + vsflx(i)=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n)) else !water flux treated as a virtual salt flux vsflx(i)=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n)) endif @@ -667,3 +670,4 @@ subroutine mxkrtmbj(m,n, sdot, j) !> May 2014 - use land/sea masks (e.g. ip) to skip land !> Aug. 2018 - added wtrflx, salflx now only actual salt flux !> Nov. 2018 - allow for wtrflx in buoyancy flux +!> May 2024 - added epmass=2 for river only mass exchange diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxpwp.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxpwp.F90 index 46983090..ba05e783 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mxpwp.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mxpwp.F90 @@ -289,9 +289,13 @@ subroutine mxpwpaij(m,n, i,j) sflux1=surflx(i,j)-sswflx(i,j) dtemp=(sflux1+(1.-swfrac(k+1))*sswflx(i,j))* & delt1*g*qoneta/(spcifh*max(onemm,dp1d(k))) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln= salflx(i,j)* & delt1*g*qoneta/max(onemm,dp1d(k)) + elseif (epmass.eq.2) then !river only is mass flux + dsaln=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g*qoneta/max(onemm,dp1d(k)) else !water flux treated as a virtual salt flux dsaln=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & delt1*g*qoneta/max(onemm,dp1d(k)) @@ -305,9 +309,13 @@ subroutine mxpwpaij(m,n, i,j) else !.not.pensol dtemp=surflx(i,j)* & delt1*g*qoneta/(spcifh*max(onemm,dp1d(k))) - if (epmass) then !only actual salt flux + if (epmass.eq.1) then !only actual salt flux dsaln= salflx(i,j)* & delt1*g*qoneta/max(onemm,dp1d(k)) + elseif (epmass.eq.2) then !river only is mass flux + dsaln=(salflx(i,j)- & + (wtrflx(i,j)-rivflx(i,j))*saln(i,j,1,n))* & + delt1*g*qoneta/max(onemm,dp1d(k)) else !water flux treated as a virtual salt flux dsaln=(salflx(i,j)-wtrflx(i,j)*saln(i,j,1,n))* & delt1*g*qoneta/max(onemm,dp1d(k)) @@ -811,3 +819,4 @@ subroutine mlbdep(t1d,s1d,th1d,tr1d,u1d,v1d,p1d,dp1d,kmlb,kmax) !> May 2014 - use land/sea masks (e.g. ip) to skip land !> Aug. 2018 - added wtrflx, salflx now only actual salt flux !> Nov. 2018 - allow for oneta in swfrac and surface fluxes +!> May 2024 - added epmass=2 for river only mass exchange diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 index e48b0739..5bda32af 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 @@ -40,14 +40,14 @@ subroutine thermf_oi(m,n) real*8 d1,d2,d3,d4,ssum(2),s1(2) ! if (ishelf.ne.0) then !sea ice and an ice shelf -!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then if (ishlf(i,j).eq.1) then !standard ocean point if ( cpl_swflx .and. cpl_lwmdnflx .and. cpl_lwmupflx & .and. cpl_precip .and. cesmbeta ) then - sstflx(i,j) = (1.0-covice(i,j))*sstflx(i,j) !relax over ocean surflx(i,j) = surflx(i,j) + flxice(i,j) !ocn/ice frac. in coupler salflx(i,j) = sflice(i,j) + & !ocn/ice frac. in coupler @@ -55,7 +55,6 @@ subroutine thermf_oi(m,n) wtrflx(i,j) = wtrflx(i,j) + wflice(i,j) + & !ocn/ice frac. in coupler rivflx(i,j) !rivers evrywhere else - sswflx(i,j) = (1.0-covice(i,j))*sswflx(i,j) + & !ocean fswice(i,j) !ice cell average surflx(i,j) = (1.0-covice(i,j))*surflx(i,j) + & !ocean @@ -79,7 +78,7 @@ subroutine thermf_oi(m,n) endif !ip enddo !i enddo !j -!$OMP END PARALLEL DO +!$OMP END PARALLEL DO elseif (iceflg.ne.0) then ! switch off the sssflx under sea ice if (sss_underice==.false.) then @@ -99,8 +98,9 @@ subroutine thermf_oi(m,n) endif if ( cpl_swflx .and. cpl_lwmdnflx .and. cpl_lwmupflx & .and. cpl_precip .and. cesmbeta ) then -! --- allow for sea ice -!$OMP PARALLEL DO PRIVATE(j,i) +! --- allow for sea ice +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then @@ -115,9 +115,10 @@ subroutine thermf_oi(m,n) endif enddo !i enddo -!$OMP END PARALLEL DO +!$OMP END PARALLEL DO else -!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then @@ -134,12 +135,13 @@ subroutine thermf_oi(m,n) util1(i,j) = surflx(i,j)*scp2(i,j) util2(i,j) = wtrflx(i,j)*scp2(i,j) endif - enddo !i - enddo !j -!$OMP END PARALLEL DO + enddo !i + enddo !j +!$OMP END PARALLEL DO endif ! .not. cpl else !no sea ice -!$OMP PARALLEL DO PRIVATE(j,i) +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then !sswflx, surflx and sstflx unchanged @@ -150,10 +152,42 @@ subroutine thermf_oi(m,n) endif !ip enddo !i enddo !j -!$OMP END PARALLEL DO +!$OMP END PARALLEL DO endif !ishlf:iceflg ! - if (epmass .and. empbal.eq.1) then + if (epmass.eq.2) then +! +! --- balance rivflx, via a region-wide offset: +! +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) + do j=1,jj + do i=1,ii + if (SEA_P) then + util3(i,j) = rivflx(i,j)*scp2(i,j) + endif + enddo !i + enddo !j + call xcsum(d2, util3,ipa) + d2 = d2/area + q = -d2 +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) + do j=1,jj + do i=1,ii + if (SEA_P) then + rivflx(i,j) = rivflx(i,j) + q + endif + enddo !i + enddo !j + if (ldebug_empbal .and. & + mnproc.eq.1 .and. mod(nstep,100).le.1) then + write (lp,'(i9,a,1pe12.5)') & + nstep,' offset rivflx by ',q + endif !master + endif +! + if (epmass.eq.1 .and. empbal.eq.1) then ! ! --- balance wtrflx to emptgt, via a region-wide offset: ! @@ -177,16 +211,23 @@ subroutine thermf_oi(m,n) nstep,' offset wtrflx by ',q endif !master endif !not already balanced - elseif (.not.epmass .and. empbal.eq.1) then + elseif (epmass.ne.1 .and. empbal.eq.1) then ! ! --- balance wtrflx to emptgt, via a region-wide offset: ! --- virtual salt flux case, balance sss*wtrflx ! +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then - util3(i,j) = wtrflx(i,j)*saln(i,j,1,n)*scp2(i,j) - util4(i,j) = saln(i,j,1,n)*scp2(i,j) + if (epmass.eq.0) then + util3(i,j) = wtrflx(i,j)*saln(i,j,1,n)*scp2(i,j) + else !river only as mass exchange + util3(i,j) = (wtrflx(i,j)-rivflx(i,j))* & + saln(i,j,1,n)*scp2(i,j) + endif + util4(i,j) = saln(i,j,1,n)*scp2(i,j) endif enddo !i enddo !j @@ -211,13 +252,15 @@ subroutine thermf_oi(m,n) nstep,' offset wtrflx by ',q endif !master endif !not already balanced - elseif (epmass .and. empbal.eq.2) then + elseif (epmass.eq.1 .and. empbal.eq.2) then ! ! --- balance wtrflx to emptgt, by scaling down +ve or -ve anomalies ! call xcsum(d2, util2,ipa) !total d2 = d2/area !basin-wide average if (d2.ne.emptgt) then !not balanced at emptgt +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then @@ -270,18 +313,25 @@ subroutine thermf_oi(m,n) endif !master endif !reduce -ve:reduce +ve endif !not already balanced - elseif (.not.epmass .and. empbal.eq.2) then + elseif (epmass.ne.1 .and. empbal.eq.2) then ! ! --- balance wtrflx to emptgt, by scaling down +ve or -ve anomalies ! --- virtual salt flux case, balance sss*wtrflx ! +!$OMP PARALLEL DO PRIVATE(j,i) & +!$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then - util2(i,j) = wtrflx(i,j)*saln(i,j,1,n)*scp2(i,j) - util3(i,j) = max(wtrflx(i,j)-emptgt,0.0)* & - saln(i,j,1,n)*scp2(i,j) - util4(i,j) = saln(i,j,1,n)*scp2(i,j) + if (epmass.eq.0) then + util5(i,j) = wtrflx(i,j) + else !river only as mass exchange + util5(i,j) = wtrflx(i,j)-rivflx(i,j) + endif + util2(i,j) = util5(i,j)*saln(i,j,1,n)*scp2(i,j) + util3(i,j) = max(util5(i,j)-emptgt,0.0)* & + saln(i,j,1,n)*scp2(i,j) + util4(i,j) = saln(i,j,1,n)*scp2(i,j) endif enddo !i enddo !j @@ -299,8 +349,8 @@ subroutine thermf_oi(m,n) do j=1,jj do i=1,ii if (SEA_P) then - if (wtrflx(i,j).lt.emptgt) then - wtrflx(i,j) = emptgt + q*(wtrflx(i,j)-emptgt) + if (util5(i,j).lt.emptgt) then + wtrflx(i,j) = emptgt + q*(util5(i,j)-emptgt) endif util2(i,j) = wtrflx(i,j)*scp2(i,j) endif @@ -319,8 +369,8 @@ subroutine thermf_oi(m,n) do j=1,jj do i=1,ii if (SEA_P) then - if (wtrflx(i,j).gt.emptgt) then - wtrflx(i,j) = emptgt + q*(wtrflx(i,j)-emptgt) + if (util5(i,j).gt.emptgt) then + wtrflx(i,j) = emptgt + q*(util5(i,j)-emptgt) endif util2(i,j) = wtrflx(i,j)*scp2(i,j) endif @@ -333,13 +383,13 @@ subroutine thermf_oi(m,n) endif !master endif !reduce -ve:reduce +ve endif !empbal - +! !!Alex New calculation of epmass, with E-P applied to the top layer !!AJW Modified new epmass calculation to use actual dp rather than h !!AJW updates pbavg, dp and S.1 !!AJW When E-P>0 and layer 1 is thin, some may get merged deeper ! - if (epmass) then !requires btrlfr=.true., see blkdat.F + if (epmass.gt.0) then !requires btrlfr=.true., see blkdat.F !$OMP PARALLEL DO PRIVATE(j,i,k,k1,emnp,dpemnp,dpemnp1,dplay1, & !$OMP onetanew,onetaold,pbanew,pbaold,q) do j=1,jj @@ -350,7 +400,11 @@ subroutine thermf_oi(m,n) ! --- This only works if pbavg and dp have just been integrated from ! --- the same time, hence btrlfr=false is not allowed ! - emnp = -wtrflx(i,j)*svref !m/s + if (epmass.eq.1) then + emnp = -wtrflx(i,j)*svref !m/s + else !river only + emnp = -rivflx(i,j)*svref !m/s + endif pbaold = pbavg(i,j,n) pbanew = pbaold - emnp*delt1*onem !emnp mass = rho0*vol if (.false. .and. i.eq.itest.and.j.eq.jtest) then @@ -361,9 +415,11 @@ subroutine thermf_oi(m,n) pbaold - emaxfr*(pbot(i,j) + pbaold) ) onetaold = 1.0 + pbaold/pbot(i,j) onetanew = 1.0 + pbanew/pbot(i,j) - !some evap may have been discarded - emnp = (pbaold - pbanew)/(delt1*onem) - wtrflx(i,j) = -emnp/svref + if (epmass.eq.1) then + !some evap may have been discarded + emnp = (pbaold - pbanew)/(delt1*onem) + wtrflx(i,j) = -emnp/svref + endif if (.false. .and. i.eq.itest.and.j.eq.jtest) then write(lp,'(i9,a,1p5g12.4)') & nstep,' pba,emnp = ',pbanew*qonem, & @@ -1324,7 +1380,6 @@ subroutine thermfj(m,n,dtime, j) +imp_lwdflx(i,j,1) & +imp_lwuflx(i,j,1) elseif (natm.eq.2) then - radfl=(radflx(i,j,l0)*w0+radflx(i,j,l1)*w1) else radfl=(radflx(i,j,l0)*w0+radflx(i,j,l1)*w1 & @@ -1341,7 +1396,6 @@ subroutine thermfj(m,n,dtime, j) else radfl = radfl - sb_cst*(temp(i,j,1,n)+tzero)**4 + swfl endif - sstflx(i,j) = 0.0 elseif (lwflag.gt.0) then ! --- over-ocean longwave correction to radfl (Qsw+Qlw). @@ -1521,7 +1575,6 @@ subroutine thermfj(m,n,dtime, j) if (empflg.lt.0) then evape = slat*clh*wind*(0.97*qsatur(esst)-vpmx) endif - ! --- Latent Heat flux (W/m2) if(cesmbeta .and. cpl_latflx) then evap=imp_latflx(i,j,1) @@ -1789,7 +1842,7 @@ subroutine thermfj(m,n,dtime, j) ce_n10 = 34.6 *cd_n10_rt*1.e-3 !L-Y eqn. 6b stab = 0.5 + sign(0.5,airt-temp(i,j,1,n)) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt*1.e-3 !L-Y eqn. 6c - +! cd10 = cd_n10 !first guess for exchange coeff's at z ch10 = ch_n10 ce10 = ce_n10 @@ -1827,7 +1880,7 @@ subroutine thermfj(m,n,dtime, j) ce10 = ce_n10/(1.0+ce_n10*xx/cd_n10_rt) * & sqrt(cd10/cd_n10) !L-Y 10c end do !it_a - +! ! --- Latent Heat flux slat = evaplh*rair if (empflg.lt.0) then @@ -1838,7 +1891,7 @@ subroutine thermfj(m,n,dtime, j) else evap = slat*ce10*wind*(qsatur5(temp(i,j,1,n),qrair)-vpmx) endif - +! ! --- Sensible Heat flux ssen = cpcore*rair if(cesmbeta .and. cpl_sensflx) then @@ -1846,7 +1899,7 @@ subroutine thermfj(m,n,dtime, j) else snsibl = ssen*ch10*wind*(temp(i,j,1,n)-airt) endif - +! ! --- Total surface fluxes surflx(i,j) = radfl - snsibl - evap elseif (flxflg.eq.3) then @@ -2386,3 +2439,4 @@ subroutine swfrml_ij(akpar,hbl,bot,zzscl,jerlov,swfrml) !> Dec. 2023 - add cesmbeta as a master switch to cpl_ !> Jan. 2024 - evap with epmass==1 can extend below a thin enough layer 1 !> Jan. 2024 - evap with epmass==1 may be clipped for small oneta +!> May 2024 - added epmass=2 for river only mass exchange From 4dbfeb36a745e5f033df8cd1c832c45a57b56e69 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 13 Oct 2024 12:42:37 +0200 Subject: [PATCH 3/4] Update blkdat to align with changes --- TP5a0.06/expt_02.3/blkdat.input | 2 ++ TP5a0.06/expt_02.3/blkdat.input_2.3 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/TP5a0.06/expt_02.3/blkdat.input b/TP5a0.06/expt_02.3/blkdat.input index 79cce4d3..bc8dee9d 100644 --- a/TP5a0.06/expt_02.3/blkdat.input +++ b/TP5a0.06/expt_02.3/blkdat.input @@ -106,6 +106,7 @@ S-Z(15-11): dp00/f/x/i=3m/1.125/12m/1m; ds=1m/1.125/4m; src_2.2.12; 8.0 'hybrlx' = HYBGEN: inverse relaxation coefficient (time steps) 0.01 'hybiso' = HYBGEN: Use PCM if layer is within hybiso of target density 1.0 'hybthn' = HYBGEN: ratio of layer thicknesses to select the thiner + 0.0 'hybthk' = HYBGEN: thick-thin-thick ratio to expand the thin layer 3 'hybmap' = hybrid remapper flag (0=PCM, 1=PLM, 2=PPM) 0 'hybflg' = hybrid generator flag (0=T&S, 1=th&S, 2=th&T) 0 'advflg' = thermal advection flag (0=T&S, 1=th&S, 2=th&T) @@ -205,6 +206,7 @@ S-Z(15-11): dp00/f/x/i=3m/1.125/12m/1m; ds=1m/1.125/4m; src_2.2.12; 0.0 'tid_t0' = TIDES: origin for ramp time (model day) 12 'clmflg' = climatology frequency flag (6=bimonthly, 12=monthly) 4 'wndflg' = wind stress input flag (0=none,1=u/v-grid,2,3=p-grid,4,5=wnd10m) + 1.0 'ocnscl' = scale factor for Uocn in relative wind (0.0: absolute wind) 2 'ustflg' = ustar forcing flag (3=input,1,2=wndspd,4=stress) 6 'flxflg' = thermal forcing flag (0=none,3=net_flux,1-2,4-6=sst-based) 6 'empflg' = E-P forcing flag (0=none,3=net_E-P, 1-2,4-6=sst-based_E) diff --git a/TP5a0.06/expt_02.3/blkdat.input_2.3 b/TP5a0.06/expt_02.3/blkdat.input_2.3 index f44a21ed..3b785b77 100755 --- a/TP5a0.06/expt_02.3/blkdat.input_2.3 +++ b/TP5a0.06/expt_02.3/blkdat.input_2.3 @@ -106,6 +106,7 @@ S-Z(15-11): dp00/f/x/i=3m/1.125/12m/1m; ds=1m/1.125/4m; src_2.2.12; 8.0 'hybrlx' = HYBGEN: inverse relaxation coefficient (time steps) 0.01 'hybiso' = HYBGEN: Use PCM if layer is within hybiso of target density 1.0 'hybthn' = HYBGEN: ratio of layer thicknesses to select the thiner + 0.0 'hybthk' = HYBGEN: thick-thin-thick ratio to expand the thin layer 3 'hybmap' = hybrid remapper flag (0=PCM, 1=PLM, 2=PPM) 0 'hybflg' = hybrid generator flag (0=T&S, 1=th&S, 2=th&T) 0 'advflg' = thermal advection flag (0=T&S, 1=th&S, 2=th&T) @@ -205,6 +206,7 @@ S-Z(15-11): dp00/f/x/i=3m/1.125/12m/1m; ds=1m/1.125/4m; src_2.2.12; 0.0 'tid_t0' = TIDES: origin for ramp time (model day) 12 'clmflg' = climatology frequency flag (6=bimonthly, 12=monthly) 4 'wndflg' = wind stress input flag (0=none,1=u/v-grid,2,3=p-grid,4,5=wnd10m) + 1.0 'ocnscl' = scale factor for Uocn in relative wind (0.0: absolute wind) 2 'ustflg' = ustar forcing flag (3=input,1,2=wndspd,4=stress) 6 'flxflg' = thermal forcing flag (0=none,3=net_flux,1-2,4-6=sst-based) 6 'empflg' = E-P forcing flag (0=none,3=net_E-P, 1-2,4-6=sst-based_E) From 06caf7a15ff0580b03a5964f328fb6c7348ac065 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Mon, 9 Dec 2024 10:13:35 +0100 Subject: [PATCH 4/4] Remove bug in thermf. Now bit for bit with old version --- hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 | 28 +++++-------------- .../HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 | 4 +-- .../RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 | 21 +++++++------- hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 | 2 +- 4 files changed, 21 insertions(+), 34 deletions(-) diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 index 362f0be3..cf80808e 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90 @@ -1844,9 +1844,8 @@ subroutine blkdat(linit) ! ! --- 'clmflg' = climatology frequency flag (6=bimonthly,12=monthly) ! --- 'wndflg' = wind stress input flag (0=none,1=uv-grid,2,3=p-grid,4,5=wnd10m) -! --- (=3 wind speed from wind stress; =4,5 wind stress from wind) +! --- (=3 wind speed from wind stress; =4,5 relative wind stress from wind U10-Uocn) ! --- (=4 for COARE 3.0; =5 for COREv2 bulk parameterization) -! --- (4,5 use relative wind U10-Uocn; -4,-5 use absolute wind U10) ! --- 'ocnscl' = scale factor for Uocn in relative wind (0.0: absolute wind) ! --- 'ustflg' = ustar forcing flag (3=input,1,2=wndspd,4=stress) ! --- 'flxflg' = thermal forcing flag (0=none,3=net_flux,1-2,4-6=sst-based) @@ -1877,11 +1876,11 @@ subroutine blkdat(linit) call blkini(clmflg,'clmflg') call blkini(wndflg,'wndflg') call blkinr(ocnscl,'ocnscl','(a6," =",f10.4," ")') - if (ocnscl.eq.0.0 .and. wndflg.ge.4) then + if (wndflg.lt.0 .or. wndflg.gt.5) then if (mnproc.eq.1) then - write(lp,'(/ a /)') & - &'error - ocnscl must be >0.0 for relative winds (wndflg=4,5)' - call flush(lp) + write(lp,'(/ a /)') & + &'error - wndflg must be between 0 and 5' + call flush(lp) endif !1st tile call xcstop('(blkdat)') stop '(blkdat)' @@ -1895,12 +1894,7 @@ subroutine blkdat(linit) call xcstop('(blkdat)') stop '(blkdat)' endif - wndflg = abs(wndflg) - if (ocnscl.eq.0.0) then - amoflg = 0 !U10 - else - amoflg = 1 !U10-Uocn - endif + wndflg = wndflg call blkini(ustflg,'ustflg') call blkini(flxflg,'flxflg') call blkini(empflg,'empflg') @@ -1968,15 +1962,6 @@ subroutine blkdat(linit) call xcstop('(blkdat)') stop '(blkdat)' endif - if (wndflg.lt.0 .or. wndflg.gt.5) then - if (mnproc.eq.1) then - write(lp,'(/ a /)') & - &'error - wndflg must be between 0 and 5' - call flush(lp) - endif !1st tile - call xcstop('(blkdat)') - stop '(blkdat)' - endif if (ustflg.lt.1 .or. ustflg.gt.4) then if (mnproc.eq.1) then write(lp,'(/ a /)') & @@ -2960,3 +2945,4 @@ subroutine blkins(svar,cvar,sdefault,append_dot) !> May 2024 - added epmass=2 for river only mass exchange !> Aug. 2024 - added ocnscl !> Sep. 2024 - added hybthk +!> Dec. 2024 - Removed negative wndflg and amoflg due to inclusion of ocnscl diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 index 1c146f73..265a9344 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_cb_arrays.F90 @@ -769,7 +769,6 @@ module mod_cb_arrays ! --- 'iceflg' = sea ice model flag (0=none,1=energy loan,2=coupled/esmf) ! --- 'ishelf' = ice shelf flag (0=none,1=ice shelf over ocean) ! --- 'wndflg' = wind stress input flag (0=none,1=uv-grid,2,3=p-grid,4,5=wnd10m) -! --- 'amoflg' = relative wind flag (0=U10:wndflg=-4,-5;1=U10-Uocn:wndflg=4,5) ! --- 'ustflg' = ustar forcing flag (3=input,1=wndspd,2=stress) ! --- 'flxflg' = thermal forcing flag (0=none,3=net-flux,1,2,4=sst-based) ! --- 'empflg' = E-P forcing flag (0=none,3=net_E-P, 1,2,4=sst-based_E) @@ -816,7 +815,7 @@ module mod_cb_arrays ntracr,mtracr,trcflg(mxtrcr),itracr(801:899), & clmflg,dypflg,iniflg,lbflag,mapflg,yrflag,sshflg, & iversn,iexpt,jerlv0, & - iceflg,ishelf,icmflg,wndflg,amoflg,ustflg, & + iceflg,ishelf,icmflg,wndflg,ustflg, & flxflg,empflg,dswflg,albflg,lwflag,sstflg,sssflg, & epmass,empbal,sssbal, & difsmo,disp_count @@ -1929,3 +1928,4 @@ end module mod_cb_arrays !> May 2024 - added epmass=2 for river only mass exchange !> Aug. 2024 - added ocnscl !> Sep. 2024 - added hybthk +!> Dec. 2024 - removed amoflg diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 index fe44b0f6..06f09256 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/mod_momtum.F90 @@ -370,7 +370,7 @@ subroutine momtum_hs(m,n) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_P) then - if (amoflg.ne.0 .or. & + if (ocnscl.ne.0.0 .or. & (iceflg.eq.2 .and. si_c(i,j).gt.0.0)) then ! --- average currents over top thkcdw meters thksur = onem*min( thkcdw, depths(i,j) ) @@ -495,7 +495,7 @@ subroutine momtum_hs(m,n) endif ! if (wndflg.eq.4 .and. flxflg.eq.6) then - if (amoflg.ne.0) then + if (ocnscl.ne.0.0) then ! --- use wind-current in place of wind for everything ! --- set ocnscl to 1.0 for full relative wind samo = sqrt( (wndx-ocnscl*usur)**2 +& @@ -508,14 +508,14 @@ subroutine momtum_hs(m,n) else ! --- use wind for everything cdw = 1.0e-3*cd_coarep(wind,vpmx,airt,pair, & - temp(i,j,1,n)) + temp(i,j,1,n)) surtx( i,j) = rair*cdw*wind*wndx surty( i,j) = rair*cdw*wind*wndy - endif !amoflg + endif !ocnscl elseif (wndflg.eq.4) then cdw = 1.0e-3*cd_coare(wind,vpmx,airt, & - temp(i,j,1,n)) - if (amoflg.ne.0) then + temp(i,j,1,n)) + if (ocnscl.ne.0.0) then ! --- use wind-current magnitude and direction for stress ! --- set ocnscl to 1.0 for full relative wind samo = sqrt( (wndx-ocnscl*usur)**2 +& @@ -526,12 +526,12 @@ subroutine momtum_hs(m,n) ! --- use wind for everything surtx(i,j) = rair*cdw*wind*wndx surty(i,j) = rair*cdw*wind*wndy - endif !amoflg + endif !ocnscl else ! wndflg.eq.5 ! --- vpmx assumed to contain specific humidity cdw = 1.0e-3*cd_core2(wind,vpmx,airt, & temp(i,j,1,n)) - if (amoflg.ne.0) then + if (ocnscl.ne.0.0) then ! --- use wind-current magnitude and direction for stress ! --- set ocnscl to 1.0 for full relative wind samo = sqrt( (wndx-ocnscl*usur)**2 +& @@ -539,10 +539,10 @@ subroutine momtum_hs(m,n) surtx(i,j) = rair*cdw*samo*(wndx-ocnscl*usur) surty(i,j) = rair*cdw*samo*(wndy-ocnscl*vsur) else -! --- use U10 magnitude and direction for stress +! --- use U10 magnitude and direction for stress surtx(i,j) = rair*cdw*wind*wndx surty(i,j) = rair*cdw*wind*wndy - endif !amoflg + endif endif endif !wndflg @@ -5733,3 +5733,4 @@ end module mod_momtum !> Mar. 2023 - added momtum_cfl in a CPP macro !> Dec. 2023 - add cesmbeta as a master switch to cpl_ !> Aug. 2024 - replace U10-Uocn with U10-ocnscl*Uocn +!> Dec. 2024 - Replace amoflg with ocnscl diff --git a/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 b/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 index 5bda32af..0e052ab1 100644 --- a/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 +++ b/hycom/RELO/HYCOM_NERSC_src_v2.3/thermf.F90 @@ -1314,7 +1314,7 @@ subroutine thermfj(m,n,dtime, j) if (SEA_P) then if (flxflg.gt.0) then ! --- wind = wind, or wind-ocean, speed (m/s) - if (flxflg.eq.6 .and. amoflg.ne.0) then + if (flxflg.eq.6 .and. ocnscl.ne.0.0) then wind=wndocn(i,j) !magnitude of wind minus ocean current elseif(cesmbeta .and. cpl_wndspd) then wind=imp_wndspd(i,j,1)