Line data Source code
1 : !===============================================================================
2 : module calc_pressure
3 :
4 : implicit none
5 :
6 : public :: init_pressure, & ! Procedure(s)
7 : calculate_thvm
8 :
9 : private ! default scope
10 :
11 : contains
12 :
13 : !=============================================================================
14 0 : subroutine init_pressure( gr, thvm, p_sfc, &
15 0 : p_in_Pa, exner, p_in_Pa_zm, exner_zm )
16 :
17 : ! Description:
18 : ! Calculates the initial pressure according to the hydrostatic
19 : ! approximation. Combining the moist equation of state and the hydrostatic
20 : ! approximation, the change of pressure with respect to height can be
21 : ! calculated based on theta_v, such that:
22 : !
23 : ! dp/dz = - p * grav / ( Rd * theta_v * exner );
24 : !
25 : ! where exner = ( p / p0 )^(Rd/Cp);
26 : !
27 : ! and where p0 is a reference pressure of 100000 Pa.
28 : !
29 : ! The integral equation is set up to integrate over p on the left-hand side
30 : ! and integrate over z on the right-hand side. The equation is:
31 : !
32 : ! INT(p1:p2) p^(Rd/Cp-1) dp
33 : ! = - p0^(Rd/Cp) * ( grav / Rd ) * INT(z1:z2) (1/thvm) dz.
34 : !
35 : ! The value of mean theta_v (thvm) is calculated at each thermodynamic grid
36 : ! level, and linear interpolation is used in the integral equation for all
37 : ! altitude in-between successive thermodynamic levels, such that:
38 : !
39 : ! thvm(z) = ( ( thvm2 - thvm1 ) / ( z2 - z1 ) ) * ( z - z1 ) + thvm1.
40 : !
41 : ! The integrals are solved, and the results for pressure can be rewritten
42 : ! in terms of exner, such that:
43 : !
44 : ! exner2 - exner1
45 : ! | - ( grav / Cp )
46 : ! | * ( ( z2 - z1 ) / ( thvm2 - thvm1 ) ) * ln( thvm2 / thvm1 );
47 : ! = | where thvm2 /= thvm1;
48 : ! |
49 : ! | - ( grav / Cp ) * ( z2 - z1 ) / thvm; where thvm2 = thvm1 (= thvm).
50 : !
51 : ! The value of pressure (exner) can be calculated using the above equation
52 : ! at all vertical levels once the value of pressure is known at one level.
53 : ! Since the surface pressure is known at the initial time, that allows
54 : ! pressure to be calculated for the rest of the vertical profile.
55 :
56 : ! References:
57 : !-----------------------------------------------------------------------
58 :
59 : use grid_class, only: &
60 : grid, & ! Type
61 : zt2zm ! Procedure(s)
62 :
63 : use constants_clubb, only: &
64 : one, & ! 1
65 : Cp, & ! Specific heat of dry air [J/(kg K)]
66 : kappa, & ! Rd/Cp [-]
67 : p0, & ! Reference pressure of 100000 Pa [Pa]
68 : grav, & ! Acceleration of gravity (9.81 m/s^2) [m/s^2]
69 : zero_threshold
70 :
71 : use clubb_precision, only: &
72 : core_rknd ! Variable(s)
73 :
74 : implicit none
75 :
76 : type (grid), target, intent(in) :: gr
77 :
78 : ! Input Variables
79 : real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
80 : thvm ! Mean theta_v (thermodynamic levels) [K]
81 :
82 : real( kind = core_rknd ), intent(in) :: &
83 : p_sfc ! Surface pressure [Pa]
84 :
85 : ! Output Variables
86 : real( kind = core_rknd ), dimension(gr%nz), intent(out) :: &
87 : p_in_Pa, & ! Pressure (thermodynamic levels) [Pa]
88 : exner, & ! Exner function (thermodynamic levels) [-]
89 : p_in_Pa_zm, & ! Pressure on momentum levels [Pa]
90 : exner_zm ! Exner function on momentum levels [-]
91 :
92 : ! Local Variables
93 : real( kind = core_rknd ), dimension(gr%nz) :: &
94 0 : thvm_zm ! Mean theta_v interpolated to momentum grid levels [K]
95 :
96 : real( kind = core_rknd ), parameter :: &
97 : g_ov_Cp = grav / Cp ! g / Cp [K/m]
98 :
99 : integer :: k ! Vertical level index
100 :
101 :
102 : ! The pressure (and exner) at the lowest momentum level is the surface
103 : ! pressure (and exner based on the surface pressure).
104 0 : p_in_Pa_zm(1) = p_sfc
105 0 : exner_zm(1) = ( p_sfc / p0 )**kappa
106 :
107 : ! Set the pressure (and exner) at the lowest thermodynamic level, which is
108 : ! below the model lower boundary, to their values at the model lower
109 : ! boundary or surface.
110 0 : p_in_Pa(1) = p_in_Pa_zm(1)
111 0 : exner(1) = exner_zm(1)
112 :
113 : ! Interpolate theta_v to momentum levels.
114 0 : thvm_zm = zt2zm( gr, thvm, zero_threshold )
115 :
116 : ! Calculate exner at all other thermodynamic and momentum grid levels.
117 : ! exner2
118 : ! = exner1
119 : ! | ( grav / Cp )
120 : ! | * ( ( z2 - z1 ) / ( thvm2 - thvm1 ) ) * ln( thvm2 / thvm1 );
121 : ! - | where thvm2 /= thvm1;
122 : ! |
123 : ! | ( grav / Cp ) * ( z2 - z1 ) / thvm; where thvm2 = thvm1 (= thvm).
124 :
125 : ! Calculate exner at thermodynamic level 2 (first thermodynamic grid level
126 : ! that is above the lower boundary).
127 0 : if ( abs( thvm(2) - thvm_zm(1) ) > epsilon( thvm ) * thvm(2) ) then
128 :
129 : exner(2) &
130 : = exner_zm(1) &
131 0 : - g_ov_Cp * ( gr%zt(1,2) - gr%zm(1,1) ) / ( thvm(2) - thvm_zm(1) ) &
132 0 : * log( thvm(2) / thvm_zm(1) )
133 :
134 : else ! thvm(2) = thvm_zm(1)
135 :
136 0 : exner(2) = exner_zm(1) - g_ov_Cp * ( gr%zt(1,2) - gr%zm(1,1) ) / thvm(2)
137 :
138 : endif
139 :
140 : ! Calculate pressure on thermodynamic levels.
141 0 : p_in_Pa(2) = p0 * exner(2)**(one/kappa)
142 :
143 : ! Loop over all other thermodynamic vertical grid levels.
144 0 : do k = 3, gr%nz, 1
145 :
146 : ! Calculate exner on thermodynamic levels.
147 0 : if ( abs( thvm(k) - thvm(k-1) ) > epsilon( thvm ) * thvm(k) ) then
148 :
149 : exner(k) &
150 : = exner(k-1) &
151 0 : - g_ov_Cp * ( gr%zt(1,k) - gr%zt(1,k-1) ) / ( thvm(k) - thvm(k-1) ) &
152 0 : * log( thvm(k) / thvm(k-1) )
153 :
154 : else ! thvm(k+1) = thvm(k)
155 :
156 0 : exner(k) = exner(k-1) - g_ov_Cp * ( gr%zt(1,k) - gr%zt(1,k-1) ) / thvm(k)
157 :
158 : endif
159 :
160 : ! Calculate pressure on thermodynamic levels.
161 0 : p_in_Pa(k) = p0 * exner(k)**(one/kappa)
162 :
163 : enddo ! k = 2, gr%nz, 1
164 :
165 : ! Loop over all momentum grid levels.
166 0 : do k = 2, gr%nz, 1
167 :
168 : ! Calculate exner on momentum levels.
169 0 : if ( abs( thvm(k) - thvm_zm(k) ) > epsilon( thvm ) * thvm(k) ) then
170 :
171 : exner_zm(k) &
172 : = exner(k) &
173 0 : - g_ov_Cp * ( gr%zm(1,k) - gr%zt(1,k) ) / ( thvm_zm(k) - thvm(k) ) &
174 0 : * log( thvm_zm(k) / thvm(k) )
175 :
176 : else ! thvm(k) = thvm_zm(k)
177 :
178 : exner_zm(k) &
179 0 : = exner(k) - g_ov_Cp * ( gr%zm(1,k) - gr%zt(1,k) ) / thvm_zm(k)
180 :
181 : endif
182 :
183 : ! Calculate pressure on momentum levels.
184 0 : p_in_Pa_zm(k) = p0 * exner_zm(k)**(one/kappa)
185 :
186 : enddo ! k = 2, gr%nz, 1
187 :
188 :
189 0 : return
190 :
191 : end subroutine init_pressure
192 :
193 : !=============================================================================
194 8935056 : subroutine calculate_thvm( nz, ngrdcol, &
195 8935056 : thlm, rtm, rcm, exner, thv_ds_zt, &
196 8935056 : thvm)
197 : ! Description:
198 : ! Calculates mean theta_v based on a linearized approximation to the theta_v
199 : ! equation. This equation also includes liquid water loading.
200 :
201 : ! References:
202 : !-----------------------------------------------------------------------
203 :
204 : use constants_clubb, only: &
205 : Lv, & ! Latent Heat of Vaporizaion [J/kg]
206 : Cp, & ! Specific Heat of Dry Air [J/(kg K)]
207 : ep1, & ! Rv/Rd - 1 [-]
208 : ep2 ! Rv/Rd [-]
209 :
210 : use clubb_precision, only: &
211 : core_rknd
212 :
213 : implicit none
214 :
215 : ! ---------------------------- Input Variables ----------------------------
216 : integer, intent(in) :: &
217 : nz, &
218 : ngrdcol
219 :
220 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
221 : thlm, & ! Mean theta_l (thermodynamic levels) [K]
222 : rtm, & ! Mean total water (thermodynamic levels) [kg/kg]
223 : rcm, & ! Mean cloud water (thermodynamic levels) [kg/kg]
224 : exner, & ! Exner function (thermodynamic levels) [-]
225 : thv_ds_zt ! Reference theta_v on thermodynamic levels [K]
226 :
227 : ! ---------------------------- Return Variable ----------------------------
228 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
229 : thvm ! Mean theta_v (thermodynamic levels) [K]
230 :
231 : ! ---------------------------- Local Variables ----------------------------
232 : integer :: i, k
233 :
234 : ! ---------------------------- Begin Code ----------------------------
235 :
236 : !$acc data copyin( thlm, rtm, rcm, exner, thv_ds_zt ) &
237 : !$acc copyout( thvm )
238 :
239 : ! Calculate mean theta_v
240 : !$acc parallel loop gang vector collapse(2) default(present)
241 768414816 : do k = 1, nz
242 12690480816 : do i = 1, ngrdcol
243 23844132000 : thvm(i,k) = thlm(i,k) + ep1 * thv_ds_zt(i,k) * rtm(i,k) &
244 36525677760 : + ( Lv / ( Cp * exner(i,k) ) - ep2 * thv_ds_zt(i,k) ) * rcm(i,k)
245 : end do
246 : end do
247 : !$acc end parallel loop
248 :
249 : !$acc end data
250 :
251 8935056 : return
252 :
253 : end subroutine calculate_thvm
254 : !=============================================================================
255 :
256 : end module calc_pressure
|