diff --git a/src/Mars.f90 b/src/Mars.f90 index a19d10b7..215d26a2 100644 --- a/src/Mars.f90 +++ b/src/Mars.f90 @@ -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 diff --git a/src/ModHwm14.f90 b/src/ModHwm14.f90 index fd3837b5..66f13f78 100644 --- a/src/ModHwm14.f90 +++ b/src/ModHwm14.f90 @@ -68,10 +68,10 @@ subroutine hwm14(iyd, sec, alt, glat, glon, stl, f107a, f107, ap, path, w) use hwm implicit none integer(4), intent(in) :: iyd - real(4), intent(in) :: sec, alt, glat, glon, stl, f107a, f107 - real(4), intent(in) :: ap(2) - real(4), intent(out) :: w(2) - real(4) :: dw(2) + real(8), intent(in) :: sec, alt, glat, glon, stl, f107a, f107 + real(8), intent(in) :: ap(2) + real(8), intent(out) :: w(2) + real(8) :: dw(2) character(*), intent(in) :: path pathdefault = path @@ -263,14 +263,14 @@ module dwm integer(4) :: nvshterm ! # of VSH basis functions integer(4), allocatable :: termarr(:, :) ! 3 x nterm index of coupled terms - real(4), allocatable :: coeff(:) ! Model coefficients - real(4), allocatable :: vshterms(:, :) ! VSH basis values - real(4), allocatable :: termval(:, :) ! Term values to which coefficients are applied + real(8), allocatable :: coeff(:) ! Model coefficients + real(8), allocatable :: vshterms(:, :) ! VSH basis values + real(8), allocatable :: termval(:, :) ! Term values to which coefficients are applied real(8), allocatable :: dpbar(:, :) ! Associated lengendre fns real(8), allocatable :: dvbar(:, :) real(8), allocatable :: dwbar(:, :) real(8), allocatable :: mltterms(:, :) ! MLT Fourier terms - real(4) :: twidth ! Transition width of high-lat mask + real(8) :: twidth ! Transition width of high-lat mask real(8), parameter :: pi = 3.1415926535897932 real(8), parameter :: dtor = pi/180.d0 @@ -512,9 +512,9 @@ subroutine hwmqt(IYD, SEC, ALT, GLAT, GLON, STL, F107A, F107, AP, W) implicit none integer, intent(in) :: IYD - real(4), intent(in) :: SEC, ALT, GLAT, GLON, STL, F107A, F107 - real(4), intent(in) :: AP(2) - real(4), intent(out) :: W(2) + real(8), intent(in) :: SEC, ALT, GLAT, GLON, STL, F107A, F107 + real(8), intent(in) :: AP(2) + real(8), intent(out) :: W(2) ! Local variables @@ -760,8 +760,8 @@ subroutine hwmqt(IYD, SEC, ALT, GLAT, GLON, STL, F107A, F107, AP, W) enddo - w(1) = sngl(v) - w(2) = sngl(u) + w(1) = dble(v) + w(2) = dble(u) return @@ -923,6 +923,8 @@ subroutine initdwm(nmaxout, mmaxout) implicit none integer(4), intent(out) :: nmaxout, mmaxout + real(4), allocatable :: coeff4(:) + real(4) :: twidth4 call findandopen(dwmdefault, 23) if (allocated(termarr)) deallocate(termarr, coeff) @@ -930,8 +932,12 @@ subroutine initdwm(nmaxout, mmaxout) allocate(termarr(0:2, 0:nterm - 1)) read(23) termarr allocate(coeff(0:nterm - 1)) - read(23) coeff - read(23) twidth + allocate(coeff4(0:nterm - 1)) + read(23) coeff4 + coeff = dble(coeff4) + read(23) twidth4 + twidth = dble(twidth4) + deallocate(coeff4) close(23) if (allocated(termval)) deallocate(termval, dpbar, dvbar, dwbar, mltterms, vshterms) @@ -960,18 +966,18 @@ subroutine dwm07(IYD, SEC, ALT, GLAT, GLON, AP, DW) implicit none INTEGER, intent(in) :: IYD - REAL(4), intent(in) :: SEC, ALT, GLAT, GLON - REAL(4), intent(in) :: AP(2) - REAL(4), intent(out) :: DW(2) + REAL(8), intent(in) :: SEC, ALT, GLAT, GLON + REAL(8), intent(in) :: AP(2) + REAL(8), intent(out) :: DW(2) - real(4), save :: day, ut, mlat, mlon, mlt, kp - real(4) :: mmpwind, mzpwind - real(4), save :: f1e, f1n, f2e, f2n - real(4), save :: glatlast = 1.0e16, glonlast = 1.0e16 - real(4), save :: daylast = 1.0e16, utlast = 1.0e16, aplast = 1.0e16 - real(4), parameter :: talt = 125.0 !, twidth=5.0 + real(8), save :: day, ut, mlat, mlon, mlt, kp + real(8) :: mmpwind, mzpwind + real(8), save :: f1e, f1n, f2e, f2n + real(8), save :: glatlast = 1.0d16, glonlast = 1.0d16 + real(8), save :: daylast = 1.0d16, utlast = 1.0d16, aplast = 1.0d16 + real(8), parameter :: talt = 125.0d0 !, twidth=5.0 - real(4), external :: ap2kp, mltcalc + real(8), external :: ap2kp, mltcalc !CONVERT AP TO KP if (ap(2) .ne. aplast) then @@ -984,8 +990,8 @@ subroutine dwm07(IYD, SEC, ALT, GLAT, GLON, AP, DW) endif !COMPUTE QD MAGNETIC LOCAL TIME (LOW-PRECISION) - day = real(mod(iyd, 1000)) - ut = sec/3600.0 + day = dble(mod(iyd, 1000)) + ut = sec/3600.0d0 if ((day .ne. daylast) .or. (ut .ne. utlast) .or. & (glat .ne. glatlast) .or. (glon .ne. glonlast)) then mlt = mltcalc(mlat, mlon, day, ut) @@ -999,7 +1005,7 @@ subroutine dwm07(IYD, SEC, ALT, GLAT, GLON, AP, DW) dw(2) = f2e*mmpwind + f1e*mzpwind !APPLY HEIGHT PROFILE - dw = dw/(1 + exp(-(alt - talt)/twidth)) + dw = dw/(1.0d0 + exp(-(alt - talt)/twidth)) glatlast = glat glonlast = glon @@ -1018,22 +1024,22 @@ subroutine dwm07b(mlt, mlat, kp, mmpwind, mzpwind) use alf, only: alfbasis implicit none - real(4), intent(in) :: mlt !Magnetic local time (hours) - real(4), intent(in) :: mlat !Magnetic latitude (degrees) - real(4), intent(in) :: kp !3-hour Kp + real(8), intent(in) :: mlt !Magnetic local time (hours) + real(8), intent(in) :: mlat !Magnetic latitude (degrees) + real(8), intent(in) :: kp !3-hour Kp - real(4), intent(out) :: mmpwind !Mer. disturbance wind (+north, QD coordinates) - real(4), intent(out) :: mzpwind !Zon. disturbance wind (+east, QD coordinates) + real(8), intent(out) :: mmpwind !Mer. disturbance wind (+north, QD coordinates) + real(8), intent(out) :: mzpwind !Zon. disturbance wind (+east, QD coordinates) ! Local variables integer(4) :: iterm, ivshterm, n, m - real(4) :: termvaltemp(0:1) - real(4), save :: kpterms(0:2) - real(4) :: latwgtterm - real(4), save :: mltlast = 1.e16, mlatlast = 1.e16, kplast = 1.e16 + real(8) :: termvaltemp(0:1) + real(8), save :: kpterms(0:2) + real(8) :: latwgtterm + real(8), save :: mltlast = 1.d16, mlatlast = 1.d16, kplast = 1.d16 real(8) :: theta, phi, mphi - real(4), external :: latwgt2 + real(8), external :: latwgt2 !LOAD MODEL PARAMETERS IF NECESSARY if (dwminit) call initdwm(nmaxdwm, mmaxdwm) @@ -1058,17 +1064,17 @@ subroutine dwm07b(mlt, mlat, kp, mmpwind, mzpwind) if ((mlat .ne. mlatlast) .or. (mlt .ne. mltlast)) then ivshterm = 0 do n = 1, nmax - vshterms(0, ivshterm) = -sngl(dvbar(n, 0)*mltterms(0, 0)) - vshterms(0, ivshterm + 1) = sngl(dwbar(n, 0)*mltterms(0, 0)) + vshterms(0, ivshterm) = -dble(dvbar(n, 0)*mltterms(0, 0)) + vshterms(0, ivshterm + 1) = dble(dwbar(n, 0)*mltterms(0, 0)) vshterms(1, ivshterm) = -vshterms(0, ivshterm + 1) vshterms(1, ivshterm + 1) = vshterms(0, ivshterm) ivshterm = ivshterm + 2 do m = 1, mmax if (m .gt. n) cycle - vshterms(0, ivshterm) = -sngl(dvbar(n, m)*mltterms(m, 0)) - vshterms(0, ivshterm + 1) = sngl(dvbar(n, m)*mltterms(m, 1)) - vshterms(0, ivshterm + 2) = sngl(dwbar(n, m)*mltterms(m, 1)) - vshterms(0, ivshterm + 3) = sngl(dwbar(n, m)*mltterms(m, 0)) + vshterms(0, ivshterm) = -dble(dvbar(n, m)*mltterms(m, 0)) + vshterms(0, ivshterm + 1) = dble(dvbar(n, m)*mltterms(m, 1)) + vshterms(0, ivshterm + 2) = dble(dwbar(n, m)*mltterms(m, 1)) + vshterms(0, ivshterm + 3) = dble(dwbar(n, m)*mltterms(m, 0)) vshterms(1, ivshterm) = -vshterms(0, ivshterm + 2) vshterms(1, ivshterm + 1) = -vshterms(0, ivshterm + 3) vshterms(1, ivshterm + 2) = vshterms(0, ivshterm) @@ -1088,7 +1094,7 @@ subroutine dwm07b(mlt, mlat, kp, mmpwind, mzpwind) !GENERATE COUPLED TERMS do iterm = 0, nterm - 1 - termvaltemp = (/1.0, 1.0/) + termvaltemp = (/1.0d0, 1.0d0/) if (termarr(0, iterm) .ne. 999) termvaltemp = termvaltemp*vshterms(0:1, termarr(0, iterm)) if (termarr(1, iterm) .ne. 999) termvaltemp = termvaltemp*kpterms(termarr(1, iterm)) if (termarr(2, iterm) .ne. 999) termvaltemp = termvaltemp*latwgtterm @@ -1113,18 +1119,18 @@ end subroutine dwm07b function ap2kp(ap0) - real(4), parameter :: apgrid(0:27) = (/0., 2., 3., 4., 5., 6., 7., 9., 12., 15., 18., & - 22., 27., 32., 39., 48., 56., 67., 80., 94., & - 111., 132., 154., 179., 207., 236., 300., 400./) - real(4), parameter :: kpgrid(0:27) = (/0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., & - 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., & - 22., 23., 24., 25., 26., 27./)/3.0 - real(4) :: ap0, ap, ap2kp + real(8), parameter :: apgrid(0:27) = (/0.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0, 7.d0, 9.d0, 12.d0, 15.d0, 18.d0, & + 22.d0, 27.d0, 32.d0, 39.d0, 48.d0, 56.d0, 67.d0, 80.d0, 94.d0, & + 111.d0, 132.d0, 154.d0, 179.d0, 207.d0, 236.d0, 300.d0, 400.d0/) + real(8), parameter :: kpgrid(0:27) = (/0.d0, 1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0, 7.d0, 8.d0, 9.d0, 10.d0, 11.d0, & + 12.d0, 13.d0, 14.d0, 15.d0, 16.d0, 17.d0, 18.d0, 19.d0, 20.d0, 21.d0, & + 22.d0, 23.d0, 24.d0, 25.d0, 26.d0, 27.d0/)/3.0d0 + real(8) :: ap0, ap, ap2kp integer(4) :: i ap = ap0 - if (ap .lt. 0) ap = 0 - if (ap .gt. 400) ap = 400 + if (ap .lt. 0.0d0) ap = 0.0d0 + if (ap .gt. 400.0d0) ap = 400.0d0 i = 1 do while (ap .gt. apgrid(i)) @@ -1133,7 +1139,7 @@ function ap2kp(ap0) if (ap .eq. apgrid(i)) then ap2kp = kpgrid(i) else - ap2kp = kpgrid(i - 1) + (ap - apgrid(i - 1))/(3.0*(apgrid(i) - apgrid(i - 1))) + ap2kp = kpgrid(i - 1) + (ap - apgrid(i - 1))/(3.0d0*(apgrid(i) - apgrid(i - 1))) endif return @@ -1162,7 +1168,7 @@ module gd2qdc real(8), allocatable :: shgradtheta(:) !Array to hold spherical harmonic gradients real(8), allocatable :: shgradphi(:) !Array to hold spherical harmonic gradients real(8), allocatable :: normadj(:) !Adjustment to VSH normalization factor - real(4) :: epoch, alt + real(8) :: epoch, alt real(8), parameter :: pi = 3.1415926535897932d0 real(8), parameter :: dtor = pi/180.0d0 @@ -1180,11 +1186,14 @@ subroutine initgd2qd() character(250) :: datafile integer(4) :: iterm, n integer(4) :: j + real(4) :: epoch4, alt4 datafile = trim(pathdefault)//'gd2qd.dat' call findandopen(datafile, 23) - read(23) nmax, mmax, nterm, epoch, alt + read(23) nmax, mmax, nterm, epoch4, alt4 + epoch = dble(epoch4) + alt = dble(alt4) if (allocated(coeff)) then deallocate(coeff, xcoeff, ycoeff, zcoeff, sh, shgradtheta, shgradphi, normadj) endif @@ -1229,9 +1238,9 @@ subroutine gd2qd(glatin, glon, qlat, qlon, f1e, f1n, f2e, f2n) implicit none - real(4), intent(in) :: glatin, glon - real(4), intent(out) :: qlat, qlon - real(4), intent(out) :: f1e, f1n, f2e, f2n + real(8), intent(in) :: glatin, glon + real(8), intent(out) :: qlat, qlon + real(8), intent(out) :: f1e, f1n, f2e, f2n integer(4) :: n, m, i real(8) :: glat, theta, phi @@ -1244,13 +1253,13 @@ subroutine gd2qd(glatin, glon, qlat, qlon, f1e, f1n, f2e, f2n) if (gd2qdinit) call initgd2qd() - glat = dble(glatin) + glat = glatin if (glat .ne. glatalf) then theta = (90.d0 - glat)*dtor call alfbasis(nmax, mmax, theta, gpbar, gvbar, gwbar) glatalf = glat endif - phi = dble(glon)*dtor + phi = glon*dtor i = 0 do n = 0, nmax @@ -1283,8 +1292,8 @@ subroutine gd2qd(glatin, glon, qlat, qlon, f1e, f1n, f2e, f2n) sinqlon = dsin(qlonrad) cosqlat = x*cosqlon + y*sinqlon - qlat = sngl(datan2(z, cosqlat)/dtor) - qlon = sngl(qlonrad/dtor) + qlat = dble(datan2(z, cosqlat)/dtor) + qlon = dble(qlonrad/dtor) xgradtheta = dot_product(shgradtheta, xcoeff) ygradtheta = dot_product(shgradtheta, ycoeff) @@ -1294,10 +1303,10 @@ subroutine gd2qd(glatin, glon, qlat, qlon, f1e, f1n, f2e, f2n) ygradphi = dot_product(shgradphi, ycoeff) zgradphi = dot_product(shgradphi, zcoeff) - f1e = sngl(-zgradtheta*cosqlat + (xgradtheta*cosqlon + ygradtheta*sinqlon)*z) - f1n = sngl(-zgradphi*cosqlat + (xgradphi*cosqlon + ygradphi*sinqlon)*z) - f2e = sngl(ygradtheta*cosqlon - xgradtheta*sinqlon) - f2n = sngl(ygradphi*cosqlon - xgradphi*sinqlon) + f1e = dble(-zgradtheta*cosqlat + (xgradtheta*cosqlon + ygradtheta*sinqlon)*z) + f1n = dble(-zgradphi*cosqlat + (xgradphi*cosqlon + ygradphi*sinqlon)*z) + f2e = dble(ygradtheta*cosqlon - xgradtheta*sinqlon) + f2n = dble(ygradphi*cosqlon - xgradphi*sinqlon) return @@ -1315,8 +1324,8 @@ function mltcalc(qlat, qlon, day, ut) implicit none - real(4), intent(in) :: qlat, qlon, day, ut - real(4) :: mltcalc + real(8), intent(in) :: qlat, qlon, day, ut + real(8) :: mltcalc integer(4) :: n, m, i real(8) :: asunglat, asunglon, asunqlon @@ -1329,7 +1338,7 @@ function mltcalc(qlat, qlon, day, ut) if (gd2qdinit) call initgd2qd() !COMPUTE GEOGRAPHIC COORDINATES OF ANTI-SUNWARD DIRECTION (LOW PRECISION) - asunglat = -asin(sin((dble(day) + dble(ut)/24.0d0 - 80.0d0)*dtor)*sineps)/dtor + asunglat = -asin(sin((day + ut/24.0d0 - 80.0d0)*dtor)*sineps)/dtor asunglon = -ut*15.d0 !COMPUTE MAGNETIC COORDINATES OF ANTI-SUNWARD DIRECTION @@ -1353,10 +1362,10 @@ function mltcalc(qlat, qlon, day, ut) enddo x = dot_product(sh, xcoeff) y = dot_product(sh, ycoeff) - asunqlon = sngl(datan2(y, x)/dtor) + asunqlon = dble(datan2(y, x)/dtor) !COMPUTE MLT - mltcalc = (qlon - asunqlon)/15.0 + mltcalc = (qlon - asunqlon)/15.0d0 return @@ -1370,20 +1379,20 @@ subroutine kpspl3(kp, kpterms) implicit none - real(4), intent(in) :: kp - real(4), intent(out) :: kpterms(0:2) + real(8), intent(in) :: kp + real(8), intent(out) :: kpterms(0:2) integer(4) :: i, j - real(4) :: x, kpspl(0:6) - real(4), parameter :: node(0:7) = (/-10., -8., 0., 2., 5., 8., 18., 20./) + real(8) :: x, kpspl(0:6) + real(8), parameter :: node(0:7) = (/-10.d0, -8.d0, 0.d0, 2.d0, 5.d0, 8.d0, 18.d0, 20.d0/) - x = max(kp, 0.0) - x = min(x, 8.0) + x = max(kp, 0.0d0) + x = min(x, 8.0d0) - kpterms(0:2) = 0.0 + kpterms(0:2) = 0.0d0 do i = 0, 6 - kpspl(i) = 0.0 - if ((x .ge. node(i)) .and. (x .lt. node(i + 1))) kpspl(i) = 1.0 + kpspl(i) = 0.0d0 + if ((x .ge. node(i)) .and. (x .lt. node(i + 1))) kpspl(i) = 1.0d0 enddo do j = 2, 3 do i = 0, 8 - j - 1 @@ -1407,24 +1416,24 @@ function latwgt2(mlat, mlt, kp0, twidth) implicit none - real(4) :: latwgt2 - real(4) :: mlat, mlt, kp0, kp, twidth - real(4) :: mltrad, sinmlt, cosmlt, tlat + real(8) :: latwgt2 + real(8) :: mlat, mlt, kp0, kp, twidth + real(8) :: mltrad, sinmlt, cosmlt, tlat - real(4), parameter :: coeff(0:5) = (/65.7633, -4.60256, -3.53915, & - -1.99971, -0.752193, 0.972388/) + real(8), parameter :: coeff(0:5) = (/65.7633d0, -4.60256d0, -3.53915d0, & + -1.99971d0, -0.752193d0, 0.972388d0/) - real(4), parameter :: pi = 3.141592653590 - real(4), parameter :: dtor = pi/180.d0 + real(8), parameter :: pi = 3.141592653590d0 + real(8), parameter :: dtor = pi/180.d0 - mltrad = mlt*15.0*dtor + mltrad = mlt*15.0d0*dtor sinmlt = sin(mltrad) cosmlt = cos(mltrad) - kp = max(kp0, 0.0) - kp = min(kp, 8.0) + kp = max(kp0, 0.0d0) + kp = min(kp, 8.0d0) tlat = coeff(0) + coeff(1)*cosmlt + coeff(2)*sinmlt + & kp*(coeff(3) + coeff(4)*cosmlt + coeff(5)*sinmlt) - latwgt2 = 1.0/(1 + exp(-(abs(mlat) - tlat)/twidth)) + latwgt2 = 1.0d0/(1.0d0 + exp(-(abs(mlat) - tlat)/twidth)) return diff --git a/src/Venus.f90 b/src/Venus.f90 index e5868a5c..ed71d9f7 100644 --- a/src/Venus.f90 +++ b/src/Venus.f90 @@ -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 diff --git a/src/apex_more.f90 b/src/apex_more.f90 index b42c5c33..c8732bba 100644 --- a/src/apex_more.f90 +++ b/src/apex_more.f90 @@ -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 @@ -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) diff --git a/src/apex_routines.f b/src/apex_routines.f index 244e1d79..3c236a56 100644 --- a/src/apex_routines.f +++ b/src/apex_routines.f @@ -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 @@ -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) @@ -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 @@ -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) @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/apexsh.f90 b/src/apexsh.f90 index 3be8f72a..5464adf7 100644 --- a/src/apexsh.f90 +++ b/src/apexsh.f90 @@ -148,7 +148,7 @@ 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) @@ -156,8 +156,8 @@ module apxshmodule ! 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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -576,13 +576,13 @@ 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 @@ -590,13 +590,13 @@ subroutine apxq2g(qlat0, qlon0, alt, prec, glatout, glonout, error) 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 @@ -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 @@ -780,4 +780,3 @@ subroutine alfbasis(nmax, mmax, theta, P, V, W) end subroutine alfbasis !*************************************************************************************************** - diff --git a/src/calc_electrodynamics.f90 b/src/calc_electrodynamics.f90 index 30b98bf3..3e6ea63d 100644 --- a/src/calc_electrodynamics.f90 +++ b/src/calc_electrodynamics.f90 @@ -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 @@ -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 diff --git a/src/init_b0.f90 b/src/init_b0.f90 index 332acf2d..e8c7a455 100644 --- a/src/init_b0.f90 +++ b/src/init_b0.f90 @@ -28,7 +28,7 @@ subroutine init_b0 AltMinIono = (2*RadialDistance_GB(1, 1, -1, 1) - & RadialDistance_GB(1, 1, 1, 1) - RBody)/1000.0 - date = iStartTime(1) + float(iJulianDay)/float(jday(iStartTime(1), 12, 31)) + date = dble(iStartTime(1)) + dble(iJulianDay)/dble(jday(iStartTime(1), 12, 31)) call loadapxsh(apexfile, date) do iBlock = 1, nBlocks @@ -179,7 +179,7 @@ subroutine get_magfield_all(GeoLat, GeoLon, GeoAlt, alat, alon, xmag, ymag, zmag real :: alatp, alatm, alonp, alonm, sinIm, londiff integer, external :: jday !-------------------------------------------------------------------------- - date = iStartTime(1) + float(iJulianDay)/float(jday(iStartTime(1), 12, 31)) + date = dble(iStartTime(1)) + dble(iJulianDay)/dble(jday(iStartTime(1), 12, 31)) twodegrees = 2.0*pi/180.0 @@ -596,4 +596,3 @@ subroutine test_mag_point(rBelow, LShell, RBody) endif end subroutine test_mag_point - diff --git a/src/init_msis.Earth.f90 b/src/init_msis.Earth.f90 index ee59106c..ac998fb9 100644 --- a/src/init_msis.Earth.f90 +++ b/src/init_msis.Earth.f90 @@ -276,8 +276,8 @@ subroutine init_msis real :: ffactor, h, no real, dimension(7) :: ap = 10.0 - real*4 :: hwm_utime, hwm_alt, hwm_lat, hwm_lon, hwm_lst - real*4 :: hwm_f107a, hwm_f107, hwm_ap(2), qw(2) + real(8) :: hwm_utime, hwm_alt, hwm_lat, hwm_lon, hwm_lst + real(8) :: hwm_f107a, hwm_f107, hwm_ap(2), qw(2) character(250) :: path = './DataIn/LowerBCs/' @@ -534,8 +534,8 @@ subroutine msis_bcs(iJulianDay, UTime, Alt, LatIn, LonIn, Lst, & real :: AP_I(7), ffactor, no integer :: iyd, iJulianDayOffset - real*4 :: hwm_utime, hwm_alt, hwm_lat, hwm_lon, hwm_lst - real*4 :: hwm_f107a, hwm_f107, hwm_ap(2), qw(2) + real(8) :: hwm_utime, hwm_alt, hwm_lat, hwm_lon, hwm_lst + real(8) :: hwm_f107a, hwm_f107, hwm_ap(2), qw(2) character(250) :: path = './DataIn/LowerBCs/' diff --git a/src/tides.f90 b/src/tides.f90 index 764f5d7d..b2a9df3f 100644 --- a/src/tides.f90 +++ b/src/tides.f90 @@ -685,7 +685,7 @@ subroutine update_tides iLonHalf = pi/dLonGswm - rfac = float(iTimeArray(5))/60.0 + float(iTimeArray(6))/3600.0 + rfac = dble(iTimeArray(5))/60.0d0 + dble(iTimeArray(6))/3600.0d0 iFac1 = iTimeArray(4) + 1 ! this will go from 1-24 iFac2 = mod(iTimeArray(4) + 1, 24) + 1 ! this will go from 2-25 -> 2-24 back to 1