1 | MODULE dcmip2016_supercell_mod |
---|
2 | |
---|
3 | !======================================================================= |
---|
4 | ! |
---|
5 | ! Date: April 22, 2016 |
---|
6 | ! |
---|
7 | ! Functions for setting up idealized initial conditions for the |
---|
8 | ! Klemp et al. supercell test. Before sampling the result, |
---|
9 | ! supercell_test_init() must be called. |
---|
10 | ! |
---|
11 | ! SUBROUTINE supercell_test( |
---|
12 | ! lon,lat,p,z,zcoords,u,v,t,thetav,ps,rho,q,pert) |
---|
13 | ! |
---|
14 | ! Given a point specified by: |
---|
15 | ! lon longitude (radians) |
---|
16 | ! lat latitude (radians) |
---|
17 | ! p/z pressure (Pa) / height (m) |
---|
18 | ! zcoords 1 if z is specified, 0 if p is specified |
---|
19 | ! pert 1 if thermal perturbation included, 0 if not |
---|
20 | ! |
---|
21 | ! the functions will return: |
---|
22 | ! p pressure if z is specified (Pa) |
---|
23 | ! z geopotential height if p is specified (m) |
---|
24 | ! u zonal wind (m s^-1) |
---|
25 | ! v meridional wind (m s^-1) |
---|
26 | ! t temperature (K) |
---|
27 | ! thetav virtual potential temperature (K) |
---|
28 | ! ps surface pressure (Pa) |
---|
29 | ! rho density (kj m^-3) |
---|
30 | ! q water vapor mixing ratio (kg/kg) |
---|
31 | ! |
---|
32 | ! Author: Paul Ullrich |
---|
33 | ! University of California, Davis |
---|
34 | ! Email: paullrich@ucdavis.edu |
---|
35 | ! |
---|
36 | ! Based on a code by Joseph Klemp |
---|
37 | ! (National Center for Atmospheric Research) |
---|
38 | ! |
---|
39 | !======================================================================= |
---|
40 | |
---|
41 | IMPLICIT NONE |
---|
42 | |
---|
43 | !======================================================================= |
---|
44 | ! Physical constants |
---|
45 | !======================================================================= |
---|
46 | |
---|
47 | REAL(8), PARAMETER :: & |
---|
48 | a = 6371220.0d0, & ! Reference Earth's Radius (m) |
---|
49 | Rd = 287.0d0, & ! Ideal gas const dry air (J kg^-1 K^1) |
---|
50 | g = 9.80616d0, & ! Gravity (m s^2) |
---|
51 | cp = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1) |
---|
52 | Lvap = 2.5d6, & ! Latent heat of vaporization of water |
---|
53 | Rvap = 461.5d0, & ! Ideal gas constnat for water vapor |
---|
54 | Mvap = 0.608d0, & ! Ratio of molar mass of dry air/water |
---|
55 | pi = 3.14159265358979d0, & ! pi |
---|
56 | p0 = 100000.0d0, & ! surface pressure (Pa) |
---|
57 | kappa = 2.d0/7.d0, & ! Ratio of Rd to cp |
---|
58 | omega = 7.29212d-5, & ! Reference rotation rate of the Earth (s^-1) |
---|
59 | deg2rad = pi/180.d0 ! Conversion factor of degrees to radians |
---|
60 | |
---|
61 | !======================================================================= |
---|
62 | ! Test case parameters |
---|
63 | !======================================================================= |
---|
64 | INTEGER(4), PARAMETER :: & |
---|
65 | nz = 30 , & ! number of vertical levels in init |
---|
66 | nphi = 16 ! number of meridional points in init |
---|
67 | |
---|
68 | REAL(8), PARAMETER :: & |
---|
69 | z1 = 0.0d0 , & ! lower sample altitude |
---|
70 | z2 = 50000.0d0 ! upper sample altitude |
---|
71 | |
---|
72 | REAL(8), PARAMETER :: & |
---|
73 | X = 120.d0 , & ! Earth reduction factor |
---|
74 | theta0 = 300.d0 , & ! theta at the equatorial surface |
---|
75 | theta_tr = 343.d0 , & ! theta at the tropopause |
---|
76 | z_tr = 12000.d0 , & ! altitude at the tropopause |
---|
77 | T_tr = 213.d0 , & ! temperature at the tropopause |
---|
78 | pseq = 100000.0d0 ! surface pressure at equator (Pa) |
---|
79 | !pseq = 95690.0d0 ! surface pressure at equator (Pa) |
---|
80 | |
---|
81 | REAL(8), PARAMETER :: & |
---|
82 | us = 30.d0 , & ! maximum zonal wind velocity |
---|
83 | uc = 15.d0 , & ! coordinate reference velocity |
---|
84 | zs = 5000.d0 , & ! lower altitude of maximum velocity |
---|
85 | zt = 1000.d0 ! transition distance of velocity |
---|
86 | |
---|
87 | REAL(8), PARAMETER :: & |
---|
88 | pert_dtheta = 3.d0 , & ! perturbation magnitude |
---|
89 | pert_lonc = 0.d0 , & ! perturbation longitude |
---|
90 | pert_latc = 0.d0 , & ! perturbation latitude |
---|
91 | pert_rh = 10000.d0 * X , & ! perturbation horiz. halfwidth |
---|
92 | pert_zc = 1500.d0 , & ! perturbation center altitude |
---|
93 | pert_rz = 1500.d0 ! perturbation vert. halfwidth |
---|
94 | |
---|
95 | !----------------------------------------------------------------------- |
---|
96 | ! Coefficients computed from initialization |
---|
97 | !----------------------------------------------------------------------- |
---|
98 | INTEGER(4) :: initialized = 0 |
---|
99 | |
---|
100 | REAL(8), DIMENSION(nphi) :: phicoord |
---|
101 | REAL(8), DIMENSION(nz) :: zcoord |
---|
102 | REAL(8), DIMENSION(nphi,nz) :: thetavyz |
---|
103 | REAL(8), DIMENSION(nphi,nz) :: exneryz |
---|
104 | REAL(8), DIMENSION(nz) :: qveq |
---|
105 | |
---|
106 | CONTAINS |
---|
107 | |
---|
108 | !======================================================================= |
---|
109 | ! Generate the supercell initial conditions |
---|
110 | !======================================================================= |
---|
111 | SUBROUTINE supercell_init() & |
---|
112 | BIND(c, name = "supercell_init") |
---|
113 | |
---|
114 | IMPLICIT NONE |
---|
115 | |
---|
116 | ! d/dphi and int(dphi) operators |
---|
117 | REAL(8), DIMENSION(nphi,nphi) :: ddphi, intphi |
---|
118 | |
---|
119 | ! d/dz and int(dz) operators |
---|
120 | REAL(8), DIMENSION(nz, nz) :: ddz, intz |
---|
121 | |
---|
122 | ! Buffer matrices for computing SVD of d/dphi operator |
---|
123 | REAL(8), DIMENSION(nphi,nphi) :: ddphibak |
---|
124 | REAL(8), DIMENSION(nphi,nphi) :: svdpu, svdpvt |
---|
125 | REAL(8), DIMENSION(nphi) :: svdps |
---|
126 | REAL(8), DIMENSION(5*nphi) :: pwork |
---|
127 | |
---|
128 | ! Buffer matrices for computing SVD of d/dz operator |
---|
129 | REAL(8), DIMENSION(nz, nz) :: ddzbak |
---|
130 | REAL(8), DIMENSION(nz, nz) :: svdzu, svdzvt |
---|
131 | REAL(8), DIMENSION(nz) :: svdzs |
---|
132 | REAL(8), DIMENSION(5*nz) :: zwork |
---|
133 | |
---|
134 | ! Buffer data for calculation of SVD |
---|
135 | INTEGER(4) :: lwork, info |
---|
136 | |
---|
137 | ! Sampled values of ueq**2 and d/dz(ueq**2) |
---|
138 | REAL(8), DIMENSION(nphi, nz) :: ueq2, dueq2 |
---|
139 | |
---|
140 | ! Buffer matrices for iteration |
---|
141 | REAL(8), DIMENSION(nphi, nz) :: phicoordmat, dztheta, rhs, irhs |
---|
142 | |
---|
143 | ! Buffer for sampled potential temperature at equator |
---|
144 | REAL(8), DIMENSION(nz) :: thetaeq |
---|
145 | |
---|
146 | ! Buffer for computed equatorial Exner pressure and relative humidity |
---|
147 | REAL(8), DIMENSION(nz) :: exnereq, H |
---|
148 | |
---|
149 | ! Variables for calculation of equatorial profile |
---|
150 | REAL(8) :: exnereqs, p, T, qvs, qv |
---|
151 | |
---|
152 | ! Error metric |
---|
153 | REAL(8) :: err |
---|
154 | |
---|
155 | ! Loop indices |
---|
156 | INTEGER(4) :: i, k, iter |
---|
157 | |
---|
158 | ! Chebyshev nodes in the phi direction |
---|
159 | do i = 1, nphi |
---|
160 | phicoord(i) = - cos(dble(i-1) * pi / dble(nphi-1)) |
---|
161 | phicoord(i) = 0.25d0 * pi * (phicoord(i) + 1.0d0) |
---|
162 | end do |
---|
163 | |
---|
164 | ! Matrix of phis |
---|
165 | do k = 1, nz |
---|
166 | phicoordmat(:,k) = phicoord |
---|
167 | end do |
---|
168 | |
---|
169 | ! Chebyshev nodes in the z direction |
---|
170 | do k = 1, nz |
---|
171 | zcoord(k) = - cos(dble(k-1) * pi / dble(nz-1)) |
---|
172 | zcoord(k) = z1 + 0.5d0*(z2-z1)*(zcoord(k)+1.0d0) |
---|
173 | end do |
---|
174 | |
---|
175 | ! Compute the d/dphi operator |
---|
176 | do i = 1, nphi |
---|
177 | call diff_lagrangian_polynomial_coeffs( & |
---|
178 | nphi, phicoord, ddphi(:,i), phicoord(i)) |
---|
179 | end do |
---|
180 | |
---|
181 | ! Zero derivative at pole |
---|
182 | ddphi(:,nphi) = 0.0d0 |
---|
183 | |
---|
184 | ! Compute the d/dz operator |
---|
185 | do k = 1, nz |
---|
186 | call diff_lagrangian_polynomial_coeffs( & |
---|
187 | nz, zcoord, ddz(:,k), zcoord(k)) |
---|
188 | end do |
---|
189 | |
---|
190 | ! Compute the int(dphi) operator via pseudoinverse |
---|
191 | lwork = 5*nphi |
---|
192 | |
---|
193 | ddphibak = ddphi |
---|
194 | call DGESVD('A', 'A', & |
---|
195 | nphi, nphi, ddphibak, nphi, & |
---|
196 | svdps, svdpu, nphi, svdpvt, nphi, & |
---|
197 | pwork, lwork, info) |
---|
198 | |
---|
199 | if (info .ne. 0) then |
---|
200 | write(*,*) 'Unable to compute SVD of d/dphi matrix' |
---|
201 | stop |
---|
202 | end if |
---|
203 | |
---|
204 | do i = 1, nphi |
---|
205 | if (abs(svdps(i)) .le. 1.0d-12) then |
---|
206 | call DSCAL(nphi, 0.0d0, svdpu(1,i), 1) |
---|
207 | else |
---|
208 | call DSCAL(nphi, 1.0d0 / svdps(i), svdpu(1,i), 1) |
---|
209 | end if |
---|
210 | end do |
---|
211 | call DGEMM('T', 'T', & |
---|
212 | nphi, nphi, nphi, 1.0d0, svdpvt, nphi, svdpu, nphi, 0.0d0, & |
---|
213 | intphi, nphi) |
---|
214 | |
---|
215 | ! Compute the int(dz) operator via pseudoinverse |
---|
216 | lwork = 5*nz |
---|
217 | |
---|
218 | ddzbak = ddz |
---|
219 | call DGESVD('A', 'A', & |
---|
220 | nz, nz, ddzbak, nz, & |
---|
221 | svdzs, svdzu, nz, svdzvt, nz, & |
---|
222 | zwork, lwork, info) |
---|
223 | |
---|
224 | if (info .ne. 0) then |
---|
225 | write(*,*) 'Unable to compute SVD of d/dz matrix' |
---|
226 | stop |
---|
227 | end if |
---|
228 | |
---|
229 | do i = 1, nz |
---|
230 | if (abs(svdzs(i)) .le. 1.0d-12) then |
---|
231 | call DSCAL(nz, 0.0d0, svdzu(1,i), 1) |
---|
232 | else |
---|
233 | call DSCAL(nz, 1.0d0 / svdzs(i), svdzu(1,i), 1) |
---|
234 | end if |
---|
235 | end do |
---|
236 | call DGEMM('T', 'T', & |
---|
237 | nz, nz, nz, 1.0d0, svdzvt, nz, svdzu, nz, 0.0d0, & |
---|
238 | intz, nz) |
---|
239 | |
---|
240 | ! Sample the equatorial velocity field and its derivative |
---|
241 | do k = 1, nz |
---|
242 | ueq2(1,k) = zonal_velocity(zcoord(k), 0.0d0) |
---|
243 | ueq2(1,k) = ueq2(1,k)**2 |
---|
244 | end do |
---|
245 | do k = 1, nz |
---|
246 | dueq2(1,k) = dot_product(ddz(:,k), ueq2(1,:)) |
---|
247 | end do |
---|
248 | do i = 2, nphi |
---|
249 | ueq2(i,:) = ueq2(1,:) |
---|
250 | dueq2(i,:) = dueq2(1,:) |
---|
251 | end do |
---|
252 | |
---|
253 | ! Initialize potential temperature at equator |
---|
254 | do k = 1, nz |
---|
255 | thetaeq(k) = equator_theta(zcoord(k)) |
---|
256 | H(k) = equator_relative_humidity(zcoord(k)) |
---|
257 | end do |
---|
258 | thetavyz(1,:) = thetaeq |
---|
259 | |
---|
260 | ! Exner pressure at the equatorial surface |
---|
261 | exnereqs = (pseq / p0)**(Rd/cp) |
---|
262 | |
---|
263 | ! Iterate on equatorial profile |
---|
264 | do iter = 1, 12 |
---|
265 | |
---|
266 | ! Calculate Exner pressure in equatorial column (p0 at surface) |
---|
267 | rhs(1,:) = - g / cp / thetavyz(1,:) |
---|
268 | do k = 1, nz |
---|
269 | exnereq(k) = dot_product(intz(:,k), rhs(1,:)) |
---|
270 | end do |
---|
271 | do k = 2, nz |
---|
272 | exnereq(k) = exnereq(k) + (exnereqs - exnereq(1)) |
---|
273 | end do |
---|
274 | exnereq(1) = exnereqs |
---|
275 | |
---|
276 | ! Calculate new pressure and temperature |
---|
277 | do k = 1, nz |
---|
278 | p = p0 * exnereq(k)**(cp/Rd) |
---|
279 | T = thetaeq(k) * exnereq(k) |
---|
280 | |
---|
281 | qvs = saturation_mixing_ratio(p, T) |
---|
282 | qveq(k) = qvs * H(k) |
---|
283 | |
---|
284 | thetavyz(1,k) = thetaeq(k) * (1.d0 + 0.61d0 * qveq(k)) |
---|
285 | end do |
---|
286 | end do |
---|
287 | |
---|
288 | !do k = 1, nz |
---|
289 | ! write(*,*) exnereq(k) * thetaeq(k) |
---|
290 | !end do |
---|
291 | |
---|
292 | ! Iterate on remainder of domain |
---|
293 | do iter = 1, 12 |
---|
294 | |
---|
295 | ! Compute d/dz(theta) |
---|
296 | do i = 1, nphi |
---|
297 | do k = 1, nz |
---|
298 | dztheta(i,k) = dot_product(ddz(:,k), thetavyz(i,:)) |
---|
299 | end do |
---|
300 | end do |
---|
301 | |
---|
302 | ! Compute rhs |
---|
303 | rhs = sin(2.0d0*phicoordmat)/(2.0d0*g) & |
---|
304 | * (ueq2 * dztheta - thetavyz * dueq2) |
---|
305 | |
---|
306 | ! Integrate |
---|
307 | do k = 1, nz |
---|
308 | do i = 1, nphi |
---|
309 | irhs(i,k) = dot_product(intphi(:,i), rhs(:,k)) |
---|
310 | end do |
---|
311 | end do |
---|
312 | |
---|
313 | ! Apply boundary conditions (fixed Dirichlet condition at equator) |
---|
314 | do i = 2, nphi |
---|
315 | irhs(i,:) = irhs(i,:) + (thetavyz(1,:) - irhs(1,:)) |
---|
316 | end do |
---|
317 | irhs(1,:) = thetavyz(1,:) |
---|
318 | |
---|
319 | ! Compute difference after iteration |
---|
320 | !err = sum(irhs - thetavyz) |
---|
321 | !write(*,*) iter, err |
---|
322 | |
---|
323 | ! Update iteration |
---|
324 | thetavyz = irhs |
---|
325 | end do |
---|
326 | |
---|
327 | ! Calculate pressure through remainder of domain |
---|
328 | rhs = - ueq2 * sin(phicoordmat) * cos(phicoordmat) / cp / thetavyz |
---|
329 | |
---|
330 | do k = 1, nz |
---|
331 | do i = 1, nphi |
---|
332 | exneryz(i,k) = dot_product(intphi(:,i), rhs(:,k)) |
---|
333 | end do |
---|
334 | do i = 2, nphi |
---|
335 | exneryz(i,k) = exneryz(i,k) + (exnereq(k) - exneryz(1,k)) |
---|
336 | end do |
---|
337 | |
---|
338 | exneryz(1,k) = exnereq(k) |
---|
339 | end do |
---|
340 | |
---|
341 | ! Initialization successful |
---|
342 | initialized = 1 |
---|
343 | |
---|
344 | END SUBROUTINE supercell_init |
---|
345 | |
---|
346 | !----------------------------------------------------------------------- |
---|
347 | ! Evaluate the supercell initial conditions |
---|
348 | !----------------------------------------------------------------------- |
---|
349 | SUBROUTINE supercell_test(lon,lat,p,z,zcoords,u,v,t,thetav,ps,rho,q,pert) & |
---|
350 | BIND(c, name = "supercell_test") |
---|
351 | |
---|
352 | IMPLICIT NONE |
---|
353 | |
---|
354 | !------------------------------------------------ |
---|
355 | ! Input / output parameters |
---|
356 | !------------------------------------------------ |
---|
357 | REAL(8), INTENT(IN) :: & |
---|
358 | lon, & ! Longitude (radians) |
---|
359 | lat ! Latitude (radians) |
---|
360 | |
---|
361 | REAL(8), INTENT(INOUT) :: & |
---|
362 | p, & ! Pressure (Pa) |
---|
363 | z ! Altitude (m) |
---|
364 | |
---|
365 | INTEGER, INTENT(IN) :: zcoords ! 1 if z coordinates are specified |
---|
366 | ! 0 if p coordinates are specified |
---|
367 | |
---|
368 | REAL(8), INTENT(OUT) :: & |
---|
369 | u, & ! Zonal wind (m s^-1) |
---|
370 | v, & ! Meridional wind (m s^-1) |
---|
371 | t, & ! Temperature (K) |
---|
372 | thetav, & ! Virtual potential Temperature (K) |
---|
373 | ps, & ! Surface Pressure (Pa) |
---|
374 | rho, & ! density (kg m^-3) |
---|
375 | q ! water vapor mixing ratio (kg/kg) |
---|
376 | |
---|
377 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
378 | ! 0 if no perturbation should be included |
---|
379 | |
---|
380 | !------------------------------------------------ |
---|
381 | ! Local variables |
---|
382 | !------------------------------------------------ |
---|
383 | |
---|
384 | ! Absolute latitude |
---|
385 | REAL(8) :: nh_lat |
---|
386 | |
---|
387 | ! Check that we are initialized |
---|
388 | if (initialized .ne. 1) then |
---|
389 | write(*,*) 'supercell_init() has not been called' |
---|
390 | stop |
---|
391 | end if |
---|
392 | |
---|
393 | !------------------------------------------------ |
---|
394 | ! Begin sampling |
---|
395 | !------------------------------------------------ |
---|
396 | |
---|
397 | ! Northern hemisphere latitude |
---|
398 | if (lat .le. 0.0d0) then |
---|
399 | nh_lat = -lat |
---|
400 | else |
---|
401 | nh_lat = lat |
---|
402 | end if |
---|
403 | |
---|
404 | ! Sample surface pressure |
---|
405 | CALL supercell_z(lon, lat, 0.d0, ps, thetav, rho, q, pert) |
---|
406 | |
---|
407 | ! Calculate dependent variables |
---|
408 | if (zcoords .eq. 1) then |
---|
409 | CALL supercell_z(lon, lat, z, p, thetav, rho, q, pert) |
---|
410 | else |
---|
411 | CALL supercell_p(lon, lat, p, z, thetav, rho, q, pert) |
---|
412 | end if |
---|
413 | |
---|
414 | ! Sample the zonal velocity |
---|
415 | u = zonal_velocity(z, lat) |
---|
416 | |
---|
417 | ! Zero meridional velocity |
---|
418 | v = 0.d0 |
---|
419 | |
---|
420 | ! Temperature |
---|
421 | t = thetav / (1.d0 + 0.61d0 * q) * (p / p0)**(Rd/cp) |
---|
422 | |
---|
423 | END SUBROUTINE supercell_test |
---|
424 | |
---|
425 | !----------------------------------------------------------------------- |
---|
426 | ! Calculate pointwise pressure and temperature |
---|
427 | !----------------------------------------------------------------------- |
---|
428 | SUBROUTINE supercell_z(lon, lat, z, p, thetav, rho, q, pert) |
---|
429 | |
---|
430 | REAL(8), INTENT(IN) :: & |
---|
431 | lon, & ! Longitude (radians) |
---|
432 | lat, & ! Latitude (radians) |
---|
433 | z ! Altitude (m) |
---|
434 | |
---|
435 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
436 | ! 0 if no perturbation should be included |
---|
437 | |
---|
438 | ! Evaluated variables |
---|
439 | REAL(8), INTENT(OUT) :: p, thetav, rho, q |
---|
440 | |
---|
441 | ! Northern hemisphere latitude |
---|
442 | REAL(8) :: nh_lat |
---|
443 | |
---|
444 | ! Pointwise Exner pressure |
---|
445 | REAL(8) :: exner |
---|
446 | |
---|
447 | ! Assembled variable values in a column |
---|
448 | REAL(8), DIMENSION(nz) :: varcol |
---|
449 | |
---|
450 | ! Coefficients for computing a polynomial fit in each coordinate |
---|
451 | REAL(8), DIMENSION(nphi) :: fitphi |
---|
452 | REAL(8), DIMENSION(nz) :: fitz |
---|
453 | |
---|
454 | ! Loop indices |
---|
455 | INTEGER(4) :: k |
---|
456 | |
---|
457 | ! Northern hemisphere latitude |
---|
458 | if (lat .le. 0.0d0) then |
---|
459 | nh_lat = -lat |
---|
460 | else |
---|
461 | nh_lat = lat |
---|
462 | end if |
---|
463 | |
---|
464 | ! Perform fit |
---|
465 | CALL lagrangian_polynomial_coeffs(nz, zcoord, fitz, z) |
---|
466 | CALL lagrangian_polynomial_coeffs(nphi, phicoord, fitphi, nh_lat) |
---|
467 | |
---|
468 | ! Obtain exner pressure of background state |
---|
469 | do k = 1, nz |
---|
470 | varcol(k) = dot_product(fitphi, exneryz(:,k)) |
---|
471 | end do |
---|
472 | exner = dot_product(fitz, varcol) |
---|
473 | p = p0 * exner**(cp/Rd) |
---|
474 | |
---|
475 | ! Sample the initialized fit at this point for theta_v |
---|
476 | do k = 1, nz |
---|
477 | varcol(k) = dot_product(fitphi, thetavyz(:,k)) |
---|
478 | end do |
---|
479 | thetav = dot_product(fitz, varcol) |
---|
480 | |
---|
481 | ! Sample water vapor mixing ratio |
---|
482 | q = dot_product(fitz, qveq) |
---|
483 | |
---|
484 | ! Fixed density |
---|
485 | rho = p / (Rd * exner * thetav) |
---|
486 | |
---|
487 | ! Modified virtual potential temperature |
---|
488 | if (pert .ne. 0) then |
---|
489 | thetav = thetav & |
---|
490 | + thermal_perturbation(lon, lat, z) * (1.d0 + 0.61d0 * q) |
---|
491 | end if |
---|
492 | |
---|
493 | ! Updated pressure |
---|
494 | p = p0 * (rho * Rd * thetav / p0)**(cp/(cp-Rd)) |
---|
495 | |
---|
496 | END SUBROUTINE supercell_z |
---|
497 | |
---|
498 | !----------------------------------------------------------------------- |
---|
499 | ! Calculate pointwise z and temperature given pressure |
---|
500 | !----------------------------------------------------------------------- |
---|
501 | SUBROUTINE supercell_p(lon, lat, p, z, thetav, rho, q, pert) |
---|
502 | |
---|
503 | REAL(8), INTENT(IN) :: & |
---|
504 | lon, & ! Longitude (radians) |
---|
505 | lat, & ! Latitude (radians) |
---|
506 | p ! Pressure (Pa) |
---|
507 | |
---|
508 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
509 | ! 0 if no perturbation should be included |
---|
510 | |
---|
511 | ! Evaluated variables |
---|
512 | REAL(8), INTENT(OUT) :: z, thetav, rho, q |
---|
513 | |
---|
514 | ! Bounding interval and sampled values |
---|
515 | REAL(8) :: za, zb, zc, pa, pb, pc |
---|
516 | |
---|
517 | ! Iterate |
---|
518 | INTEGER(4) :: iter |
---|
519 | |
---|
520 | za = z |
---|
521 | zb = z2 |
---|
522 | |
---|
523 | CALL supercell_z(lon, lat, za, pa, thetav, rho, q, pert) |
---|
524 | CALL supercell_z(lon, lat, zb, pb, thetav, rho, q, pert) |
---|
525 | |
---|
526 | if (pa .lt. p) then |
---|
527 | write(*,*) 'Requested pressure out of range on bottom, adjust sample interval' |
---|
528 | write(*,*) pa, p |
---|
529 | stop |
---|
530 | end if |
---|
531 | if (pb .gt. p) then |
---|
532 | write(*,*) 'Requested pressure out of range on top, adjust sample interval' |
---|
533 | write(*,*) pb, p |
---|
534 | stop |
---|
535 | end if |
---|
536 | |
---|
537 | ! Iterate using fixed point method |
---|
538 | do iter = 1, 100 |
---|
539 | |
---|
540 | zc = (za * (pb - p) - zb * (pa - p)) / (pb - pa) |
---|
541 | |
---|
542 | CALL supercell_z(lon, lat, zc, pc, thetav, rho, q, pert) |
---|
543 | |
---|
544 | !write(*,*) pc |
---|
545 | |
---|
546 | if (abs((pc - p) / p) .lt. 1.d-10) then |
---|
547 | exit |
---|
548 | end if |
---|
549 | |
---|
550 | if (pc .gt. p) then |
---|
551 | za = zc |
---|
552 | pa = pc |
---|
553 | else |
---|
554 | zb = zc |
---|
555 | pb = pc |
---|
556 | end if |
---|
557 | end do |
---|
558 | |
---|
559 | if (iter .eq. 101) then |
---|
560 | write(*,*) 'Iteration failed to converge' |
---|
561 | stop |
---|
562 | end if |
---|
563 | |
---|
564 | z = zc |
---|
565 | |
---|
566 | END SUBROUTINE supercell_p |
---|
567 | |
---|
568 | !----------------------------------------------------------------------- |
---|
569 | ! Calculate pointwise z and temperature given pressure |
---|
570 | !----------------------------------------------------------------------- |
---|
571 | REAL(8) FUNCTION thermal_perturbation(lon, lat, z) |
---|
572 | |
---|
573 | REAL(8), INTENT(IN) :: & |
---|
574 | lon, & ! Longitude (radians) |
---|
575 | lat, & ! Latitude (radians) |
---|
576 | z ! Altitude (m) |
---|
577 | |
---|
578 | ! Great circle radius from the perturbation centerpoint |
---|
579 | REAL(8) :: gr |
---|
580 | |
---|
581 | ! Approximately spherical radius from the perturbation centerpoint |
---|
582 | REAL(8) :: Rtheta |
---|
583 | |
---|
584 | gr = a*acos(sin(pert_latc*deg2rad)*sin(lat) + & |
---|
585 | (cos(pert_latc*deg2rad)*cos(lat)*cos(lon-pert_lonc*deg2rad))) |
---|
586 | |
---|
587 | Rtheta = sqrt((gr/pert_rh)**2 + ((z - pert_zc) / pert_rz)**2) |
---|
588 | |
---|
589 | if (Rtheta .le. 1.d0) then |
---|
590 | thermal_perturbation = pert_dtheta * (cos(0.5d0 * pi * Rtheta))**2 |
---|
591 | else |
---|
592 | thermal_perturbation = 0.0d0 |
---|
593 | end if |
---|
594 | |
---|
595 | END FUNCTION thermal_perturbation |
---|
596 | |
---|
597 | !----------------------------------------------------------------------- |
---|
598 | ! Calculate the reference zonal velocity |
---|
599 | !----------------------------------------------------------------------- |
---|
600 | REAL(8) FUNCTION zonal_velocity(z, lat) |
---|
601 | |
---|
602 | IMPLICIT NONE |
---|
603 | |
---|
604 | REAL(8), INTENT(IN) :: z, lat |
---|
605 | |
---|
606 | if (z .le. zs - zt) then |
---|
607 | zonal_velocity = us * (z / zs) - uc |
---|
608 | elseif (abs(z - zs) .le. zt) then |
---|
609 | zonal_velocity = & |
---|
610 | (-4.0d0/5.0d0 + 3.0d0*z/zs - 5.0d0/4.0d0*(z**2)/(zs**2)) * us - uc |
---|
611 | else |
---|
612 | zonal_velocity = us - uc |
---|
613 | end if |
---|
614 | |
---|
615 | zonal_velocity = zonal_velocity * cos(lat) |
---|
616 | |
---|
617 | END FUNCTION zonal_velocity |
---|
618 | |
---|
619 | !----------------------------------------------------------------------- |
---|
620 | ! Calculate pointwise theta at the equator at the given altitude |
---|
621 | !----------------------------------------------------------------------- |
---|
622 | REAL(8) FUNCTION equator_theta(z) |
---|
623 | |
---|
624 | IMPLICIT NONE |
---|
625 | |
---|
626 | REAL(8), INTENT(IN) :: z |
---|
627 | |
---|
628 | if (z .le. z_tr) then |
---|
629 | equator_theta = & |
---|
630 | theta0 + (theta_tr - theta0) * (z / z_tr)**(1.25d0) |
---|
631 | else |
---|
632 | equator_theta = & |
---|
633 | theta_tr * exp(g/cp/T_tr * (z - z_tr)) |
---|
634 | end if |
---|
635 | |
---|
636 | END FUNCTION equator_theta |
---|
637 | |
---|
638 | !----------------------------------------------------------------------- |
---|
639 | ! Calculate pointwise relative humidity (in %) at the equator at the |
---|
640 | ! given altitude |
---|
641 | !----------------------------------------------------------------------- |
---|
642 | REAL(8) FUNCTION equator_relative_humidity(z) |
---|
643 | |
---|
644 | IMPLICIT NONE |
---|
645 | |
---|
646 | REAL(8), INTENT(IN) :: z |
---|
647 | |
---|
648 | if (z .le. z_tr) then |
---|
649 | equator_relative_humidity = 1.0d0 - 0.75d0 * (z / z_tr)**(1.25d0) |
---|
650 | else |
---|
651 | equator_relative_humidity = 0.25d0 |
---|
652 | end if |
---|
653 | |
---|
654 | END FUNCTION equator_relative_humidity |
---|
655 | |
---|
656 | !----------------------------------------------------------------------- |
---|
657 | ! Calculate saturation mixing ratio (in kg/kg) in terms of pressure |
---|
658 | ! (in Pa) and temperature (in K) |
---|
659 | !----------------------------------------------------------------------- |
---|
660 | REAL(8) FUNCTION saturation_mixing_ratio(p, T) |
---|
661 | |
---|
662 | IMPLICIT NONE |
---|
663 | |
---|
664 | REAL(8), INTENT(IN) :: & |
---|
665 | p, & ! Pressure in Pa |
---|
666 | T ! Temperature |
---|
667 | |
---|
668 | saturation_mixing_ratio = & |
---|
669 | 380.d0 / p * exp(17.27d0 * (T - 273.d0) / (T - 36.d0)) |
---|
670 | |
---|
671 | if (saturation_mixing_ratio > 0.014) then |
---|
672 | saturation_mixing_ratio = 0.014 |
---|
673 | end if |
---|
674 | |
---|
675 | END FUNCTION saturation_mixing_ratio |
---|
676 | |
---|
677 | !----------------------------------------------------------------------- |
---|
678 | ! Calculate coefficients for a Lagrangian polynomial |
---|
679 | !----------------------------------------------------------------------- |
---|
680 | SUBROUTINE lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
681 | |
---|
682 | IMPLICIT NONE |
---|
683 | |
---|
684 | ! Number of points to fit |
---|
685 | INTEGER(4), INTENT(IN) :: npts |
---|
686 | |
---|
687 | ! Sample points to fit |
---|
688 | REAL(8), DIMENSION(npts), INTENT(IN) :: x |
---|
689 | |
---|
690 | ! Computed coefficients |
---|
691 | REAL(8), DIMENSION(npts), INTENT(OUT) :: coeffs |
---|
692 | |
---|
693 | ! Point at which sample is taken |
---|
694 | REAL(8), INTENT(IN) :: xs |
---|
695 | |
---|
696 | ! Loop indices |
---|
697 | INTEGER(4) :: i, j |
---|
698 | |
---|
699 | ! Compute the Lagrangian polynomial coefficients |
---|
700 | do i = 1, npts |
---|
701 | coeffs(i) = 1.0d0 |
---|
702 | do j = 1, npts |
---|
703 | if (i .eq. j) then |
---|
704 | cycle |
---|
705 | end if |
---|
706 | coeffs(i) = coeffs(i) * (xs - x(j)) / (x(i) - x(j)) |
---|
707 | end do |
---|
708 | end do |
---|
709 | |
---|
710 | END SUBROUTINE lagrangian_polynomial_coeffs |
---|
711 | |
---|
712 | !----------------------------------------------------------------------- |
---|
713 | ! Calculate coefficients of the derivative of a Lagrangian polynomial |
---|
714 | !----------------------------------------------------------------------- |
---|
715 | SUBROUTINE diff_lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
716 | |
---|
717 | IMPLICIT NONE |
---|
718 | |
---|
719 | ! Number of points to fit |
---|
720 | INTEGER(4), INTENT(IN) :: npts |
---|
721 | |
---|
722 | ! Sample points to fit |
---|
723 | REAL(8), DIMENSION(npts), INTENT(IN) :: x |
---|
724 | |
---|
725 | ! Computed coefficients |
---|
726 | REAL(8), DIMENSION(npts), INTENT(OUT) :: coeffs |
---|
727 | |
---|
728 | ! Point at which sample is taken |
---|
729 | REAL(8), INTENT(IN) :: xs |
---|
730 | |
---|
731 | ! Loop indices |
---|
732 | INTEGER(4) :: i, j, imatch |
---|
733 | |
---|
734 | ! Buffer sum |
---|
735 | REAL(8) :: coeffsum, differential |
---|
736 | |
---|
737 | ! Check if xs is equivalent to one of the values of x |
---|
738 | imatch = (-1) |
---|
739 | do i = 1, npts |
---|
740 | if (abs(xs - x(i)) < 1.0d-14) then |
---|
741 | imatch = i |
---|
742 | exit |
---|
743 | end if |
---|
744 | end do |
---|
745 | |
---|
746 | ! Equivalence detected; special treatment required |
---|
747 | if (imatch .ne. (-1)) then |
---|
748 | do i = 1, npts |
---|
749 | coeffs(i) = 1.0d0 |
---|
750 | coeffsum = 0.0d0 |
---|
751 | |
---|
752 | do j = 1, npts |
---|
753 | if ((j .eq. i) .or. (j .eq. imatch)) then |
---|
754 | cycle |
---|
755 | end if |
---|
756 | |
---|
757 | coeffs(i) = coeffs(i) * (xs - x(j)) / (x(i) - x(j)) |
---|
758 | coeffsum = coeffsum + 1.0 / (xs - x(j)) |
---|
759 | end do |
---|
760 | |
---|
761 | if (i .ne. imatch) then |
---|
762 | coeffs(i) = coeffs(i) & |
---|
763 | * (1.0 + (xs - x(imatch)) * coeffsum) & |
---|
764 | / (x(i) - x(imatch)) |
---|
765 | else |
---|
766 | coeffs(i) = coeffs(i) * coeffsum |
---|
767 | end if |
---|
768 | end do |
---|
769 | |
---|
770 | ! No equivalence; simply differentiate Lagrangian fit |
---|
771 | else |
---|
772 | call lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
773 | |
---|
774 | do i = 1, npts |
---|
775 | differential = 0.0d0 |
---|
776 | do j = 1, npts |
---|
777 | if (i .eq. j) then |
---|
778 | cycle |
---|
779 | end if |
---|
780 | differential = differential + 1.0 / (xs - x(j)) |
---|
781 | end do |
---|
782 | coeffs(i) = coeffs(i) * differential |
---|
783 | end do |
---|
784 | end if |
---|
785 | |
---|
786 | END SUBROUTINE |
---|
787 | |
---|
788 | END MODULE dcmip2016_supercell_mod |
---|