Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updates from upstream including relative winds #270

Open
wants to merge 6 commits into
base: develop
Choose a base branch
from
Open

Updates from upstream including relative winds #270

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f015361
Deleted redundant line
TillRasmussen Oct 7, 2024
111ee98
Merge pull request #10 from nansencenter/develop
TillRasmussen Oct 13, 2024
9675206
Merge remote-tracking branch 'refs/remotes/origin/develop' into develop
TillRasmussen Oct 13, 2024
a6314b4
Updates from upstream+ a few bugfix. change epmass, ocnscl and hybthk
TillRasmussen Oct 13, 2024
4dbfeb3
Update blkdat to align with changes
TillRasmussen Oct 13, 2024
06caf7a
Remove bug in thermf. Now bit for bit with old version
TillRasmussen Dec 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions TP5a0.06/expt_02.3/blkdat.input
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions TP5a0.06/expt_02.3/blkdat.input_2.3
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
72 changes: 49 additions & 23 deletions hycom/RELO/HYCOM_NERSC_src_v2.3/blkdat.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1830,9 +1844,9 @@ 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)
! --- (=1 MICOM bulk parameterization)
Expand Down Expand Up @@ -1861,12 +1875,26 @@ subroutine blkdat(linit)
endif !1st tile
call blkini(clmflg,'clmflg')
call blkini(wndflg,'wndflg')
if (wndflg.lt.0) then
wndflg = -wndflg
amoflg = 0 !U10
else
amoflg = 1 !U10-Uocn
call blkinr(ocnscl,'ocnscl','(a6," =",f10.4," ")')
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 (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 = wndflg
call blkini(ustflg,'ustflg')
call blkini(flxflg,'flxflg')
call blkini(empflg,'empflg')
Expand Down Expand Up @@ -1934,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 /)') &
Expand Down Expand Up @@ -2180,7 +2199,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
Expand All @@ -2202,7 +2220,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
Expand All @@ -2217,20 +2238,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)')
Expand All @@ -2240,7 +2261,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'
Expand Down Expand Up @@ -2920,3 +2942,7 @@ 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
!> Dec. 2024 - Removed negative wndflg and amoflg due to inclusion of ocnscl
2 changes: 1 addition & 1 deletion hycom/RELO/HYCOM_NERSC_src_v2.3/forfun.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,7 +325,7 @@ subroutine forfuna
enddo !j
endif !stroff:else
!
endif ! windf = .true.
endif ! windf = .true.
!
if (thermo) then
!
Expand Down
133 changes: 129 additions & 4 deletions hycom/RELO/HYCOM_NERSC_src_v2.3/hybgen.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
!
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)') &
Expand Down Expand Up @@ -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
2 changes: 1 addition & 1 deletion hycom/RELO/HYCOM_NERSC_src_v2.3/mod_barotp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))* &
Expand Down
Loading