forked from NCAR/GLOW
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrcolum.f90
164 lines (124 loc) · 4.55 KB
/
rcolum.f90
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
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
! Subroutine RCOLUM
! 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, 1988, 1991
! Stan Solomon, 2016: removed problematic extrapolation below lower
! 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
! vertical column density, and then uses a fit to the Chapman Grazing
! Incidence Integral [Smith and Smith, JGR 77, 3592, 1972] to calculate
! the slant column density. If CHI is less than 90 degrees, column
! densities are calculated directly; if CHI is greater than 90 degrees
! the column density at grazing height for 90 degrees is calculated and
! doubled, and the column density above ZZ(J) is subtracted. If the
! grazing height is lower than the bottom of the atmosphere supplied,
! column densities are set to 'infinity', i.e., 1.0e30.
subroutine rcolum (chi, zz, zmaj, tn, zcol, zvcd, jmax, nmaj)
implicit none
integer,intent(in) :: jmax, nmaj
real,intent(in) :: chi, zz(jmax), zmaj(nmaj,jmax), tn(jmax)
real,intent(out) :: zcol(nmaj,jmax), zvcd(nmaj,jmax)
integer,parameter :: nm=3
real,parameter :: pi=3.1415926536
real,parameter :: re=6.37e8
integer :: i, j, k
real :: zcg(nm), ghrg, ghz, tng
real,external :: chap
call vcd (zz, zmaj, zvcd, jmax, nmaj)
if (chi >= 2.) then
do i=1,nmaj
do j=1,jmax
zcol(i,j) = 1.0e30
enddo
enddo
return
endif
if (chi <= pi/2.) then
do i=1,nmaj
do j=1,jmax
zcol(i,j) = zvcd(i,j) * chap(chi,zz(j),tn(j),i)
enddo
enddo
else
do j=1,jmax
ghrg=(re+zz(j))*sin(chi)
ghz=ghrg-re
if (ghz <= zz(1)) then
do i=1,nmaj
zcol(i,j) = 1.0e30
enddo
else
do k=1,j-1
if (zz(k) <= ghz .and. zz(k+1) > ghz) then
tng = tn(k)+(tn(k+1)-tn(k))*(ghz-zz(k))/(zz(k+1)-zz(k))
do i=1,nmaj
zcg(i) = zvcd(i,k) * (zvcd(i,k+1) / zvcd(i,k)) ** &
((ghz-zz(k)) / (zz(k+1)-zz(k)))
enddo
endif
enddo
do i=1,nmaj
zcol(i,j) = 2. * zcg(i) * chap(pi/2.,ghz,tng,i) &
- zvcd(i,j) * chap(pi-chi,zz(j),tn(j),i)
enddo
endif
enddo
endif
return
end subroutine rcolum
!----------------------------------------------------------------------
real function chap (chi, z, t, i)
implicit none
real,intent(in) :: chi, z, t
integer,intent(in) :: i
integer,parameter :: nmaj=3
real,parameter :: pi=3.1415926536
real,parameter :: re=6.37e8
real,parameter :: g=978.1
real :: am(nmaj), gr, hn, hg, hf, sqhf
real,external :: sperfc
data am/16., 32., 28./
gr=g*(re/(re+z))**2
hn=1.38e-16*t/(am(i)*1.662e-24*gr)
hg=(re+z)/hn
hf=0.5*hg*(cos(chi)**2)
sqhf=sqrt(hf)
chap=sqrt(0.5*pi*hg)*sperfc(sqhf)
return
end function chap
!----------------------------------------------------------------------
real function sperfc(dummy)
implicit none
real,intent(in) :: dummy
if (dummy <= 8.) then
sperfc = (1.0606963+0.55643831*dummy) / &
(1.0619896+1.7245609*dummy+dummy*dummy)
else
sperfc=0.56498823/(0.06651874+dummy)
endif
return
end function sperfc
!----------------------------------------------------------------------
subroutine vcd(zz,zmaj,zvcd,jmax,nmaj)
implicit none
integer,intent(in) :: jmax,nmaj
real,intent(in) :: zz(jmax), zmaj(nmaj,jmax)
real,intent(out) :: zvcd(nmaj,jmax)
integer :: i, j
real :: rat
do i=1,nmaj
zvcd(i,jmax) = zmaj(i,jmax) &
* (zz(jmax)-zz(jmax-1)) &
/ alog(zmaj(i,jmax-1)/zmaj(i,jmax))
do j=jmax-1,1,-1
rat = zmaj(i,j+1) / zmaj(i,j)
zvcd(i,j) = zvcd(i,j+1)+zmaj(i,j)*(zz(j)-zz(j+1))/alog(rat)*(1.-rat)
enddo
enddo
return
end subroutine vcd