module readtgcm

! 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.

! Ben Foster and Stan Solomon, 1/15

! Stan Solomon, 12/15, 1/16, minor updates

! Stan Solomon, 3/16, addition of helium

! Stan Solomon, 11/16, bug fix to prevent reading variables not found on file

! Reads TIE-GCM or TIME-GCM history file, obtaining needed inputs for GLOW.

! Usage (see program glowdriver.f90):

! First call read_tgcm_coords to get coordinates from history file

! Next call find_mtimes to find indices of model times on history file

! Then call read_tgcm (up to 100 times) to read individual history data

! read_tgcm_coords returns:

! nlon Number of longitudes

! nlat Number of latitudes

! nlev Number of midpoint levels

! nilev Number of interface levels

! glon Array of longitudes

! glat Array of latitudes

! zlev Array of midpoints

! zilev Array of interfaces

! read_tgcm returns:

! 3D global fields:

! zg Geometric height of pressure surface, km

! tn Neutral temperature, K

! ti Ion temperature, K

! te Electron temperature, K

! ne Electron number density, cm-3

! o2 O2 number density, cm-3

! o1 O number density, cm-3

! n2 N2 number density, cm-3

! no NO number density, cm-3

! n4s N(4S) number density, cm-3

! n2d N(2D) number density, cm-3 (optional - used to initialize GLOW but not necessary)

! he He number density, cm-3 (optional - not used by GLOW)

! z Geopotential height, km (not used by GLOW)

! un Zonal wind, cm/s (not used by GLOW but passed through to output)

! vn Meridional winds, cm/s (not used by GLOW but passed through to output)

! 2D global fields:

! eflux Energy flux of auroral electrons, erg cm-2 s-1

! nflux Number flux of auroral electrons, cm-2 s-1

! alfa Characteristic energy of auroral electron flux, keV

! cusp Multiplier for cusp electron flux

! drizzle Multiplier for polar drizzle electron flux

! Vectors as a function of time:

! ec Cusp maximum energy flux, erg cm-2 s-1

! alfac Characteristic energy of cusp electron flux, keV

! ed Drizzle maximum energy flux, erg cm-2 s-1

! alfad Characteristic energy of drizzle electron flux, keV

! hpower Hemispheric power of auroral electron flux, GW

! f107d Solar F10.7 daily index

! f107a Solar F10.7 81-day centered average

! e1 Auroral oval parameter

! e2 Auroral oval parameter

! h1 Auroral oval parameter

! h2 Auroral oval parameter

! Scalar values for each call:

! iyear Model year, yyyy

! iday Model day, ddd

! ut Model universal time, hours

use netcdf

implicit none

integer :: iyear,iday

real :: ut,dlev

!

! Grid dimensions and coordinates:

!

integer :: nlon,nlat,nlev,nilev ! grid dimensions

integer :: ntime ! number of histories on file

real,allocatable,dimension(:),save :: & ! grid coordinates

glon,glat,zlev,zilev

integer,parameter :: mxtimes=100

!

! 2d and 3d global fields at a single time (history):

!

real,allocatable,dimension(:,:,:),save :: & ! 3d fields (nlon,nlat,nlev)

tn,un,vn,o2,o1,n2,he,no,ti,te,z,zg,ne,n2d,n4s

real,allocatable,dimension(:,:),save :: & ! 2d fields (nlon,nlat)

cusp,drizzle,alfa,nflux,eflux

!

! Scalars (1d in time):

!

real,allocatable,dimension(:),save :: & ! scalars (ntime)

f107d,f107a,hpower,e1,e2,h1,h2,alfac,ec,alfad,ed

integer,parameter :: nf=20 ! number of fields to read

character(len=8) :: fnames(nf) ! field names on file

!

! Print to stdout flag:

!

integer,save :: iprint=1

contains

!-----------------------------------------------------------------------

subroutine read_tgcm(ncfile,itime)

!

! Read needed fields from tgcm history file. Note sub read_tgcm_coords

! must be called before this is called.

!

! Args:

character(len=*),intent(in) :: ncfile ! file path

integer,intent(in) :: itime ! time index to read

!

! Local:

integer :: n,istat,ncid,k,i,ndims

integer :: id,idv_lon,idv_lat,idv_lev,idv_ilev,idv_time

character(len=1024) :: msg

character(len=NF90_MAX_NAME) :: varname

integer :: iyear1(1),iday1(1)

real :: ut1(1)

real,allocatable,save :: f3d(:,:,:) ! global 3d field (lon,lat,lev)

real,allocatable,save :: f2d(:,:) ! global 2d field (lon,lat)

real,allocatable,dimension(:,:,:),save :: o2mmr,o1mmr,n2mmr,hemmr

integer :: mtime(3)

logical :: found(nf)

!

! These fields with these names must be on the file as 3d+time variables:

! Note that CUSP,DRIZZLE,ALFA,NFLUX,EFLUX are 2d (lon,lat).

!

fnames = (/'TN ','UN ','VN ','O2 ','O1 ', &

'NO ','TI ','TE ','Z ','ZG ', &

'NE ','N2 ','HE ','N2D ','N4S ', &

'CUSP ','DRIZZLE ','ALFA ','NFLUX ','EFLUX '/)

istat = nf90_open(ncfile,NF90_NOWRITE,ncid)

if (istat /= NF90_NOERR) then

write(msg,"('Error opening file ',a)") trim(ncfile)

call handle_ncerr(istat,trim(msg),1)

else

if (iprint > 0) &

write(6,"(/,'Opened file ',a)") trim(ncfile)

endif

!

! Allocate arrays to read 2d and 3d spatial variables:

!

if (.not.allocated(f3d)) allocate(f3d(nlon,nlat,nlev)) ! 3d global read array

if (.not.allocated(f2d)) allocate(f2d(nlon,nlat)) ! 2d global read array

!

! Read itime history (model time mtime(:,itime)):

!

if (itime <= 0.or.itime > ntime) then

write(6,"('>>> read_tgcm: bad itime=',i4,' (ntime=',i4,')')") &

itime,ntime

stop 'itime'

endif

istat = nf90_inq_varid(ncid,'year',id)

istat = nf90_get_var(ncid,id,iyear1,(/itime/),(/1/))

iyear = iyear1(1)

istat = nf90_inq_varid(ncid,'day',id)

istat = nf90_get_var(ncid,id,iday1,(/itime/),(/1/))

iday = iday1(1)

istat = nf90_inq_varid(ncid,'ut',id)

istat = nf90_get_var(ncid,id,ut1,(/itime/),(/1/))

ut = ut1(1)

istat = nf90_inq_varid(ncid,'mtime',id)

istat = nf90_get_var(ncid,id,mtime,(/1,itime/),(/3,1/))

write(6,"(/,'read_tgcm: itime=',i4,' ntime=',i4,' iyear=',i5,' iday=',i4,' ut=',f8.2,' mtime=',3i4)") &

itime,ntime,iyear,iday,ut,mtime

!

! Read vars at current history/time:

!

found(:) = .false.

do n=1,nf

istat = nf90_inq_varid(ncid,fnames(n),id)

if (istat == NF90_NOERR) then ! read variable only if found, scs 11/16

istat = nf90_inquire_variable(ncid,id,ndims=ndims)

!

! 3d var+time: assume (nlon,nlat,nlev):

!

if (ndims==4) then

f3d = 0.

istat = nf90_get_var(ncid,id,f3d,(/1,1,1,itime/),(/nlon,nlat,nlev,1/))

if (istat /= NF90_NOERR) then

write(6,"('>>> Error reading 3d var ')") fnames(n)

call handle_ncerr(istat,'Error from nf90_get_var',1)

endif

if (iprint > 0) &

write(6,"('Read 3d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") &

itime,n,trim(fnames(n)),id,minval(f3d(:,:,1:nlev-1)),maxval(f3d(:,:,1:nlev-1))

!

! 2d var+time: assume (nlon,nlat):

!

elseif (ndims==3) then

f2d = 0.

istat = nf90_get_var(ncid,id,f2d,(/1,1,itime/),(/nlon,nlat,1/))

if (istat /= NF90_NOERR) then

write(6,"('>>> Error reading 2d var ')") fnames(n)

call handle_ncerr(istat,'Error from nf90_get_var',1)

endif

if (iprint > 0) &

write(6,"('Read 2d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") &

itime,n,trim(fnames(n)),id,minval(f2d(:,:)),maxval(f2d(:,:))

endif ! 2d or 3d var

!

! Transfer field to module data:

! Top level may be missing in tn,un,vn,ti,te, so set these to values at nlev-1 for now:

!

select case(trim(fnames(n)))

case('TN')

if (.not.allocated(tn)) allocate(tn(nlon,nlat,nlev))

tn = f3d

tn(:,:,nlev) = tn(:,:,nlev-1)

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(tn),maxval(tn)

found(n) = .true.

case('UN')

if (.not.allocated(un)) allocate(un(nlon,nlat,nlev))

un = f3d

un(:,:,nlev) = un(:,:,nlev-1)

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(un),maxval(un)

found(n) = .true.

case('VN')

if (.not.allocated(vn)) allocate(vn(nlon,nlat,nlev))

vn = f3d

vn(:,:,nlev) = vn(:,:,nlev-1)

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(vn),maxval(vn)

found(n) = .true.

case('O2')

if (.not.allocated(o2)) allocate(o2(nlon,nlat,nlev))

o2 = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(o2),maxval(o2)

found(n) = .true.

case('O1')

if (.not.allocated(o1)) allocate(o1(nlon,nlat,nlev))

o1 = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(o1),maxval(o1)

found(n) = .true.

case('N2')

if (.not.allocated(n2)) allocate(n2(nlon,nlat,nlev))

n2 = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(n2),maxval(n2)

found(n) = .true.

case('HE')

if (.not.allocated(he)) allocate(he(nlon,nlat,nlev))

he = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(he),maxval(he)

found(n) = .true.

case('NO')

if (.not.allocated(no)) allocate(no(nlon,nlat,nlev))

no = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(no),maxval(no)

found(n) = .true.

case('TI')

if (.not.allocated(ti)) allocate(ti(nlon,nlat,nlev))

ti = f3d

ti(:,:,nlev) = ti(:,:,nlev-1)

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(ti),maxval(ti)

found(n) = .true.

case('TE')

if (.not.allocated(te)) allocate(te(nlon,nlat,nlev))

te = f3d

te(:,:,nlev) = te(:,:,nlev-1)

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(te),maxval(te)

found(n) = .true.

case('NE')

if (.not.allocated(ne)) allocate(ne(nlon,nlat,nlev))

ne = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(ne),maxval(ne)

found(n) = .true.

case('Z')

if (.not.allocated(z)) allocate(z(nlon,nlat,nlev))

z = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(z),maxval(z)

found(n) = .true.

case('ZG')

if (.not.allocated(zg)) allocate(zg(nlon,nlat,nlev))

zg = f3d*1.e-5 ! cm->km

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(zg),maxval(zg)

found(n) = .true.

case('N2D')

if (.not.allocated(n2d)) allocate(n2d(nlon,nlat,nlev))

n2d = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(n2d),maxval(n2d)

found(n) = .true.

case('N4S')

if (.not.allocated(n4s)) allocate(n4s(nlon,nlat,nlev))

n4s = f3d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(n4s),maxval(n4s)

found(n) = .true.

!

! 2d fields (nlon,nlat):

!

case('CUSP')

if (.not.allocated(cusp)) allocate(cusp(nlon,nlat))

cusp = f2d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(cusp),maxval(cusp)

found(n) = .true.

case('DRIZZLE')

if (.not.allocated(drizzle)) allocate(drizzle(nlon,nlat))

drizzle = f2d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(drizzle),maxval(drizzle)

found(n) = .true.

case('ALFA')

if (.not.allocated(alfa)) allocate(alfa(nlon,nlat))

alfa = f2d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(alfa),maxval(alfa)

found(n) = .true.

case('NFLUX')

if (.not.allocated(nflux)) allocate(nflux(nlon,nlat))

nflux = f2d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(nflux),maxval(nflux)

found(n) = .true.

case('EFLUX')

if (.not.allocated(eflux)) allocate(eflux(nlon,nlat))

eflux = f2d

if (iprint > 0) &

write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") &

fnames(n),minval(eflux),maxval(eflux)

found(n) = .true.

!

! Unknown field (this should not happen):

!

case default

write(6,"('>>> read_tgcm: unknown field: ',a)") trim(fnames(n))

stop 'read tgcm 3d fields'

end select ! fnames(n)

else

write(6,"('>>> WARNING: did not find field ',a,' on the tgcm file')") trim(fnames(n))

endif ! bottom of found variable conditional, scs 11/16

enddo ! n=1,nf

!

! Define n2 mmr from o2,o,he if necessary:

!

if (.not.found(findx('N2'))) then

if (.not.allocated(N2)) allocate(n2(nlon,nlat,nlev))

if (found(findx('O1')).and.found(findx('O2')).and.found(findx('HE'))) then

n2 = 1.-o2-o1-he

if (iprint > 0) &

write(6,"('read_tgcm: Field N2 (1-O2-O-HE) min,max=',2es12.4)") &

minval(n2),maxval(n2)

else

write(6,"('>>> FATAL: N2 not found on the file, and at least one of')")

write(6,"(' o2, o, he also not found on the file.')")

stop 'read_tgcm'

endif

found(findx('N2'))=.true.

endif ! N2 not found

!

! If ZG not on the file, calculate it using sub calczg (from tgcm addiag.F):

!

if (.not.found(findx('ZG'))) then

if (.not.found(findx('TN')).or..not.found(findx('O2')).or.&

.not.found(findx('O1')).or..not.found(findx('N2')).or.&

.not.found(findx('Z' )).or..not.found(findx('HE'))) then

write(6,"('>>> FATAL: ZG not found on the file, and at least one of')")

write(6,"(' tn,o2,o,n2,he,z also not found on the file.')")

stop 'read_tgcm'

endif

if (.not.allocated(zg)) allocate(zg(nlon,nlat,nlev))

call calczg(tn,o2,o1,n2,he,z,zg,glat,nlon,nlat,nlev,dlev)

zg = zg*1.e-5 ! cm to km

if (iprint > 0) &

write(6,"('Field ZG (from calczg) min,max (km)=',2es12.4)") &

minval(zg),maxval(zg)

found(findx('ZG'))=.true.

endif ! zg not on the file

!

! If N2D was not found, set it to zero:

!

if (.not.found(findx('N2D'))) then

if (.not.allocated(n2d)) allocate(n2d(nlon,nlat,nlev))

n2d = 0.

write(6,"('read_tgcm: N2D not found on tgcm history, so set it to zero.')")

found(findx('N2D'))=.true.

endif

!

! If HE was not found, set it to zero and issue warning:

!

if (.not.found(findx('HE'))) then

if (.not.allocated(he)) allocate(he(nlon,nlat,nlev))

he = 0.

write(6,"('read_tgcm: HE not found on tgcm history, so set it to zero.')")

write(6,"('>>> WARNING: This may cause erroneous O density at high altitude')")

found(findx('HE'))=.true.

endif

!

! Auroral fields not found can be set to zero:

!

if (.not.found(findx('CUSP'))) then

if (.not.allocated(cusp)) allocate(cusp(nlon,nlat))

cusp = 0.

write(6,"('read_tgcm: CUSP not found on tgcm history, so set it to zero.')")

found(findx('CUSP'))=.true.

endif

if (.not.found(findx('DRIZZLE'))) then

if (.not.allocated(drizzle)) allocate(drizzle(nlon,nlat))

drizzle = 0.

write(6,"('read_tgcm: DRIZZLE not found on tgcm history, so set it to zero.')")

found(findx('DRIZZLE'))=.true.

endif

if (.not.found(findx('ALFA'))) then

if (.not.allocated(alfa)) allocate(alfa(nlon,nlat))

alfa = 0.

write(6,"('read_tgcm: ALFA not found on tgcm history, so set it to zero.')")

found(findx('ALFA'))=.true.

endif

if (.not.found(findx('NFLUX'))) then

if (.not.allocated(nflux)) allocate(nflux(nlon,nlat))

nflux = 0.

write(6,"('read_tgcm: NFLUX not found on tgcm history, so set it to zero.')")

found(findx('NFLUX'))=.true.

endif

if (.not.found(findx('EFLUX'))) then

if (.not.allocated(eflux)) allocate(eflux(nlon,nlat))

eflux = 0.

write(6,"('read_tgcm: EFLUX not found on tgcm history, so set it to zero.')")

found(findx('EFLUX'))=.true.

endif

!

! If any other fields were not found on the history, stop w/ fatal error:

!

do n=1,nf

if (.not.found(n)) then

write(6,"(/,'>>> FATAL: field ',a,' was not found on tgcm file ',a)") &

fnames(n),trim(ncfile)

stop 'read_tgcm'

endif

enddo

!

! Convert species from mmr to cm3 (still in the time/history loop):

!

if (.not.allocated(o2mmr)) allocate(o2mmr(nlon,nlat,nlev))

if (.not.allocated(o1mmr)) allocate(o1mmr(nlon,nlat,nlev))

if (.not.allocated(n2mmr)) allocate(n2mmr(nlon,nlat,nlev))

if (.not.allocated(hemmr)) allocate(hemmr(nlon,nlat,nlev))

o2mmr = o2

o1mmr = o1

n2mmr = n2

hemmr = he

call denconv(o2, 32.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(o1, 16.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(n2, 28.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(he, 4.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(no, 30.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(n2d,14.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

call denconv(n4s,14.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev)

if (iprint > 0) then

write(6,"('After denconv: O2 (cm3) min,max=',2es12.4)") minval(o2),maxval(o2)

write(6,"('After denconv: O1 (cm3) min,max=',2es12.4)") minval(o1),maxval(o1)

write(6,"('After denconv: N2 (cm3) min,max=',2es12.4)") minval(n2),maxval(n2)

write(6,"('After denconv: HE (cm3) min,max=',2es12.4)") minval(he),maxval(he)

write(6,"('After denconv: NO (cm3) min,max=',2es12.4)") minval(no),maxval(no)

write(6,"('After denconv: N2D(cm3) min,max=',2es12.4)") minval(n2d),maxval(n2d)

write(6,"('After denconv: N4S(cm3) min,max=',2es12.4)") minval(n4s),maxval(n4s)

endif

!

! Close the file:

!

istat = nf90_close(ncid)

end subroutine read_tgcm

!-----------------------------------------------------------------------

subroutine calczg(tn,o2,o1,n2,he,z,zg,glat,nlon,nlat,nlev,dlev)

! use params_module,only: dz,glat

! use init_module,only: istep

! use cons_module,only: boltz,avo

! use addfld_module,only: addfld

!

! Given geopotential z (calculated with the model constant gravity),

! calculate geopotential zg with varying gravity. This is taken from

! tgcmproc_f90, routines calchts and glatf in proclat.F.

!

! Args:

integer,intent(in) :: nlon,nlat,nlev

real,dimension(nlon,nlat,nlev),intent(in):: &

tn, & ! neutral temperature (deg K)

o2, & ! molecular oxygen (mmr)

o1, & ! atomic oxygen (mmr)

n2, & ! molecular nitrogen (mmr)

he, & ! atomic helium (mmr)

z ! geopotential calculated with constant gravity (from addiag)

real,intent(in) :: glat(nlat),dlev

real,dimension(nlon,nlat,nlev),intent(out) :: &

zg ! output geopotential calculated with varying gravity

!

! Local:

integer :: i,j,k

real :: g0,r0,c2

real,dimension(nlev) :: xmas,g

real,parameter :: dgtr = 1.74533E-2

real,parameter :: avo = 6.023e23

real,parameter :: boltz = 1.38E-16

!

! Latitude scan:

do j=1,nlat

c2 = cos(2.*dgtr*glat(j))

g0 = 980.616*(1.-.0026373*c2)

r0 = 2.*g0/(3.085462e-6 + 2.27e-9*c2) ! effective earth radius

!

! Longitude scan:

do i=1,nlon

g(1)=g0*(r0/(r0+0.5*(z(i,j,1)+z(i,j,2))))**2

xmas(:) = 1./(o1(i,j,:)/16.+o2(i,j,:)/32.+ &

n2(i,j,:)/28.+he(i,j,:)/4.)/avo

!

! Levels:

zg(i,j,1) = z(i,j,1)

do k=2,nlev-1

zg(i,j,k) = zg(i,j,k-1) + boltz*dlev*tn(i,j,k-1) / &

(xmas(k-1)*g(k-1))

g(k)=g0*(r0/(r0+0.5*(zg(i,j,k)+z(i,j,k+1))))**2

enddo ! k=lev0+1,lev1-1

zg(i,j,nlev) = 1.5*zg(i,j,nlev-1)-0.5*zg(i,j,nlev-2)

enddo ! i=lon0,lon1

enddo ! j=1,nlat

end subroutine calczg

!-----------------------------------------------------------------------

subroutine denconv(f,wt,tn,o2,o1,n2,he,zlev,nlon,nlat,nlev)

!

! Convert density input field f from mmr to cm3.

! Inputs o2,o1,n2,he are in mmr.

! zlev is the vertical coordinate (zp).

!

implicit none

!

! Args:

integer,intent(in) :: nlon,nlat,nlev

real,intent(inout) :: f(nlon,nlat,nlev)

real,dimension(nlon,nlat,nlev),intent(in) :: tn,o2,o1,n2,he ! mmr

real,intent(in) :: wt

!

! Assume levels coordinate is increasing from bottom to top, regular dz)

real,intent(in) :: zlev(nlev) ! levels coord (tgcm zp)

!

! Local:

integer :: i,j,k

real :: zp,dlev,mbar,pkt

real,parameter :: p0=5.e-4 ! microbars

real,parameter :: boltz=1.3805e-16

dlev = (zlev(nlev)-zlev(1))/(nlev-1)

do k=1,nlev

zp = zlev(1)+(k-1)*dlev

do j=1,nlat

do i=1,nlon

mbar = 1./(o2(i,j,k)/32.+o1(i,j,k)/16.+n2(i,j,k)/28.+he(i,j,k)/4.)

pkt = p0*exp(-zp)/(boltz*tn(i,j,k))

f(i,j,k) = f(i,j,k) * pkt * mbar / wt

enddo

enddo

enddo

end subroutine denconv

!-----------------------------------------------------------------------

integer function findx(name)

!

! Find index in fnames(nf) that matches input name.

!

character(len=*),intent(in) :: name

integer :: n

findx = 0

do n=1,nf

if (trim(fnames(n))==name) then

findx = n

exit

endif

enddo

if (findx==0) then

write(6,"('>>> findx: could not find index to field ',a)") name

stop 'findx'

endif

end function findx

!-----------------------------------------------------------------------

integer function find_mtimes(ncfile,start_mtime,stop_mtime,mtimes,itimes)

!

! Open ncfile, find start and stop model times. Return all mtimes

! from start to stop, and return in mtimes(mxtimes). Function value

! is the number of mtimes.

!

! Args:

character(len=*),intent(in) :: ncfile

integer,intent(in) :: start_mtime(3),stop_mtime(3)

integer,intent(out) :: &

mtimes(3,mxtimes), & ! model times from start to stop

itimes(mxtimes) ! time index on file from start to stop

!

! Local:

integer :: i,ncid,istat,ntimes,id,idunlim,itime0,itime1

character(len=1024) :: msg

character(len=NF90_MAX_NAME) :: varname

integer,allocatable,save :: mtime_file(:,:) ! (3,ntime)

integer :: mtime(3)

!

! Init return values:

find_mtimes = 0

mtimes = 0

itimes = 0

!

! Open the file:

istat = nf90_open(ncfile,NF90_NOWRITE,ncid)

if (istat /= NF90_NOERR) then

write(msg,"('>>> find_mtimes: Error opening file ',a)") trim(ncfile)

call handle_ncerr(istat,trim(msg),1)

endif

!

! Get number of times (unlimited variable):

istat = nf90_inq_dimid(ncid,'time',idunlim)

istat = nf90_inquire_dimension(ncid,idunlim,varname,ntimes)

!

! Get model times on the file, and close:

if (.not.allocated(mtime_file)) allocate(mtime_file(3,ntimes))

istat = nf90_inq_varid(ncid,'mtime',id)

istat = nf90_get_var(ncid,id,mtime_file,(/1,1/),(/3,ntimes/))

istat = nf90_close(ncid)

!

! Search for start time:

itime0 = 0

if (start_mtime(1) == -999) then ! read first history on the file

itime0 = 1

else

do i=1,ntimes

if (mtime_file(1,i)==start_mtime(1).and. &

mtime_file(2,i)==start_mtime(2).and. &

mtime_file(3,i)==start_mtime(3)) then

itime0 = i

write(6,"('Found requested start_mtime ',3i4,' (history ',i4,' on the file).')") &

start_mtime,itime0

exit

endif

enddo

endif

if (itime0==0) then

write(6,"('>>> Could not find tgcm start_mtime ',3i4,' mtimes on the file are as follows:')") &

start_mtime

do i=1,ntimes

write(6,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i)

enddo

stop 'start_mtime'

endif

!

! Search for stop time:

itime1 = 0

if (stop_mtime(1) == -999) then ! stop time is last history

itime1 = ntimes

else

do i=1,ntimes

if (mtime_file(1,i)==stop_mtime(1).and. &

mtime_file(2,i)==stop_mtime(2).and. &

mtime_file(3,i)==stop_mtime(3)) then

itime1 = i

write(6,"('Found requested stop_mtime ',3i4,' (history ',i4,' on the file).')") &

stop_mtime,itime1

exit

endif

enddo

endif

if (itime1==0) then

write(6,"('>>> Could not find tgcm stop_mtime ',3i4,' mtimes on the file are as follows:')") &

stop_mtime

do i=1,ntimes

write(6,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i)

enddo

stop 'start_mtime'

endif

!

! This should not happen, but you never know...

if (itime1 < itime0) then

write(6,"('>>> find_mtimes: bad itime0 must be <= itime1: itime0,1=',2i4)") &

itime0,itime1

write(6,"('>>> mtime_file(itime0)=',3i4,' mtime_file(itime1)=',3i4)") &

mtime_file(:,itime0),mtime_file(:,itime1)

stop 'itime0,1'

endif

!

! Return values:

find_mtimes = itime1-itime0+1 ! number of model times

write(6,"('Number of model times to read = ',i4)") find_mtimes

do i=itime0,itime1

mtimes(:,i-itime0+1) = mtime_file(:,i) ! model times from start to stop

itimes(i-itime0+1) = i ! file index of each model time

enddo

end function find_mtimes

!-----------------------------------------------------------------------

subroutine read_tgcm_coords(ncfile)

!

! Read dimensions, coordinates, and 1d time variables from tgcm history file:

! (this must be called before tgcm_read is called)

!

! Args:

character(len=*),intent(in) :: ncfile ! file path

!

! Local:

integer :: istat,ncid,id,idunlim

integer :: idv_lon,idv_lat,idv_lev,idv_ilev,idv_time

character(len=1024) :: msg

character(len=NF90_MAX_NAME) :: varname

real,allocatable,save :: time(:)

istat = nf90_open(ncfile,NF90_NOWRITE,ncid)

if (istat /= NF90_NOERR) then

write(msg,"('Error opening file ',a)") trim(ncfile)

call handle_ncerr(istat,trim(msg),1)

else

write(6,"(/,'Opened file ',a)") trim(ncfile)

endif

!

! Get number of times (histories) (length of unlimited variable):

!

istat = nf90_inq_dimid(ncid,'time',idunlim)

istat = nf90_inquire_dimension(ncid,idunlim,varname,ntime)

! write(6,"('read_tgcm_coords: ntime=',i4)") ntime

istat = nf90_inq_dimid(ncid,'lon',id)

istat = nf90_inquire_dimension(ncid,id,varname,nlon)

istat = nf90_inq_dimid(ncid,'lat',id)

istat = nf90_inquire_dimension(ncid,id,varname,nlat)

istat = nf90_inq_dimid(ncid,'lev',id)

istat = nf90_inquire_dimension(ncid,id,varname,nlev)

istat = nf90_inq_dimid(ncid,'ilev',id)

istat = nf90_inquire_dimension(ncid,id,varname,nilev)

!

! Allocate coordinate variables:

!

allocate(glon(nlon))

allocate(glat(nlat))

allocate(zlev(nlev))

allocate(zilev(nilev))

allocate(time(ntime))

!

! Read coordinate vars:

!

istat = nf90_inq_varid(ncid,'lon',idv_lon)

istat = nf90_get_var(ncid,idv_lon,glon,(/1/),(/nlon/))

istat = nf90_inq_varid(ncid,'lat',idv_lat)

istat = nf90_get_var(ncid,idv_lat,glat,(/1/),(/nlat/))

istat = nf90_inq_varid(ncid,'lev',idv_lev)

istat = nf90_get_var(ncid,idv_lev,zlev,(/1/),(/nlev/))

istat = nf90_inq_varid(ncid,'ilev',idv_ilev)

istat = nf90_get_var(ncid,idv_ilev,zilev,(/1/),(/nilev/))

dlev = (zilev(nilev)-zilev(1))/(nilev-1)

istat = nf90_inq_varid(ncid,'time',idv_time)

istat = nf90_get_var(ncid,idv_time,time,(/1/),(/ntime/))

! write(6,"('read_tgcm_coords: nlon =',i4,' glon =',/,(8f10.3))") nlon,glon

! write(6,"('read_tgcm_coords: nlat =',i4,' glat =',/,(8f10.3))") nlat,glat

! write(6,"('read_tgcm_coords: nlev =',i4,' zlev =',/,(8f10.3))") nlev,zlev

! write(6,"('read_tgcm_coords: nilev=',i4,' zilev=',/,(8f10.3))") nilev,zilev

! write(6,"('read_tgcm_coords: nilev=',i4,' zilev(nilev)=',f9.3,' zilev(1)=',f9.3,' dlev=',f9.3)") &

! nilev,zilev(nilev),zilev(1),dlev

write(6,"('read_tgcm_coords: ntime=',i4,' nlon=',i4,' nlat=',i4,' nlev=',i4)") ntime,nlon,nlat,nlev

!

! Allocate and read 1d time-dependent variables (or, this could be done

! on a per-history basis below, reading into scalars instead of arrays

! (like ut,iday,iyear, etc below).

!

! f107d,f107a,hpower,e1,e2,h1,h2,alfac,ec,alfad,ed

!

if (.not.allocated(f107d)) allocate(f107d(ntime))

if (.not.allocated(f107a)) allocate(f107a(ntime))

if (.not.allocated(hpower)) allocate(hpower(ntime))

if (.not.allocated(e1)) allocate(e1(ntime))

if (.not.allocated(e2)) allocate(e2(ntime))

if (.not.allocated(h1)) allocate(h1(ntime))

if (.not.allocated(h2)) allocate(h2(ntime))

if (.not.allocated(alfac)) allocate(alfac(ntime))

if (.not.allocated(alfad)) allocate(alfad(ntime))

if (.not.allocated(ec)) allocate(ec(ntime))

if (.not.allocated(ed)) allocate(ed(ntime))

istat = nf90_inq_varid(ncid,'f107d',id)

istat = nf90_get_var (ncid,id,f107d)

istat = nf90_inq_varid(ncid,'f107a',id)

istat = nf90_get_var (ncid,id,f107a)

istat = nf90_inq_varid(ncid,'hpower',id)

istat = nf90_get_var (ncid,id,hpower)

istat = nf90_inq_varid(ncid,'e1',id)

istat = nf90_get_var (ncid,id,e1)

istat = nf90_inq_varid(ncid,'e2',id)

istat = nf90_get_var (ncid,id,e2)

istat = nf90_inq_varid(ncid,'h1',id)

istat = nf90_get_var (ncid,id,h1)

istat = nf90_inq_varid(ncid,'h2',id)

istat = nf90_get_var (ncid,id,h2)

istat = nf90_inq_varid(ncid,'alfac',id)

istat = nf90_get_var (ncid,id,alfac)

istat = nf90_inq_varid(ncid,'alfad',id)

istat = nf90_get_var (ncid,id,alfad)

istat = nf90_inq_varid(ncid,'ec',id)

istat = nf90_get_var (ncid,id,ec)

istat = nf90_inq_varid(ncid,'ed',id)

istat = nf90_get_var (ncid,id,ed)

! write(6,"('read_tgcm_coords: ntime=',i4,' f107d=',/,(8f9.2))") ntime,f107d

! write(6,"('read_tgcm_coords: ntime=',i4,' f107a=',/,(8f9.2))") ntime,f107a

! write(6,"('read_tgcm_coords: ntime=',i4,' hpower= ',/,(8f9.2))") ntime,hpower

! write(6,"('read_tgcm_coords: ntime=',i4,' e1= ',/,(8f9.2))") ntime,e1

! write(6,"('read_tgcm_coords: ntime=',i4,' e2= ',/,(8f9.2))") ntime,e2

! write(6,"('read_tgcm_coords: ntime=',i4,' h1= ',/,(8f9.2))") ntime,h1

! write(6,"('read_tgcm_coords: ntime=',i4,' h2= ',/,(8f9.2))") ntime,h2

! write(6,"('read_tgcm_coords: ntime=',i4,' alfac=',/,(8f9.2))") ntime,alfac

! write(6,"('read_tgcm_coords: ntime=',i4,' alfad=',/,(8f9.2))") ntime,alfad

! write(6,"('read_tgcm_coords: ntime=',i4,' ec= ',/,(8f9.2))") ntime,ec

! write(6,"('read_tgcm_coords: ntime=',i4,' ed= ',/,(8f9.2))") ntime,ed

istat = nf90_close(ncid)

end subroutine read_tgcm_coords

!-----------------------------------------------------------------------

subroutine handle_ncerr(istat,msg,ifatal)

!

! Handle a netcdf lib error:

!

integer,intent(in) :: istat,ifatal

character(len=*),intent(in) :: msg

!

write(6,"(/72('-'))")

write(6,"('>>> Error from netcdf library:')")

write(6,"(a)") trim(msg)

write(6,"('istat=',i5)") istat

write(6,"(a)") nf90_strerror(istat)

write(6,"(72('-')/)")

if (ifatal > 0) stop('Fatal netcdf error')

end subroutine handle_ncerr

!-----------------------------------------------------------------------

end module readtgcm