Skip to content

Codex audit of single precision code + Reduce avoidable single-precis… #92

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

Open
rafaelluizmesquita wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Codex audit of single precision code + Reduce avoidable single-precis… #92

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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: 1 addition & 1 deletion src/Mars.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1875,7 +1875,7 @@ subroutine dustprofile(PSF, PTROP, PLEV, TAUCUM, TAUREF, L_LEVELS, TauTot, ConrN

!C Calculate the Reference Pressure Grid (prdst)

refpr = (5.0*psf/ptrop)**(1.0/(float(npdst) - 1.0))
refpr = (5.0d0*psf/ptrop)**(1.0d0/(dble(npdst) - 1.0d0))
prdst(1) = ptrop

do n = 2, npdst
Expand Down
195 changes: 102 additions & 93 deletions src/ModHwm14.f90
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion src/Venus.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -854,7 +854,7 @@ subroutine nlte_tcool(iBlock)
call suaviza ( aux2gcmd, n_gcm, 1, auxgcmd )

do i=1,n_gcm
q15umco2_gcm(iLon,iLat,i) = sngl( aux2gcmd(i) )
q15umco2_gcm(iLon,iLat,i) = dble(aux2gcmd(i))
enddo

enddo !-------- END OF MAIN LONGITUDE LOOP
Expand Down
10 changes: 5 additions & 5 deletions src/apex_more.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -636,21 +636,21 @@ subroutine cofrm(date)
n = n + 1
endif ! n < m
lm = ll + l
if (m == 0) f0 = f0*float(n)/2.
if (m == 0) f0 = f0*dble(n)/2.0d0
if (m == 0) f = f0/sqrt(2.0)
nn = n + 1
mm = 1

if (m /= 0) then
f = f/sqrt(float(n - m + 1)/float(n + m))
f = f/sqrt(dble(n - m + 1)/dble(n + m))
gb(l + 1) = (tc*gh(lm) + t*gh(lm + nc))*f
else
gb(l + 1) = (tc*gh(lm) + t*gh(lm + nc))*f0
endif
gv(l + 1) = gb(l + 1)/float(nn)
gv(l + 1) = gb(l + 1)/dble(nn)
if (m /= 0) then
gb(l + 2) = (tc*gh(lm + 1) + t*gh(lm + nc + 1))*f
gv(l + 2) = gb(l + 2)/float(nn)
gv(l + 2) = gb(l + 2)/dble(nn)
l = l + 2
else
l = l + 1
Expand Down Expand Up @@ -783,7 +783,7 @@ subroutine feldg(iflag, glat, glon, alt, bnrth, beast, bdown, babs)

100 continue
il = ih - i
f = 2./float(i - k + 2)
f = 2.0d0/dble(i - k + 2)
x = xi(1)*f
y = xi(2)*f
z = xi(3)*(f + f)
Expand Down
20 changes: 10 additions & 10 deletions src/apex_routines.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)

do i=mid*10.0-100,mid*10.0+100

ang = real(i)/10.0*dtor
ang = dble(i)/10.0d0*dtor
stfcpa = ste*ctp*cos(ang)-cte*stp
stfspa = sin(ang)*ste

Expand Down Expand Up @@ -166,7 +166,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)

do iLon = 0, 360, 10

gLonGuess = real(iLon)
gLonGuess = dble(iLon)
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)

Expand All @@ -179,7 +179,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)

enddo

gLon = real(iLonBest)
gLon = dble(iLonBest)

dLon = 20.0

Expand All @@ -194,7 +194,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)

do iLat = -5,5

gLatGuess = gLat + real(iLat)/10 * dLat
gLatGuess = gLat + dble(iLat)/10.0d0*dLat
gLonGuess = gLon
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)
Expand All @@ -214,7 +214,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)
do iLon = -5,5

gLatGuess = gLat
gLonGuess = gLon + real(iLon)/10 * dLat
gLonGuess = gLon + dble(iLon)/10.0d0*dLat
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)

Expand All @@ -229,7 +229,7 @@ subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)

enddo

dLon = dLon * exp((real(iCount)-5)/10)
dLon = dLon * exp((dble(iCount)-5.0d0)/10.0d0)
iCount = iCount + 1

enddo
Expand Down Expand Up @@ -2853,7 +2853,7 @@ SUBROUTINE SUBSOLR (IYR,IDAY,IHR,IMN,SEC,SBSLLAT,SBSLLON)
G0 = -2.472 + (-.2558905*(YR-4*NLEAP) - 3.79617E-2*NLEAP)
C
C Universal time in seconds:
UT = FLOAT(IHR*3600 + IMN*60) + SEC
UT = DBLE(IHR*3600 + IMN*60) + SEC
C
C Days (including fraction) since 12 UT on January 1 of IYR:
DF = (UT/86400. - 1.5) + IDAY
Expand All @@ -2877,7 +2877,7 @@ SUBROUTINE SUBSOLR (IYR,IDAY,IHR,IMN,SEC,SBSLLAT,SBSLLON)
SINLAM = SIN(LAMRAD)
C
C Days (including fraction) since 12 UT on January 1 of 2000:
N = DF + FLOAT(365*YR + NLEAP)
N = DF + DBLE(365*YR + NLEAP)
C
C Obliquity of ecliptic:
EPSILON = 23.439 - 4.E-7*N
Expand All @@ -2895,7 +2895,7 @@ SUBROUTINE SUBSOLR (IYR,IDAY,IHR,IMN,SEC,SBSLLAT,SBSLLON)
C Equation of time (degrees):
ETDEG = L - ALPHA
NROT = NINT(ETDEG/360.)
ETDEG = ETDEG - FLOAT(360*NROT)
ETDEG = ETDEG - DBLE(360*NROT)
C
C Apparent time (degrees):
APTIME = UT/240. + ETDEG
Expand All @@ -2904,7 +2904,7 @@ SUBROUTINE SUBSOLR (IYR,IDAY,IHR,IMN,SEC,SBSLLAT,SBSLLON)
C Subsolar longitude:
SBSLLON = 180. - APTIME
NROT = NINT(SBSLLON/360.)
SBSLLON = SBSLLON - FLOAT(360*NROT)
SBSLLON = SBSLLON - DBLE(360*NROT)
C
RETURN
END
83 changes: 41 additions & 42 deletions src/apexsh.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,16 +148,16 @@ module apxshmodule
real(8), parameter :: dtor = pi/180D0, pid2 = pi/2D0, twopi = 2D0*pi
real(8), parameter :: Req = 6378.1370D0, eps = 1.D0/298.257223563D0
real(8), parameter :: Re = Req*(1 - eps/3D0), ecc2 = eps*(2 - eps)
real(4), parameter :: missing = -9999E0
real(8), parameter :: missing = -9999D0

! Req = Equatorial radius of Earth in km (WGS84 value)
! eps = flatness of ellipsoidal Earth (WGS84 value)
! Re = Mean radius of Earth in km
! ecc2 = squared eccentricity of ellipsoidal Earth

character(128) :: datafile
real(4) :: epoch
real(4) :: altlastq, altlastg
real(8) :: epoch
real(8) :: altlastq, altlastg
logical :: loadflag = .true.

end module apxshmodule
Expand All @@ -171,7 +171,7 @@ subroutine loadapxsh(datafilenew, epochnew)
implicit none

character(128) :: datafilenew, datafilelast = ''
real :: epochnew, epochlast = -999.0
real(8) :: epochnew, epochlast = -999.0d0
real(8) :: we0, we1
integer(4) :: iepoch0, iepoch1, iterm, icoord

Expand Down Expand Up @@ -211,8 +211,8 @@ subroutine loadapxsh(datafilenew, epochnew)
gcoeff0(iterm, icoord) = we0*coeff0(iterm, iepoch0, icoord + 3) + we1*coeff0(iterm, iepoch1, icoord + 3)
enddo
enddo
altlastq = -999.0
altlastg = -999.0
altlastq = -999.0d0
altlastg = -999.0d0
endif

!UPDATE LOAD VARIABLES
Expand Down Expand Up @@ -302,10 +302,10 @@ subroutine apxg2q(glat, glon, alt, vecflagin, qlatout, qlonout, f1, f2, f)

implicit none

real(4), intent(in) :: glat, glon, alt
real(8), intent(in) :: glat, glon, alt
integer(4), intent(in) :: vecflagin
real(4), intent(out) :: qlatout, qlonout
real(4), intent(out) :: f1(1:2), f2(1:2), f
real(8), intent(out) :: qlatout, qlonout
real(8), intent(out) :: f1(1:2), f2(1:2), f

integer(4) :: i, l, iterm, itermsh
real(8) :: theta, phi
Expand Down Expand Up @@ -371,8 +371,8 @@ subroutine apxg2q(glat, glon, alt, vecflagin, qlatout, qlonout, f1, f2, f)
qlat = datan2(zq, dsqrt(xq*xq + yq*yq))
cosqlat = dcos(qlat)
sinqlat = dsin(qlat)
qlonout = sngl(qlon/dtor)
qlatout = sngl(qlat/dtor)
qlonout = dble(qlon/dtor)
qlatout = dble(qlat/dtor)

!BASE VECTOR CALCULATIONS
if (vecflag .ne. 0) then
Expand All @@ -399,10 +399,10 @@ subroutine apxg2q(glat, glon, alt, vecflagin, qlatout, qlonout, f1, f2, f)
qlongrad(2) = (-sinqlon*xqgrad(2) + cosqlon*yqgrad(2))/J

!RETURN QUASI-DIPOLE BASE VECTORS
f1(1) = sngl(qlatgrad(2))
f1(2) = sngl(-qlatgrad(1))
f2(1) = sngl(-qlongrad(2))
f2(2) = sngl(qlongrad(1))
f1(1) = dble(qlatgrad(2))
f1(2) = dble(-qlatgrad(1))
f2(1) = dble(-qlongrad(2))
f2(2) = dble(qlongrad(1))
f = f1(1)*f2(2) - f1(2)*f2(1)

endif
Expand All @@ -420,12 +420,12 @@ subroutine apxg2all(glat, glon, alt, hr, vecflagin, &

implicit none

real(4), intent(in) :: glat, glon, alt, hr
real(8), intent(in) :: glat, glon, alt, hr
integer(4), intent(in) :: vecflagin
real(4), intent(out) :: qlatout, qlonout, mlat, mlon
real(4), intent(out) :: f1(1:2), f2(1:2), f
real(4), intent(out) :: d1(1:3), d2(1:3), d3(1:3), d
real(4), intent(out) :: e1(1:3), e2(1:3), e3(1:3)
real(8), intent(out) :: qlatout, qlonout, mlat, mlon
real(8), intent(out) :: f1(1:2), f2(1:2), f
real(8), intent(out) :: d1(1:3), d2(1:3), d3(1:3), d
real(8), intent(out) :: e1(1:3), e2(1:3), e3(1:3)

integer(4) :: i
real(8) :: cosmlat, Rrat, denom
Expand Down Expand Up @@ -455,7 +455,7 @@ subroutine apxg2all(glat, glon, alt, hr, vecflagin, &
Rrat = (Re + dble(hr))/Reph
cosmlat = cosqlat*dsqrt(Rrat)
if (cosmlat .le. 1D0) then
mlat = sngl(dacos(cosmlat)/dtor)
mlat = dble(dacos(cosmlat)/dtor)
if (qlat .lt. 0D0) mlat = -mlat
endif

Expand All @@ -474,10 +474,10 @@ subroutine apxg2all(glat, glon, alt, hr, vecflagin, &
if (denom .le. 0) return
denom = dsqrt(denom)
do i = 1, 3
d1(i) = sngl(Rrat*dsqrt(Rrat)*qlongrad(i))
d2(i) = sngl(-2D0*Rrat*sinqlat*qlatgrad(i)/denom)
d1(i) = dble(Rrat*dsqrt(Rrat)*qlongrad(i))
d2(i) = dble(-2D0*Rrat*sinqlat*qlatgrad(i)/denom)
enddo
d2(3) = d2(3) - sngl(Rrat*cosqlat/denom)
d2(3) = d2(3) - dble(Rrat*cosqlat/denom)
e3(1) = d1(2)*d2(3) - d1(3)*d2(2)
e3(2) = d1(3)*d2(1) - d1(1)*d2(3)
e3(3) = d1(1)*d2(2) - d1(2)*d2(1)
Expand Down Expand Up @@ -507,12 +507,12 @@ subroutine apxq2g(qlat0, qlon0, alt, prec, glatout, glonout, error)

implicit none

real(4), intent(in) :: qlat0, qlon0, alt, prec
real(4), intent(out) :: glatout, glonout, error
real(8), intent(in) :: qlat0, qlon0, alt, prec
real(8), intent(out) :: glatout, glonout, error

integer(4) :: l, iterm, itermsh, vecflagin, niter
real(4) :: qlatout, qlonout, errorlast
real(4) :: mlon, f1(1:2), f2(1:2), f
real(8) :: qlatout, qlonout, errorlast
real(8) :: mlon, f1(1:2), f2(1:2), f
real(8) :: theta, phi
real(8) :: sinqlon0, cosqlon0, sinqlat0, cosqlat0, cotqlat0, zfact
real(8) :: glat, glon
Expand Down Expand Up @@ -558,8 +558,8 @@ subroutine apxq2g(qlat0, qlon0, alt, prec, glatout, glonout, error)
zg = dot_product(sh, zgcoeff)
glon = datan2(yg, xg)
glat = datan2(zg, dsqrt(xg*xg + yg*yg))
glatout = sngl(glat/dtor)
glonout = sngl(glon/dtor)
glatout = dble(glat/dtor)
glonout = dble(glon/dtor)
error = missing

!COMPUTE REFINED GEODETIC COORDINATES
Expand All @@ -576,27 +576,27 @@ subroutine apxq2g(qlat0, qlon0, alt, prec, glatout, glonout, error)
sinqlat0 = dcos(theta)
coserror = sinqlat0*sinqlat + cosqlat0*cosqlat*cos(phi - qlon)
if (coserror .gt. 1) coserror = 1D0
error = sngl(acos(coserror)/dtor)
error = dble(acos(coserror)/dtor)
niter = 0
errorlast = 9999.0
errorlast = 9999.0d0
!OUTSIDE OF QD POLES
if ((abs(qlat0) .lt. 88.0) .and. (error .gt. prec)) then
if ((abs(qlat0) .lt. 88.0d0) .and. (error .gt. prec)) then
vecflagin = 0
do while ((error .gt. prec) .and. (niter .lt. 10) .and. (error .lt. 1.3*errorlast))
do while ((error .gt. prec) .and. (niter .lt. 10) .and. (error .lt. 1.3d0*errorlast))
delqlon = phi - qlon
if (abs(delqlon) .gt. pi) delqlon = -sign(twopi - abs(delqlon), delqlon)
delqlon = cosqlat0*delqlon
delqlat = pid2 - theta - qlat
xg = xg + xggrad(1)*delqlon + xggrad(2)*delqlat
yg = yg + yggrad(1)*delqlon + yggrad(2)*delqlat
zg = zg + zggrad(1)*delqlon + zggrad(2)*delqlat
glonout = sngl(datan2(yg, xg)/dtor)
glatout = sngl(datan2(zg, dsqrt(xg*xg + yg*yg))/dtor)
glonout = dble(datan2(yg, xg)/dtor)
glatout = dble(datan2(zg, dsqrt(xg*xg + yg*yg))/dtor)
call apxg2q(glatout, glonout, alt, vecflagin, qlatout, qlonout, f1, f2, f)
coserror = sinqlat0*sinqlat + cosqlat0*cosqlat*cos(phi - qlon)
if (coserror .gt. 1) coserror = 1D0
errorlast = error
error = sngl(acos(coserror)/dtor)
error = dble(acos(coserror)/dtor)
niter = niter + 1
enddo
!NEAR QD POLES
Expand All @@ -606,20 +606,20 @@ subroutine apxq2g(qlat0, qlon0, alt, prec, glatout, glonout, error)
cotqlat0 = cosqlat0/sinqlat0
cosglat = dcos(glat)
vecflagin = 1
do while ((error .gt. prec) .and. (niter .lt. 10) .and. (error .lt. 1.3*errorlast))
do while ((error .gt. prec) .and. (niter .lt. 10) .and. (error .lt. 1.3d0*errorlast))
zfact = zq*cotqlat0
delxq = zfact*cosqlon0 - xq
delyq = zfact*sinqlon0 - yq
denom = xqgrad(1)*yqgrad(2) - xqgrad(2)*yqgrad(1)
glon = glon + (delxq*yqgrad(2) - delyq*xqgrad(2))/denom/cosglat
glat = glat + (delyq*xqgrad(1) - delxq*yqgrad(1))/denom
glonout = sngl(glon/dtor)
glatout = sngl(glat/dtor)
glonout = dble(glon/dtor)
glatout = dble(glat/dtor)
call apxg2q(glatout, glonout, alt, vecflagin, qlatout, qlonout, f1, f2, f)
coserror = sinqlat0*sinqlat + cosqlat0*cosqlat*cos(phi - qlon)
if (coserror .gt. 1) coserror = 1D0
errorlast = error
error = sngl(acos(coserror)/dtor)
error = dble(acos(coserror)/dtor)
niter = niter + 1
enddo
endif
Expand Down Expand Up @@ -780,4 +780,3 @@ subroutine alfbasis(nmax, mmax, theta, P, V, W)
end subroutine alfbasis

!***************************************************************************************************

10 changes: 5 additions & 5 deletions src/calc_electrodynamics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,10 +198,10 @@ subroutine UA_calc_electrodynamics(UAi_nMLTs, UAi_nLats)
call stop_gitm("Error allocating array DivJuAltMC")
endif

date = iStartTime(1) + float(iJulianDay)/float(jday(iStartTime(1), 12, 31))
date = dble(iStartTime(1)) + dble(iJulianDay)/dble(jday(iStartTime(1), 12, 31))

iStart = float(iProc)/nProcs*(nMagLons + 1) + 1
iEnd = float(iProc + 1)/nProcs*(nMagLons + 1)
iStart = dble(iProc)/dble(nProcs)*(nMagLons + 1) + 1
iEnd = dble(iProc + 1)/dble(nProcs)*(nMagLons + 1)

GeoLatMC = -1.0e32
GeoLonMC = -1.0e32
Expand All @@ -214,9 +214,9 @@ subroutine UA_calc_electrodynamics(UAi_nMLTs, UAi_nLats)
write(*, *) "==> Calculating Apex->Geo", i, iStart, iEnd
do j = 1, nMagLats

MagLatMC(i, j) = float(j - 1)*MagLatRes - DynamoHighLatBoundary
MagLatMC(i, j) = dble(j - 1)*MagLatRes - DynamoHighLatBoundary

MagLonMC(i, j) = 360.0*float(i - 1)/float(nMagLons)
MagLonMC(i, j) = 360.0d0*dble(i - 1)/dble(nMagLons)

if (UseApex) then

Expand Down
Loading
Loading