-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathatmos_ops_mod.f90
156 lines (88 loc) · 4.1 KB
/
atmos_ops_mod.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
module atmos_ops
use sizes
use common_arrays
use phys_const
use define_types
implicit none
contains
subroutine set_pressure_scale
integer:: i,ipatch
do ipatch= 1, npatch
patch(ipatch)%atm%index = (/ (i, i = 1, nlayers) /)
patch(ipatch)%atm%press= press
patch(ipatch)%atm%logP = log10(patch(ipatch)%atm%press)
end do
end subroutine set_pressure_scale
subroutine layer_thickness(grav)
implicit none
!declare counters here to ensure they don't clash with counters from calling
! program
integer:: i
double precision :: p1, p2, R_spec
real, intent(IN) :: grav
! The values in the pressure grid are geometric means for the layer.
! Assume these are taken as sqrt(Pbot * Ptop)
! So (see derivation in notes)
! log(P1) = (3/2)*log(Pbar1) - (1/2)*log(Pbar2)
! log(P2) = 1/2(log (Pbar1) + log(Pbar2))
! we're using the hypsometric equation for this..
! careful! this requires SPECIFIC GAS CONSTANT
do i = 1, nlayers
R_spec = K_boltz / (patch(1)%atm(i)%mu * amu)
if (i .eq. nlayers) then
p1 = exp((0.5)*(log(patch(1)%atm(i-1)%press * patch(1)%atm(i)%press)))
p2 = (patch(1)%atm(i)%press**2) / p1
patch(1)%atm(i)%dp = p2 - p1
! no need to change P units to Pa, as we want the ratio
patch(1)%atm(i)%dz = abs((R_spec * patch(1)%atm(i)%temp / grav) * log(p1 /p2))
else
p1 = exp(((1.5)*log(patch(1)%atm(i)%press)) - &
((0.5)*log(patch(1)%atm(i+1)%press)))
p2 = exp((0.5)*(log(patch(1)%atm(i)%press * patch(1)%atm(i+1)%press)))
patch(1)%atm(i)%dp = p2 - p1
! no need to change P units to Pa, as we want the ratio
patch(1)%atm(i)%dz = abs((R_spec * patch(1)%atm(i)%temp / grav) * log(p1 /p2))
end if
end do
! last layer needs special treatment
! is at the bottom of the atmosphere, so probably doesn't matter too much
! make it a mirror of layer above
end subroutine layer_thickness
subroutine set_temp_levels(leveltemp)
! DISORT wants the temperature at levels, which we here retrieve from the
! Tbars that are in our atmospheric layers.
! This is kept in atmos_ops for consistency a
integer :: i
double precision :: p1,p2
double precision,INTENT(OUT):: leveltemp(0:nlayers)
! we'll interpolate the level temps using Pbar, Tbar, and dP
do i = 1, nlayers
if (i .eq. 1) then
p1 = exp(((1.5)*log(patch(1)%atm(i)%press)) - &
((0.5)*log(patch(1)%atm(i+1)%press)))
p2 = exp((0.5)*(log(patch(1)%atm(i)%press * patch(1)%atm(i+1)%press)))
leveltemp(i-1) = patch(1)%atm(i)%temp + &
(((patch(1)%atm(i+1)%temp - patch(1)%atm(i)%temp) / &
(patch(1)%atm(i+1)%logP - patch(1)%atm(i)%logP)) * &
(log10(p1) - patch(1)%atm(i)%logP))
leveltemp(i) = patch(1)%atm(i)%temp + &
(((patch(1)%atm(i+1)%temp - patch(1)%atm(i)%temp) / &
(patch(1)%atm(i+1)%logP - patch(1)%atm(i)%logP)) * &
(log10(p2) - patch(1)%atm(i)%logP))
elseif (i .eq. nlayers) then
p1 = exp((0.5)*(log(patch(1)%atm(i-1)%press * patch(1)%atm(i)%press)))
p2 = (patch(1)%atm(i)%press**2) / p1
leveltemp(i) = patch(1)%atm(i)%temp + &
(((patch(1)%atm(i)%temp - patch(1)%atm(i-1)%temp) / &
(patch(1)%atm(i)%logP - patch(1)%atm(i-1)%logP)) * &
(log10(p2) - patch(1)%atm(i)%logP))
else
p2 = exp((0.5)*(log(patch(1)%atm(i)%press * patch(1)%atm(i+1)%press)))
leveltemp(i) = patch(1)%atm(i)%temp + &
(((patch(1)%atm(i+1)%temp - patch(1)%atm(i)%temp) / &
(patch(1)%atm(i+1)%logP - patch(1)%atm(i)%logP)) * &
(log10(p2) - patch(1)%atm(i)%logP))
end if
end do
end subroutine set_temp_levels
end module atmos_ops