MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64!modern DOUBLE PRECISION
  INTEGER :: i,j,k,n,m,p,h,s,lines
  REAL(dp) :: PI,y,month,day,time,monitor,lat,lon,alt,max,d,de,dose,c,monitors,maxs,lats,lons,alts,zs,VCO,ds,des,days,times,altss, &
  SALT,P0,P1,Q0,Q1,R0,R1,SHI,SHO,SMR,SZ,t,WHI,WLO,WTR,WX,WY,WZ,lonss,latss,z2s,dayss,cs,VCOs,z3s,PNEU,CION,CSS,dss,dess,arg,arg1,  &
  arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10,arg11,arg12,zsff,SL,CL,SF,CF,PHIF,PHIM,XM,csf,VCOsf,zssf,PIRM,PIRX,ratio1,ratio2,  &
  SLOPE
  REAL(dp), DIMENSION(20) :: ZNW,ZNS,ZSW,ZSS
  REAL(dp), DIMENSION(36,17,20) :: WJAN,WJUL
  REAL(dp), PARAMETER :: RO36 = 0.0277778
  CHARACTER(LEN=1) :: clock
END MODULE Shared

PROGRAM AIR!equations and page numbers refer to J.W. Wilson et al., NASA Reference Publication 1257, "Transport Methods and Interactions for Space Radiations," December 1991.
  USE Shared, ONLY : i,j,lines,h,m,s,PI,y,month,day,time,monitor,lat,lon,alt,max,d,de,dose,clock
  IMPLICIT NONE
  PI = ACOS(-1.0)!PI to machine precision
  y = 2005.0!input
  month = 6.0!input
  month = month - 1.0!fractional year
  day = 8.0!input
  day = day - 1.0!fractional month
  clock = 'a'!input
  OPEN(1, file='DSTB.txt', status='old')!copy/paste from spreadsheet into text file with constant decimal places as formatted in Excel
  OPEN(2, file='DSTB_AIR.txt', status='unknown')
  OPEN(3, file='DSTB_PNEU-CION.txt', status='unknown')
  I = 0!Classical FORTRAN, pg. 203
  111 I = I + 1
  READ(1,'(i2)',END=11) h
  GO TO 111
  11 REWIND(1)!Classical FORTRAN
  lines = I - 1!always over-counting
  dose = 0.0
  main: DO j = 1, lines
    READ(1,1111) h,m,s,lat,lon,alt!lat = 75.0 for Devon Island; lon = -88.0 for Devon Island; h = h - 6 converts GPS time to local time; altitude in feet
    h = h - 1!fractional day
    IF(clock == 'p') h = h + 12!convert to 24-hour clock
    m = m - 1!fractional hour
    IF(s >= 30) m = m + 1!round to the nearest minute
    alt = alt * 0.3048!convert from feet to meters
    time = y + month / 12.0 + day / 365.0 + FLOAT(h) / (365.0 * 24.0) + FLOAT(m) / (365.0 * 24.0 * 60.0)
    monitor = 3932.6038 + 360.1950 * SIN((2 * PI * time) / 11.6923 + 3.1656)!Climax, CO Neutron Monitor (periodic fit, daily averaged)
    max = 4292.7988!3932.6038 + 360.1950 * 1.0
    CALL ALDOS(monitor,max,day,lat,lon,alt,time,d,de)
    WRITE(2,'(f9.4,e15.6)') d * 1000.0,alt / 0.3048!convert back to feet
    dose = dose + (d * 1000.0) * (1.0 / 60.0) * 0.01!convert from mrad/hr to urad/hr; DSTB data logged about every minute; convert from urad to uGy
    PRINT *, ' '
    PRINT '(a10,i7,a3)', ' Altitude:',NINT(alt / 0.3048),' ft' !round
  END DO main
  !1111 FORMAT(i2,1x,i2,1x,i2,1x,f11.8,1x,f12.8,1x,e12.6)
  1111 FORMAT(i2,1x,i2,1x,i2,1x,f11.8,1x,f13.8,1x,e12.6)!DSTB; a tab in the data also looks like a single space (1X)
  PRINT *, ' '
  PRINT *, '-------------------'
  PRINT *, ' '
  PRINT '(a6,f5.1,a4)', ' Dose:',dose,' uGy'
  CLOSE(3)
  CLOSE(2)
  CLOSE(1)
END PROGRAM AIR

SUBROUTINE ALDOS(monitors,maxs,days,lats,lons,alts,times,ds,des)
  USE Shared, ONLY : dp,c,VCO,zs
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: monitors,maxs,days,lats,lons,times
  REAL(dp), INTENT(OUT) :: ds,des
  REAL(dp), INTENT(INOUT) :: alts!INTENT(IN), but also altered
  CALL PALT(days,lats,lons,alts,zs)!get interpolated pressure (mbar)
  c = (monitors / maxs) * 100.0
  VCO = -0.00336667 * times + 9.63472!Proc 27th Int Cosmic Ray Conf, 4063 (2001)
  !The first input variable is the CNM count rate in % of maximum.
  !The second input variable is the vertical cutoff in GV.
  !The third input variable is now the atmospheric depth in g/cm2 because pressure in mbar and depth in g/cm2 are the same thing.
  CALL HDOSE(c,VCO,zs,ds,des)!get dose and dose equivalent
END

!This routine performs tri-variate interpolation for pressure in time, latitude, longitude, and altitude for summer and winter pressure fields.
!It also performs sinusoidal time interpolation for a specific season.
!So: gives atmospheric pressure for an altitude.
!Data files contain pressure values for global grid latitudes and longitudes (-90 to +90 lat.; 0 to 360 long.).
!'LGJAN.dat' & 'LGJUL.dat' contain winter and summer solstice pressure fields (natural logarithm of pressure in millibars) for the aforementioned global grid.
!These files also contain altitudes ranging from pressure fields extracted from the Langley Research Center (LaRC) three-dimensional general circulation model (GCM).
SUBROUTINE PALT(dayss,latss,lonss,altss,z2s)
  USE Shared, ONLY : dp,i,j,k,m,n,p,PI,RO36,SALT,P0,P1,Q0,Q1,R0,R1,SHI,SHO,SMR,SZ,t,WHI,WJAN,WJUL,WLO,WTR,WX,WY,WZ,ZNS,ZNW,ZSS,ZSW
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: dayss,latss,lonss
  REAL(dp), INTENT(OUT) :: z2s
  REAL(dp), INTENT(INOUT) :: altss!INTENT(IN), but also altered
  REAL(dp) :: TEXP
  EXTERNAL TEXP
  OPEN(12, file='LGJAN.dat', status='old')
  READ(12,'(6f12.5)') WJAN!implied loop
  CLOSE(12)
  OPEN(13, file='LGJUL.dat', status='old')
  READ(13,'(6f12.5)') WJUL!implied loop
  CLOSE(13)
  SALT = -1000.0!for case above 2,000 meters
  IF(altss < 2000.0)THEN
    SALT = altss!Save ALTitude; for case below 2,000 meters
    altss = 2000.0!evaluate at 2,000 meters, then make a correction (below)
  END IF
  t = dayss / 365.0![0.0,1.0]
  WTR = COS(PI * t)
  SMR = SIN(PI * t)
  WX = 0.1 * lonss
  WY = 0.1 * (latss + 80.0)
  WZ = 20.0 - 0.0005 * altss + 1.0E-08
  i = INT(WX) + 1
  j = INT(WY) + 1
  k = INT(WZ) + 1
  P1 = WX - 1.0 * FLOAT(i - 1)
  Q1 = WY - 1.0 * FLOAT(j - 1)
  R1 = WZ - 1.0 * FLOAT(k - 1)
  m = i + 1
  IF(i == 36) m = 1
  n = j + 1
  p = k + 1
  P0 = 1.0 - P1
  Q0 = 1.0 - Q1
  R0 = 1.0 - R1
  IF(ABS(latss) > 80.0) THEN!you're at a pole; which one?
    IF(latss > 80.0) THEN!you're at the North Pole
      j = 17
      Q0 = 0.1 * (latss - 80.0)
      Q1 = 1.0 - Q0
      lat1: DO k = 1, 20
        ZNW(k) = 0.0!get averages
        ZNS(k) = 0.0
        long1: DO i = 1, 36
          ZNW(k) = ZNW(k) + WJAN(i,17,k)
          ZNS(k) = ZNS(k) + WJUL(i,17,k)
        END DO long1
        ZNW(k) = ZNW(k) * RO36
        ZNS(k) = ZNS(k) * RO36
      END DO lat1
      WLO = ZNW(k)
      WHI = ZNW(p)
      SHO = ZNS(k)
      SHI = ZNS(p)
    ELSE!you're at the South Pole
      j = 1
      Q1 = 0.1 * (90.0 + latss)
      Q0 = 1.0 - Q1
      lat2: DO k = 1, 20
        ZSW(k) = 0.0!get averages
        ZSS(k) = 0.0
        long2: DO i = 1, 36
          ZSW(k) = ZSW(k) + WJAN(i,1,k)
          ZSS(k) = ZSS(k) + WJUL(i,1,k)
        END DO long2
        ZSW(k) = ZSW(k) * RO36
        ZSS(k) = ZSS(k) * RO36
      END DO lat2
      WLO = ZSW(k)
      WHI = ZSW(p)
      SHO = ZSS(k)
      SHI = ZSS(p)
    END IF
    WZ = P0 * Q0 * R0 * WJAN(i,j,k) + P1 * Q0 * R0 * WJAN(m,j,k) + P0 * Q0 * R1 * WJAN(i,j,p) + P1 * Q0 * R1 * WJAN(m,j,p) + (P0 * Q1 * R0 + P1 * Q1 * R0) * WLO +   &
    (P0 * Q1 * R1 + P1 * Q1 * R1) * WHI
    SZ = P0 * Q0 * R0 * WJUL(i,j,k) + P1 * Q0 * R0 * WJUL(m,j,k) + P0 * Q0 * R1 * WJUL(i,j,p) + P1 * Q0 * R1 * WJUL(m,j,p) + (P0 * Q1 * R0 + P1 * Q1 * R0) * SHO +   &
    (P0 * Q1 * R1 + P1 * Q1 * R1) * SHI
  ELSE!you're not at a pole
    WZ = P0 * Q0 * R0 * WJAN(i,j,k) + P1 * Q0 * R0 * WJAN(m,j,k) + P0 * Q1 * R0 * WJAN(i,n,k) + P1 * Q1 * R0 * WJAN(m,n,k) + P0 * Q0 * R1 * WJAN(i,j,p) +        &
    P1 * Q0 * R1 * WJAN(m,j,p) + P0 * Q1 * R1 * WJAN(i,n,p) + P1 * Q1 * R1 * WJAN(m,n,p)
    SZ = P0 * Q0 * R0 * WJUL(i,j,k) + P1 * Q0 * R0 * WJUL(m,j,k) + P0 * Q1 * R0 * WJUL(i,n,k) + P1 * Q1 * R0 * WJUL(m,n,k) + P0 * Q0 * R1 * WJUL(i,j,p) +        &
    P1 * Q0 * R1 * WJUL(m,j,p) + P0 * Q1 * R1 * WJUL(i,n,p) + P1 * Q1 * R1 * WJUL(m,n,p)
  END IF
  z2s = 1.02 * SQRT((TEXP(WZ) * WTR)**2 + (TEXP(SZ) * SMR)**2)!sinusoidal variation between solstices; evaluate at 2,000 meters
  IF(INT(SALT) /= INT(-1000.0)) THEN!for case below 2,000 meters
    z2s = 1050.0 * EXP(0.0005 * LOG(z2s / 1050.0) * SALT)!implement correction to account for water vapor abundance
    altss = SALT!restore original value and RETURN
  END IF
END

SUBROUTINE HDOSE(cs,VCOs,z3s,dss,dess)
  USE Shared, ONLY : dp,PNEU,arg1,arg2,CION,CSS
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: cs,VCOs,z3s
  REAL(dp), INTENT(OUT) :: dss,dess
  REAL(dp) :: PHIN,AIRIP,TEXP
  EXTERNAL PHIN,AIRIP,TEXP
  PNEU = PHIN(cs,VCOs,z3s)!neutrons/cm2-sec
  CION = AIRIP(cs,VCOs,z3s)!ion pairs/cm3-sec
  WRITE(3,'(2f7.2)') PNEU,CION
  dss = 0.00172 * CION + 0.05 * PNEU!0.00172 mrad-cm3-sec/hr; 0.05 mrad-cm2-sec/hr, from pg. 532; D rate (mrad/hr)
  arg1 = -z3s / 416.0
  arg2 = -z3s / 65.0
  CSS = 1.0 + 0.35 * TEXP(arg1) - 0.194 * TEXP(arg2)!eq. (13.9), pg. 535
  dess = CSS * 0.00172 * CION + (1.0 + 0.6) * 0.314 * PNEU!0.00172 mrem-cm3-sec/hr; 0.314 mrem-cm2-sec/hr, from pg. 532; H rate (mrem/hr)
END

REAL(dp) FUNCTION PHIN(csf,VCOsf,zsff)!This function embodies the neutron environment model by generating the 1-10 MeV ("fast") neutron flux in units of neutrons/cm2-sec.
  USE Shared, ONLY : dp,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10,arg11,arg12,SL,CL,SF,CF,PHIF,PHIM,XM
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: csf,VCOsf,zsff
  REAL(dp) :: TEXP
  EXTERNAL TEXP
  arg1 = (-VCOsf * VCOsf) / 81.0
  arg2 = (csf - 100.0) / 3.73
  arg3 = (-VCOsf * VCOsf) / 12.96
  PHIM = 0.23 + (1.1 + 0.0167 * (csf - 100.0)) * TEXP(arg1) + (0.991 + 0.051 * (csf - 100.0) + 0.4 * TEXP(arg2)) * TEXP(arg3)!"PHIM" is the fast (1-10 MeV) neutron flux at the transition Maximum altitude; eq. (13.1), pg. 530
  arg4 = (-VCOsf * VCOsf) / 25.0
  arg5 = (csf - 100.0) / 3.73
  arg6 = (-VCOsf * VCOsf) / 139.2
  PHIF = 0.17 + (0.787 + 0.035 * (csf - 100.0)) * TEXP(arg4) + (-0.107 - 0.0265 * (csf - 100.0) + 0.612 * TEXP(arg5)) * TEXP(arg6)!"PHIF" is the fast neutron Flux at a depth of 250 g/cm2; eq.(13.2), pg. 532
  SL = 165.0 + 2.0 * VCOsf!Small Lambda (SL): eq. (13.3), pg. 532; at depths < 250 g/cm2, neutrons attenuate according to this (units: g/cm2)
  arg7 = -2.0 * (csf - 100.0)
  XM = 50.0 + LOG(2000.0 + TEXP(arg7))!transition Maximum altitude (XM): eq. (13.8), pg. 532
  arg8 = 250.0 / SL
  SF = TEXP(arg8) * PHIF!Small F (SF): eq. (13.5), pg. 532
  arg9 = XM / SL
  CL = SL * (1.0 - PHIM * TEXP(arg9) / SF)!Capital Lambda (CL): eq.(13.6), pg.532
  arg10 = (XM / CL) - (XM / SL)
  CF = CL * SF * TEXP(arg10) / SL!Capital F (CF): eq. (13.7), pg. 532
  arg11 = -zsff / SL
  arg12 = -zsff / CL
  PHIN = SF * TEXP(arg11) - CF * TEXP(arg12)!The neutron flux at all altitudes is approximated as "PHIN"; eq. (13.4), pg. 532
END

REAL(dp) FUNCTION AIRIP(csf,VCOsf,zssf)!Neher air ionization data; also: GREP data
  USE Shared, ONLY : dp,i,j,n,m,PIRM,PIRX,ratio1,ratio2,SLOPE,arg
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: csf,VCOsf,zssf
  REAL(dp), DIMENSION(12) :: RT!"RT" in GV
  REAL(dp), DIMENSION(14) :: XT,PIR!"XT" in g/cm2
  DATA XT /30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0,120.0,140.0,200.0,245.0,300.0,700.0/!"700.0" should not be "1034.0" like in the tables--reason unknown
  DATA RT /0.0,0.01,0.16,0.49,1.31,1.97,2.56,5.17,8.44,11.7,14.11,17.0/!"1.31" not in original tables; added later
  REAL(dp), DIMENSION(12,14) :: PIMIN,PIMAX
  DATA PIMIN /445.0,445.0,444.0,412.0,390.0,325.0,300.0,185.0,128.0,85.0,70.0,66.0,430.0,430.0,430.0,404.0,385.0,333.0,305.0,195.0,&
  137.0,92.0,75.0,74.0,414.0,414.0,414.0,394.0,380.0,340.0,310.0,208.0,145.0,98.0,82.0,80.0,399.0,399.0,399.0,382.0,370.0,335.0,   &
  305.0,208.0,150.0,100.0,85.0,85.0,383.0,383.0,383.0,369.0,365.0,330.0,300.0,208.0,153.0,102.0,89.0,89.0,366.0,366.0,366.0,354.0, &
  350.0,312.0,290.0,208.0,156.0,105.0,94.0,91.0,349.0,349.0,349.0,339.0,335.0,308.0,285.0,208.0,156.0,107.0,95.0,93.0,332.0,332.0, &
  332.0,325.0,320.0,300.0,280.0,208.0,154.0,110.0,100.0,94.0,298.0,298.0,298.0,292.0,290.0,285.0,255.0,195.0,150.0,108.0,98.0,93.0,&
  266.0,266.0,266.0,261.0,264.0,264.0,230.0,185.0,142.0,105.0,95.0,91.0,181.0,181.0,181.0,181.0,181.0,181.0,173.0,135.0,111.0,80.0,&
  78.0,75.0,136.0,136.0,136.0,136.0,135.0,134.0,126.0,103.0,87.0,77.0,69.0,62.0,95.0,95.0,95.0,95.0,95.0,95.0,95.0,75.0,67.0,60.0, &
  50.0,48.0,11.4,11.4,11.4,11.4,11.4,11.4,11.4,10.6,10.4,10.0,10.0,10.0/!Table 13.2, pg. 533; also Wallace & Sondhaus, Table i. (units: ion pairs/cm3-sec)
  DATA PIMAX /265.0,265.0,265.0,264.0,264.0,264.0,235.0,163.0,95.0,78.0,66.0,63.0,268.0,268.0,265.0,265.0,265.0,265.0,238.0,168.0, &
  104.0,85.0,71.0,70.0,267.0,267.0,265.0,265.0,265.0,265.0,240.0,179.0,112.0,91.0,78.0,76.0,265.0,265.0,264.0,262.0,262.0,262.0,   &
  240.0,182.0,118.0,93.0,81.0,81.0,258.0,258.0,257.0,256.0,253.0,252.0,239.0,178.0,118.0,95.0,84.0,84.0,252.0,251.0,250.0,249.0,   &
  247.0,245.0,238.0,175.0,119.0,98.0,89.0,88.0,243.0,243.0,243.0,242.0,241.0,241.0,237.0,174.0,120.0,100.0,91.0,89.0,235.0,235.0,  &
  233.0,231.0,231.0,231.0,230.0,174.0,122.0,100.0,95.0,90.0,216.0,216.0,216.0,213.0,213.0,212.0,209.0,170.0,118.0,101.0,94.0,90.0, &
  197.0,197.0,197.0,197.0,197.0,197.0,197.0,160.0,117.0,98.0,91.0,87.0,145.0,145.0,145.0,145.0,145.0,145.0,145.0,159.0,100.0,75.0, &
  74.0,72.0,109.0,109.0,109.0,109.0,108.0,108.0,101.0,88.0,78.0,72.0,66.0,60.0,79.0,79.0,79.0,79.0,79.0,79.0,79.0,65.0,60.0,56.0,  &
  48.0,47.0,11.4,11.4,11.4,11.4,11.4,11.4,11.4,10.6,10.4,10.0,10.0,10.0/!Table 13.3, pg. 534 (units: ion pairs/cm3-sec)
  REAL(dp) :: TEXP
  EXTERNAL TEXP
  n = 1
  interpolate: DO i = 1, 14
    IF(zssf > XT(i)) n = i
    cutoff: DO j = 2, 11
      m = j
      IF(VCOsf <= RT(j)) EXIT
    END DO cutoff
    ratio1 = (PIMAX(m,i) - PIMAX(m - 1,i)) / (RT(m) - RT(m - 1))
    PIRX = PIMAX(m - 1,i) + ratio1 * (VCOsf - RT(m - 1))
    ratio2 = (PIMIN(m,i) - PIMIN(m - 1,i)) / (RT(m) - RT(m - 1))
    PIRM = PIMIN(m - 1,i) + ratio2 * (VCOsf - RT(m - 1))
    SLOPE = (PIRM - PIRX) / (98.0 - 80.0)!Wilson says "98.3" (pg. 533)
    PIR(i) = PIRX + SLOPE * (csf - 80.0)
  END DO interpolate
  IF(n >= 13) THEN
    arg = (XT(14) - XT(13)) / LOG(PIR(13) / PIR(14))!this way works best
    AIRIP = PIR(13) * TEXP(XT(13) / arg) * TEXP(-zssf / arg)
    RETURN
  ELSE
    n = n + 1
    SLOPE = (PIR(n) - PIR(n - 1)) / (XT(n) - XT(n - 1))
    AIRIP = PIR(n - 1) + SLOPE * (zssf - XT(n - 1))
    RETURN
  END IF
END

REAL(dp) FUNCTION TEXP(input)
  USE Shared, ONLY : dp
  IMPLICIT NONE
  REAL(dp), INTENT(IN) :: input
  REAL(dp) :: xx
  IF(input > 69.0) THEN
    xx = 69.0
  ELSE IF(input < -69.0) THEN
    xx = -69.0
  ELSE!do nothing
    xx = input
  END IF
  TEXP = EXP(xx)
END