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

Fix reproducibility issues for certain irregular layout when running in repro 64bit mode #57

Merged
merged 4 commits into from
Dec 2, 2024
Merged

Fix reproducibility issues for certain irregular layout when running in repro 64bit mode #57

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
191539e
rearrange vpert to avoid multiple scaling divisions
JosephMouallem Nov 27, 2024
91d6ecd
enforce DP parameters
JosephMouallem Nov 27, 2024
174f695
remove unnecessary temp arrays
JosephMouallem Nov 27, 2024
8c5cdda
ensure DP cutoff values
JosephMouallem Nov 27, 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
56 changes: 25 additions & 31 deletions gsmphys/satmedmfvdiff.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
real(kind=kind_phys) dtdz1(im), gdx(im),
& phih(im), phim(im), prn(im,km-1),
& rbdn(im), rbup(im), thermal(im),
& ustar(im), wstar(im), hpblx(im),
& ust3(im), wst3(im),
& ustar(im), hpblx(im),
& z0(im), crb(im),
& hgamt(im), hgamq(im),
& wscale(im),vpert(im),
& zol(im), sflux(im), radj(im),
& tx1(im), tx2(im)
Expand Down Expand Up @@ -181,28 +179,28 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
!
real(kind=kind_phys) h1
!!
parameter(gravi=1.0/grav)
parameter(gravi=1.d0/grav)
parameter(g=grav)
parameter(gocp=g/cp)
parameter(cont=cp/g,conq=hvap/g,conw=1.0/g) ! for del in pa
parameter(cont=cp/g,conq=hvap/g,conw=1.d0/g) ! for del in pa
! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa
parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
parameter(wfac=7.0,cfac=4.5)
parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
parameter(vk=0.4,rimin=-100.)
parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
parameter(wfac=7.d0,cfac=4.5d0)
parameter(gamcrt=3.d0,gamcrq=0.d0,sfcfrac=0.1d0)
parameter(vk=0.4d0,rimin=-100.0d0)
parameter(rbcr=0.25d0,zolcru=-0.02d0,tdzmin=1.d-3)
!parameter(rlmn=30.,rlmx=500.,elmx=500.)
parameter(prmin=0.25,prmax=4.0,prtke=1.0,prscu=0.67)
parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
parameter(tkmin=1.e-9,dspfac=0.5,dspmax=10.0)
parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
parameter(aphi5=5.,aphi16=16.)
parameter(elmfac=1.0,elefac=1.0,cql=100.)
parameter(dw2min=1.e-4,dkmax=1000.)
parameter(qlcr=3.5e-5,zstblmax=2500.) !,xkzinv=0.15)
parameter(h1=0.33333333)
parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15,ce0=0.4)
parameter(rchck=1.5,cdtn=25.)
parameter(prmin=0.25d0,prmax=4.d0,prtke=1.d0,prscu=0.67d0)
parameter(f0=1.d-4,crbmin=0.15d0,crbmax=0.35d0)
parameter(tkmin=1.d-9,dspfac=0.5d0,dspmax=10.0d0)
parameter(qmin=1.d-8,qlmin=1.d-12,zfmin=1.d-8)
parameter(aphi5=5.0d0,aphi16=16.d0)
parameter(elmfac=1.d0,elefac=1.d0,cql=100.0d0)
parameter(dw2min=1.d-4,dkmax=1000.d0)
parameter(qlcr=3.5d-5,zstblmax=2500.d0) !,xkzinv=0.15)
parameter(h1=1.d0/3.d0)
parameter(ck0=0.4d0,ck1=0.15d0,ch0=0.4d0,ch1=0.15d0,ce0=0.4d0)
parameter(rchck=1.5d0,cdtn=25.d0)

elmx = rlmx
!
Expand Down Expand Up @@ -424,14 +422,14 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k)

if (cap_k0_land) then
if(tem1 > 1.e-5) then
if(tem1 > 1.d-5) then
xkzo(i,k) = min(xkzo(i,k),xkzinv)
xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
endif
else
! kgao note: do not apply upper-limiter over land and sea ice points
! (consistent with change in satmedmfdifq.f in Jun 2020)
if(tem1 > 0. .and. islimsk(i) == 0 ) then
if(tem1 > 0.d0 .and. islimsk(i) == 0 ) then
xkzo(i,k) = min(xkzo(i,k), xkzinv)
xkzmo(i,k) = min(xkzmo(i,k), xkzinv)
endif
Expand Down Expand Up @@ -616,10 +614,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
if(zol(i) < zolcru) then
pcnvflg(i) = .true.
endif
wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i)
wstar(i)= wst3(i)**h1
ust3(i) = ustar(i)**3.
wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1
tem = gotvx(i,1)*sflux(i)*hpbl(i)
wscale(i)=(ustar(i)**3+wfac*vk*tem*sfcfrac)**h1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good optimization!!

ptem = ustar(i)/aphi5
wscale(i) = max(wscale(i),ptem)
endif
Expand All @@ -629,9 +625,7 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
!
do i = 1,im
if(pcnvflg(i)) then
hgamt(i) = heat(i)/wscale(i)
hgamq(i) = evap(i)/wscale(i)
vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1)
vpert(i) = (heat(i) + evap(i)*fv*theta(i,1))/wscale(i)
vpert(i) = max(vpert(i),0.)
tem = min(cfac*vpert(i),gamcrt)
thermal(i)= thermal(i) + tem
Expand Down Expand Up @@ -889,11 +883,11 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke,
do k = 1, km1
do i = 1, im
tem = vk * zl(i,k)
if (zol(i) < 0.) then
if (zol(i) < 0.d0) then
ptem = 1. - 100. * zol(i)
ptem1 = ptem**0.2
zk = tem * ptem1
elseif (zol(i) >= 1.) then
elseif (zol(i) >= 1.d0) then
zk = tem / 3.7
else
ptem = 1. + 2.7 * zol(i)
Expand Down