glow_run.F90 Source File

submodule (ionization) glow_mod

use, intrinsic :: iso_fortran_env, only: sp=>real32
use phys_consts, only : wp

!> subroutine to allocate use-associated variables
use cglow,only: cglow_init, &
  jmax,nbins,lmax,nmaj,nei,nex,nw,nc,nst, &
  idate,ut,glat,glong,f107a,f107,f107p,ap,ef,ec, &
  iscale,jlocal,kchem,xuvfac, &
  sza,dip,efrac,ierr, &
  zz,zo,zn2,zo2,zns,znd,zno,ztn,ze,zti,zte, &
  ener,del,phitop,wave1,wave2,sflux,pespec,sespec,uflx,dflx,sion, &
  photoi,photod,phono,aglw,ecalc,zxden,zeta,zceta,eheat,vcb, &
  data_dir

implicit none

logical :: first_call = .true.

contains

module procedure glow_run

! This software is part of the GLOW model.  Use is governed by the Open Source
! Academic Research License Agreement contained in the file glowlicense.txt.
! For more information see the file glow.txt.

! Stan Solomon and Ben Foster, 1/15
! Stan Solomon, 12/15, 1/16
! Stan Solomon, 3/16, MPI parallel version
! Guy Grubbs II, 4/17, modified for GEMINI integration

! For definitions of use-associated variables, see subroutine GLOW and module CGLOW.
! For definitions of TGCM input variables see module READTGCM
! For definitions of output arrays see module OUTPUT

! Other definitions:
! f107p   Solar 10.7 cm flux for previous day
! ap      Ap index of geomagnetic activity
! z       altitude array, km

! Array dimensions:
! jmax    number of altitude levels
! nbins   number of energetic electron energy bins
! lmax    number of wavelength intervals for solar flux
! nmaj    number of major` species
! nst     number of states produced by photoionization/dissociation
! nei     number of states produced by electron impact
! nex     number of ionized/excited species
! nw      number of airglow emission wavelengths
! nc      number of component production terms for each emission


real(sp), dimension(:), allocatable :: phitoptmp
integer :: j
character(len=1024) :: iri90_dir

data_dir = DATADIR
iri90_dir = trim(data_dir) // '/iri90/'


!! Execute:  Allocate arrays in other modules (formerly in common blocks):

if(first_call) then
  first_call = .false.
  jmax=size(alt,1)
  nbins = 190     !< MZ - note nbins no longer set to a default value in glow; user MUST assign this!!!
  call cglow_init
end if

!! Set electron energy grid:

call egrid (ener, del, nbins)

!! Set Maxwellian distribution into phitop array

!! Hard coded solution, future = pass ec and ef array to maxt assuming > 2 populations
allocate(phitoptmp(nbins))
phitoptmp=0

phitop=phitoptmp
do j = 1, size(PhiWmWm2,1)    !this index loops over population number
  call maxt(real(PhiWmWm2(j), sp),real(W0(j), sp),ener,del,nbins,0,0,0,phitoptmp)
  phitop = phitop + phitoptmp
end do


!! Set variables given from GEMINI
glat = real(xlat, sp)
glong = real(xlon, sp)
idate = date_doy
ut = real(UTsec, sp)
f107 = real(xf107, sp)
f107p = real(xf107, sp)
f107a = real(xf107a, sp)
ap = 5.
kchem = 4.
jlocal = 0

!! Convert densities and altitudes into

zz(:)  = real(alt(:)*1e2_wp, sp)
zo(:)  = real(nn(:,1)/1e6_wp, sp)
zo2(:) = real(nn(:,3)/1e6_wp, sp)
zn2(:) = real(nn(:,2)/1e6_wp, sp)
zno(:) = real(nn(:,6)/1e6_wp, sp)
zns(:) = real(nn(:,5)/1e6_wp, sp)
znd(:) = real(0, sp)
ztn(:) = real(Tn(:), sp)

! ZXDEN   array of excited and and/or ionized state densities at each altitude:
!           O+(2P), O+(2D), O+(4S), N+, N2+, O2+, NO+, N2(A), N(2P),
!           N(2D), O(1S), O(1D); cm-3
!ions (ns): 1=O+, 2=NO+, 3=N2+, 4=O2+, 5=N+, 6=H+
!neutrals (nn): O,N2,O2,H,N,NO
zxden(1, :) = 0
zxden(2, :) = 0
zxden(3, :) = real(ns(:,1)/1e6_wp, sp)
zxden(4, :) = real(ns(:,5)/1e6_wp, sp)
zxden(5, :) = real(ns(:,3)/1e6_wp, sp)
zxden(6, :) = real(ns(:, sp)/1e6_wp, sp)
zxden(7, :) = real(ns(:,2)/1e6_wp, sp)
zxden(8, :) = 0
zxden(9, :) = 0
zxden(10,:) = 0
zxden(11,:) = 0
zxden(12,:) = 0
zti(:) = real((Ts(:,1)*ns(:,1)+Ts(:,2)*ns(:,2)+Ts(:,3)*ns(:,3)+Ts(:,4)*ns(:,4)+Ts(:,5)*ns(:,5)) &
        /(ns(:,1)+ns(:,2)+ns(:,3)+ns(:,4)+ns(:,5)), sp)
zte(:) = real(Ts(:,7), sp)
ze(:)  = real(ns(:,7)/1e6_wp, sp)

!! Call GLOW to calculate ionized and excited species, and airglow emission rates:

call glow

ionrate(:,1) = (real(SION(1,:),wp)+real(SION(2,:),wp)*0.3_wp)*1e6_wp !O+
ionrate(:,4) = (real(SION(2,:),wp)*0.7_wp)*1e6_wp !O2+
ionrate(:,3) = (real(SION(3,:),wp)*0.84_wp)*1e6_wp !N2+
ionrate(:,5) = (real(SION(3,:),wp)*0.16_wp)*1e6_wp !N+
ionrate(:,2) = real(0., wp)*1e6_wp !NO+
ionrate(:,6) = real(0., wp)*1e6_wp !H+

eheating = real(eheat,wp)*1e6_wp

iver = real(vcb,wp)


!  do j = 1, jmax
!    IDENS(j) = sum((dflx(:,j)-uflx(:,j))*del(:))/6.241509d14
!  enddo

!  if (first_out) then
!    first_out = .false.
!    write(6,"(1x,i7,11f10.3)") idate,ut,glat,glong,f107a,f107,f107p,ap,ef,ec,ef1,ec1
!    write(6,"('   Z     Tn   Ti   Te      O        N2        NO      Ne(in)    Ne(out)   Ionrate     &
!          O+        O2+       NO+       N2+     EHeat    Jz')")
!    write(6,"(1x,0p,f5.1,3f6.0,1p,12e10.2)") (alt(j)/1.0d3,ztn(j),zti(j),zte(j),zo(j),zn2(j),zno(j),ze(j), &
!      ecalc(j),tir(j),PO(j),PO2(j),PNO(j),PN2(j),EHEATING(j),IDENS(j),j=1,jmax)
!  end if

deallocate(phitoptmp)

end procedure glow_run

end submodule glow_mod