Release of GLOW v.0.981, Stan Solomon, 6/29/17

NCAR · Jun 29, 2017 · 6cb880d · 6cb880d
1 parent 58f3b40
commit 6cb880d
Show file tree

Hide file tree

Showing 16 changed files with 200 additions and 53 deletions.
diff --git a/Glow.txt b/Glow.txt
@@ -1,6 +1,6 @@
                   Preliminary Documentation for the GLOW Model
 
-                               Version 0.98
+                               Version 0.981
 
                ***  see Releasenotes.txt for release notes  ***
 
@@ -368,6 +368,17 @@ empirical model inputs, MSIS and IRI provide full altitude profiles, but NOEM
 only extends from 100-150 km, so a constant value is used for nitric oxide
 below 100 km, and a scale height assumption above 150 km.
 
+   For TIE- or TIME-GCM runs that will be used for input to GLOW, secondary
+history files should be employed.  It is possible to use primary histories
+but not recommended, since no auroral parameters will be included.  The
+following secondary history output fields should be specified: ZG, TN, TI, TE,
+NE, O1, O2, N2, NO, N4S, N2D, HE, UN, VN, EFLUX, NFLUX, ALFA, CUSP, DRIZZLE.
+Note that HE, UN, and VN are not used, but are passed through to
+output.  The auroral parameters EFLUX, NFLUX, ALFA, CUSP, and DRIZZLE are
+optional but should be included.  (CUSP and DRIZZLE are not yet used,
+but may be soon.)  If ZG or N2 are ommitted, they will be calculated by
+readtgcm, provided that HE is included.
+
    The default energy grid extends to 50 keV and is appropriate for most
 photoelectron and auroral calculations; the electron energy grid can extend
 (in principle) to 1 GeV.  The altitude grid may be altered to suit the user's
@@ -433,6 +444,18 @@ from GCM fields the NO calculated by the GCM is used.
    Electron impact cross sections generated by EXSECT can be obtained from
 variables use associated through the CGLOW module.
 
+   A note on units:  GLOW uses cgs in its internal calculations.  Wavelengths
+are described in nanometers in this publication, but Angstroms are used for 
+solar and airglow emission wavelengths inside the code.  Energies are in 
+electron volts (eV), except for the input auroral total energy flux, which is 
+expressed in erg cm-2 (equivalent to mW m-2).  Other units employed for 
+convenience are megabarns (10-18 cm2) and Rayleighs (106 photons cm-2).  
+Conductivities are output in Siemens/meter (mhos).  Kilometers are used in 
+driver program altitude grids, but are converted to centimeters internally.  
+Dates are in yyyyddd or optionally yyddd format, universal time 
+is in seconds, and local solar time is in hours.  The solar radio 
+flux index F10.7 is in expressed in units of 10-22 W m-2. 
+
    References/Citations:
 Nagy, A. F., and P. M. Banks, Photoelectron fluxes in the ionosphere,
    J. Geophys. Res., 75, 6260, 1970.
@@ -446,3 +469,5 @@ Solomon, S. C., S. M. Bailey, and T. N. Woods, Effect of solar soft X-rays
    on the lower atmosphere, Geophys. Res. Lett., 28, 2149, 2001.
 Bailey, S. M., C. A. Barth, and S. C. Solomon, A model of nitric oxide in the
    lower thermosphere, J. Geophys. Res., 107, 1205, 2002.
+Solomon, S. C., Global modeling of thermospheric airglow in the
+   far-ultraviolet, J. Geophys. Res., submitted, 2017.
diff --git a/Makefile.glowbasic.gfort b/Makefile.glowbasic.gfort
@@ -0,0 +1,30 @@
+#
+FC = gfortran
+
+ FFLAGS = -O3 -ffree-line-length-none
+#FFLAGS = -g $(DBGFLAGS)
+
+DBGFLAGS = -debug full -traceback
+DBGFLAGS += -check bounds -check format -check output_conversion -check pointers -check uninit
+DBGFLAGS += -fpe-all=0 # this traps all floating point exceptions
+
+.SUFFIXES: .o .F .F90 .f .f90
+
+%.o: %.f90
+	$(FC) $(FFLAGS) -c  -o $@ $<
+%.o: %.f
+	$(FC) $(FFLAGS) -c  -o $@ $<
+#
+# Sources (in order of dependency):
+#
+SOURCES = cglow.f90 glowbasic.f90 glow.f90 bands.f90 conduct.f90 egrid.f90 ephoto.f90 etrans.f90 exsect.f fieldm.f gchem.f90 geomag.f90 maxt.f90 mzgrid.f90 qback.f90 rcolum.f90 rout.f90 snoem.f90 snoemint.f90 solzen.f90 ssflux.f90 iri90.f nrlmsise00.f
+
+OBJS := $(addsuffix .o, $(basename $(SOURCES)))
+EXEC = glow.exe
+
+$(EXEC): $(OBJS)
+	$(FC) -o $@ $(OBJS)
+
+clean:
+	rm -f *.o *.mod $(EXEC)
+
diff --git a/Quickstart.txt b/Quickstart.txt
@@ -4,6 +4,7 @@ Copy all the files and sub-directories into a directory, or,
   copy the file glowv0.98.tar to a directory and tar -xf glowv0.98.tar
 While in that directory:
   cp Makefile.glowbasic Makefile
+    (or, cp Makefile.glowbasic.gfort Makefile if using GNU gfortran)
   gmake
 Run the glowbasic driver with default dayglow inputs:
   glow.exe < in.basic.day > test.basic.day

diff --git a/Releasenotes.txt b/Releasenotes.txt
@@ -1,3 +1,16 @@
+GLOW version 0.981 release notes, Stan Solomon, 6/2017
+
+Very minor update to v 0.98, implemented to facilitate using GNU gfortran:
+   Changed all open statements to use action='read' instead of 'readonly'
+   Commented out two lines from nrlmsise00.f; changed 0 to 0.
+   Replaced sind and cosd functions in snoem.f90, and general cleanup.
+   Corrected rcolum.f90 to call chap with pi-chi when chi>pi/2.
+   Added a trap to etrans.f90 to assure no calculation when j=1.
+   Zeroed out all arrays in cglow.f90.
+   None of these changes should have any effect on results.
+
+-----------------------------------------------------
+
 GLOW version 0.98 release notes, Stan Solomon, 3/2017
 
 See Quickstart.txt for very basic instructions
@@ -46,7 +59,7 @@ Other significant changes:
    The quartic equation solver formerly used by GCHEM to calculate electron
 density has been abandoned because it produced minor errors and could not
 be verified.  Instead, a simple iterative method is used, which turns out
-to be more accurate and equally fast.  Thanks to Matt Hirsch for identifying
+to be more accurate and equally fast.  Thanks to Michael Hirsch for identifying
 this problem.
 
    The high-energy electron precipitation example (hexexample) is no longer

diff --git a/cglow.f90 b/cglow.f90
@@ -4,10 +4,11 @@ module cglow
 ! Academic Research License Agreement contained in the file glowlicense.txt.
 ! For more information see the file glow.txt.
 
-! Version 0.98, 1/2017
+! Version 0.981, 6/2017
 
 ! Stan Solomon and Ben Foster, 1/2015
 ! Stan Solomon, 1/2016, 3/2016, consolidated with cxglow
+! Stan Solomon, 6/2017, zeroed out arrays
 
 ! CGLOW Defines array dimensions and use-associated variables for the GLOW model.
 ! Replaces the header file glow.h and common blocks /CGLOW/, /CXSECT/, and /CXPARS/
@@ -101,6 +102,8 @@ module cglow
 
   subroutine cglow_init
 
+! Allocate variable arrays:
+
     allocate        &
       (zz   (jmax), &
        zo   (jmax), &
@@ -186,6 +189,71 @@ subroutine cglow_init
        gams(nei,nmaj), &
        gamb(nei,nmaj))
 
+! Zero all allocated variable arrays:
+
+       zz   (:)     =0.
+       zo   (:)     =0.
+       zn2  (:)     =0.
+       zo2  (:)     =0.
+       zno  (:)     =0.
+       zns  (:)     =0.
+       znd  (:)     =0.
+       zrho (:)     =0.
+       ze   (:)     =0.
+       ztn  (:)     =0.
+       zti  (:)     =0.
+       zte  (:)     =0.
+       eheat(:)     =0.
+       tez  (:)     =0.
+       tei  (:)     =0.
+       tpi  (:)     =0.
+       tir  (:)     =0.
+       ecalc(:)     =0.
+       zxden(:,:)   =0.
+       zeta(:,:)    =0.
+       zceta(:,:,:) =0.
+       zlbh(:,:)    =0.
+       phitop(:)    =0.
+       ener  (:)    =0.
+       del   (:)    =0.
+       wave1(:)     =0.
+       wave2(:)     =0.
+       sflux(:)     =0.
+       pespec(:,:)  =0.
+       sespec(:,:)  =0.
+       uflx  (:,:)  =0.
+       dflx  (:,:)  =0.
+       zmaj(:,:)    =0.
+       zcol(:,:)    =0.
+       pia (:,:)    =0.
+       sion(:,:)    =0.
+       aglw  (:,:,:)=0.
+       photoi(:,:,:)=0.
+       photod(:,:,:)=0.
+       phono(:,:)   =0.
+       sigs(:,:)    =0.
+       pe  (:,:)    =0.
+       pin (:,:)    =0.
+       sigex(:,:,:) =0.
+       sigix(:,:,:) =0.
+       siga(:,:,:)  =0.
+       sec (:,:,:)  =0.
+       iimaxx(:)    =0.
+       ww  (:,:)    =0.
+       ao  (:,:)    =0.
+       omeg(:,:)    =0.
+       anu (:,:)    =0.
+       bb  (:,:)    =0.
+       auto(:,:)    =0.
+       thi (:,:)    =0.
+       ak  (:,:)    =0.
+       aj  (:,:)    =0.
+       ts  (:,:)    =0.
+       ta  (:,:)    =0.
+       tb  (:,:)    =0.
+       gams(:,:)    =0.
+       gamb(:,:)    =0.
+
   end subroutine cglow_init
 
 !-----------------------------------------------------------------------

diff --git a/conduct.f90 b/conduct.f90
@@ -122,7 +122,7 @@ SUBROUTINE CONDUCT(lat, long, alt, nO, nO2, nN2, nOp, nO2p, nNOp, &
 
 ! Calculate Pederson and Hall conductivity using the TIE-GCM formulation:
 ! NOTE: multiplied by 1.e6 to convert densities to m^-3, so that
-! conductivities are in Sieverts/meter, a.k.a. mhos.
+! conductivities are in Siemens/meter, a.k.a. mhos.
 
 
       PedCond = (1.e6*qe/Bmagn) *                                 &

diff --git a/ephoto.f90 b/ephoto.f90
@@ -88,7 +88,7 @@ subroutine ephoto
             rion(lmax,nmaj,jmax)
 
   real,parameter :: signo = 2.0e-18
-  integer,save :: ifirst=1
+  integer :: ifirst=1
   integer :: l,n,k,i,j,m,m1,m2
   real :: aa,bb,fac,e1,e2,y,r1,r2
   character(len=1024) :: filepath
@@ -115,7 +115,7 @@ subroutine ephoto
     ifirst = 0
 
     filepath = trim(data_dir)//'ephoto_xn2.dat'
-    open(unit=1,file=filepath,status='old',readonly)
+    open(unit=1,file=filepath,status='old',action='read')
     read(1,*)
     read(1,*)
     read(1,*)
@@ -126,7 +126,7 @@ subroutine ephoto
     close(1)
 
     filepath = trim(data_dir)//'ephoto_xo2.dat'
-    open(unit=1,file=filepath,status='old',readonly)
+    open(unit=1,file=filepath,status='old',action='read')
     read(1,*)
     read(1,*)
     read(1,*)
@@ -137,7 +137,7 @@ subroutine ephoto
     close(1)
 
     filepath = trim(data_dir)//'ephoto_xo.dat'
-    open(unit=1,file=filepath,status='old',readonly)
+    open(unit=1,file=filepath,status='old',action='read')
     read(1,*)
     read(1,*)
     read(1,*)

diff --git a/etrans.f90 b/etrans.f90
@@ -262,25 +262,27 @@ subroutine etrans
 !
 ! Cascade production:
 !
-      do k = 1, j-1
-        ll = j - k
-        produa(:)=0.
-        prodda(:)=0.
-        do n = 1, nmaj
+      if (j > 1) then
+        do k = 1, j-1
+          ll = j - k
+          produa(:)=0.
+          prodda(:)=0.
+          do n = 1, nmaj
+            do i=1,jmax
+              produa(i) = produa(i) &
+                         + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phidwn(i) &
+                         + (1. - pin(n,j))*siga(n,k,j)*phiup(i))
+              prodda(i) = prodda(i) &
+                         + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phiup(i) &
+                         + (1. - pin(n,j))*siga(n,k,j)*phidwn(i))
+            enddo
+          enddo
           do i=1,jmax
-            produa(i) = produa(i) &
-                       + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phidwn(i) &
-                       + (1. - pin(n,j))*siga(n,k,j)*phiup(i))
-            prodda(i) = prodda(i) &
-                       + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phiup(i) &
-                       + (1. - pin(n,j))*siga(n,k,j)*phidwn(i))
+            produp(i,ll) = produp(i,ll) + produa(i) * rmusin
+            prodwn(i,ll) = prodwn(i,ll) + prodda(i) * rmusin
           enddo
         enddo
-        do i=1,jmax
-          produp(i,ll) = produp(i,ll) + produa(i) * rmusin
-          prodwn(i,ll) = prodwn(i,ll) + prodda(i) * rmusin
-        enddo
-      enddo
+      endif
       kk = j - 1
       if (kk > 0) then
         do i = 1, jmax

diff --git a/glow.f90 b/glow.f90
@@ -5,6 +5,7 @@
 ! For more information see the file glow.txt.
 
 ! Version 0.98, 1/2017
+! Version 0.981, 6/2017
 
 ! Stan Solomon, 1988, 1989, 1990, 1991, 1992, 1994, 2000, 2002, 2005, 2015, 2016
 !

diff --git a/glowbasic.f90 b/glowbasic.f90
@@ -4,7 +4,7 @@ program glowbasic
 ! Academic Research License Agreement contained in the file glowlicense.txt.
 ! For more information see the file glow.txt.
 
-! Version 0.98, 1/2017
+! Version 0.981, 6/2017
 
 ! Adapted from glowdriver by Stan Solomon, 2/2016
 

diff --git a/glowdriver.f90 b/glowdriver.f90
@@ -4,7 +4,7 @@ program glowdriver
 ! Academic Research License Agreement contained in the file glowlicense.txt.
 ! For more information see the file glow.txt.
 
-! Version 0.98, 1/2017
+! Version 0.981, 6/2017
 
 ! Stan Solomon and Ben Foster, 1/2015
 ! Stan Solomon, 12/2015, 1/2016

diff --git a/nrlmsise00.f b/nrlmsise00.f
@@ -5,7 +5,7 @@ SUBROUTINE GTD7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)
 C     -----------
 C        Neutral Atmosphere Empirical Model from the surface to lower
 C        exosphere
-C
+C                   -- Quick Fixes for gfortran compatibility, SCS, 2017
 C        NEW FEATURES:
 C          *Extensive satellite drag database used in model generation
 C          *Revised O2 (and O) in lower thermosphere
@@ -276,7 +276,7 @@ SUBROUTINE GTD7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)
    10 CONTINUE
       GOTO 90
    50 CONTINUE
-      DD=DENSM(ALT,1.,0,TZ,MN3,ZN3,TN3,TGN3,MN2,ZN2,TN2,TGN2)                
+      DD=DENSM(ALT,1.,0.,TZ,MN3,ZN3,TN3,TGN3,MN2,ZN2,TN2,TGN2) ! 0. SCS 2017
       T(2)=TZ
    90 CONTINUE
       ALAST=ALT
@@ -1673,8 +1673,8 @@ BLOCK DATA GTD7BK
       COMMON/DATIM7/ISDATE(3),ISTIME(2),NAME(2)
       COMMON/METSEL/IMR
       DATA IMR/0/
-      DATA ISDATE/'01-F','EB-0','2   '/,ISTIME/'15:4','9:27'/
-      DATA NAME/'MSIS','E-00'/
+C     DATA ISDATE/'01-F','EB-0','2   '/,ISTIME/'15:4','9:27'/ ! cout, SCS, 2017
+C     DATA NAME/'MSIS','E-00'/                                ! cout, SCS, 2017
 C         TEMPERATURE
       DATA PT1/
      *  9.86573E-01, 1.62228E-02, 1.55270E-02,-1.04323E-01,-3.75801E-03,

diff --git a/rcolum.f90 b/rcolum.f90
@@ -9,6 +9,7 @@
 ! boundary.  If grazing height is below lower boundary of atmosphere
 ! supplied, column density is set to 1.0e30.
 ! Stan Solomon, 2016: refactored for f90.
+! Stan Solomon, 2017: corrected chi to pi-chi at line 79.  Note this has no effect.
 
 ! Calculates the column density ZCOL for each species ZMAJ above height
 ! ZZ at zenith angle CHI.  Calls subroutine VCD to calculate the
@@ -31,7 +32,7 @@ subroutine rcolum (chi, zz, zmaj, tn, zcol, zvcd, jmax, nmaj)
       real,intent(out) :: zcol(nmaj,jmax), zvcd(nmaj,jmax)
 
       integer,parameter :: nm=3
-      real,parameter :: pi=3.1415926535
+      real,parameter :: pi=3.1415926536
       real,parameter :: re=6.37e8
 
       integer :: i, j, k
@@ -75,7 +76,7 @@ subroutine rcolum (chi, zz, zmaj, tn, zcol, zvcd, jmax, nmaj)
             enddo
             do i=1,nmaj
               zcol(i,j) = 2. * zcg(i) * chap(pi/2.,ghz,tng,i) &
-                        - zvcd(i,j) * chap(chi,zz(j),tn(j),i)
+                        - zvcd(i,j) * chap(pi-chi,zz(j),tn(j),i)
             enddo
           endif
         enddo
@@ -95,7 +96,7 @@ real function chap (chi, z, t, i)
       integer,intent(in) :: i
 
       integer,parameter ::  nmaj=3
-      real,parameter :: pi=3.1415926535
+      real,parameter :: pi=3.1415926536
       real,parameter :: re=6.37e8
       real,parameter :: g=978.1