From b68d4c0c9d3671066e72d1defba84f7d87f539b5 Mon Sep 17 00:00:00 2001 From: JosephMouallem <75327266+JosephMouallem@users.noreply.github.com> Date: Mon, 2 Dec 2024 08:44:25 -0500 Subject: [PATCH] Fix reproducibility issues for certain irregular layout when running in repro 64bit mode (#57) * rearrange vpert to avoid multiple scaling divisions * enforce DP parameters * remove unnecessary temp arrays * ensure DP cutoff values --- gsmphys/satmedmfvdiff.f | 56 ++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 31 deletions(-) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index 4c7e32f4..be5073f9 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -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) @@ -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 ! @@ -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 @@ -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 ptem = ustar(i)/aphi5 wscale(i) = max(wscale(i),ptem) endif @@ -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 @@ -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)