1 | MODULE dcmip_initial_conditions_test_1_2_3 |
---|
2 | |
---|
3 | !======================================================================= |
---|
4 | ! |
---|
5 | ! Functions for setting up initial conditions for the dynamical core tests: |
---|
6 | ! |
---|
7 | ! 11 - Deformational Advection Test |
---|
8 | ! 12 - Hadley Cell Advection Test |
---|
9 | ! 13 - Orography Advection Test |
---|
10 | ! 20 - Impact of orography on a steady-state at rest |
---|
11 | ! 21 and 22 - Non-Hydrostatic Mountain Waves Over A Schaer-Type Mountain without and with vertical wind shear |
---|
12 | ! 31 - Non-Hydrostatic Gravity Waves |
---|
13 | ! |
---|
14 | ! Given a point specified by: |
---|
15 | ! lon longitude (radians) |
---|
16 | ! lat latitude (radians) |
---|
17 | ! p/z pressure/height |
---|
18 | ! the functions will return: |
---|
19 | ! u zonal wind (m s^-1) |
---|
20 | ! v meridional wind (m s^-1) |
---|
21 | ! w vertical velocity (m s^-1) |
---|
22 | ! t temperature (K) |
---|
23 | ! phis surface geopotential (m^2 s^-2) |
---|
24 | ! ps surface pressure (Pa) |
---|
25 | ! rho density (kj m^-3) |
---|
26 | ! q specific humidity (kg/kg) |
---|
27 | ! qi tracers (kg/kg) |
---|
28 | ! p pressure if height based (Pa) |
---|
29 | ! |
---|
30 | ! |
---|
31 | ! Authors: James Kent, Paul Ullrich, Christiane Jablonowski |
---|
32 | ! (University of Michigan, dcmip@ucar.edu) |
---|
33 | ! version 4 |
---|
34 | ! July/8/2012 |
---|
35 | ! |
---|
36 | ! Change log: (v3, June/8/2012, v4 July/8/2012, v5 July/20/2012) |
---|
37 | ! |
---|
38 | ! v2: bug fixes in the tracer initialization for height-based models |
---|
39 | ! v3: test 3-1: the density is now initialized with the unperturbed background temperature (not the perturbed temperature) |
---|
40 | ! v3: constants converted to double precision |
---|
41 | ! v4: modified tracers in test 1-1, now with cutoff altitudes. Outside of the vertical domain all tracers are set to 0 |
---|
42 | ! v4: modified cos-term in vertical velocity (now cos(2 pi t/tau)) in test 1-1, now completing one full up and down cycle |
---|
43 | ! v4: added subroutine test1_advection_orography for test 1-3 |
---|
44 | ! v4: added subroutine test2_steady_state_mountain for test 2-0 |
---|
45 | ! v4: modified parameter list for tests 2-1 and 2-2 (routine test2_schaer_mountain): addition of parameters hybrid_eta, hyam, hybm |
---|
46 | ! if the logical flag hybrid_eta is true then the pressure in pressure-based model with hybrid sigma-p (eta) coordinates is |
---|
47 | ! computed internally. In that case the hybrid coefficients hyam and hybm need to be supplied via the parameter list, |
---|
48 | ! otherwise they are not used. |
---|
49 | ! v5: Change in test 11 - change to u and w, cutoff altitudes (introduced in v4) are removed again |
---|
50 | ! v5: Change in test 12 - velocities multiplies by rho0/rho, different w0 and vertical location of the initial tracer |
---|
51 | ! |
---|
52 | ! |
---|
53 | !======================================================================= |
---|
54 | |
---|
55 | USE prec |
---|
56 | |
---|
57 | IMPLICIT NONE |
---|
58 | |
---|
59 | PRIVATE |
---|
60 | PUBLIC :: test3_gravity_wave, & ! (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) |
---|
61 | test2_steady_state_mountain, & |
---|
62 | test2_schaer_mountain |
---|
63 | !----------------------------------------------------------------------- |
---|
64 | ! Physical Parameters |
---|
65 | !----------------------------------------------------------------------- |
---|
66 | |
---|
67 | REAL(rstd), parameter :: & |
---|
68 | a = 6371220.0d0, & ! Earth's Radius (m) |
---|
69 | Rd = 287.0d0, & ! Ideal gas const dry air (J kg^-1 K^1) |
---|
70 | g = 9.80616d0, & ! Gravity (m s^2) |
---|
71 | cp = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1) |
---|
72 | pi = 4.d0*atan(1.d0) ! pi |
---|
73 | |
---|
74 | !----------------------------------------------------------------------- |
---|
75 | ! Additional constants |
---|
76 | !----------------------------------------------------------------------- |
---|
77 | |
---|
78 | real(rstd), parameter :: p0 = 100000.d0 ! reference pressure (Pa) |
---|
79 | |
---|
80 | |
---|
81 | CONTAINS |
---|
82 | |
---|
83 | !========================================================================================== |
---|
84 | ! TEST CASE 11 - PURE ADVECTION - 3D DEFORMATIONAL FLOW |
---|
85 | !========================================================================================== |
---|
86 | |
---|
87 | ! The 3D deformational flow test is based on the deformational flow test of Nair and Lauritzen (JCP 2010), |
---|
88 | ! with a prescribed vertical wind velocity which makes the test truly 3D. An unscaled planet (with scale parameter |
---|
89 | ! X = 1) is selected. |
---|
90 | |
---|
91 | SUBROUTINE test1_advection_deformation (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4) |
---|
92 | |
---|
93 | IMPLICIT NONE |
---|
94 | !----------------------------------------------------------------------- |
---|
95 | ! input/output params parameters at given location |
---|
96 | !----------------------------------------------------------------------- |
---|
97 | |
---|
98 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
99 | lat, & ! Latitude (radians) |
---|
100 | z ! Height (m) |
---|
101 | |
---|
102 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
103 | |
---|
104 | integer, intent(in) :: zcoords ! 0 or 1 see below |
---|
105 | |
---|
106 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
107 | v, & ! Meridional wind (m s^-1) |
---|
108 | w, & ! Vertical Velocity (m s^-1) |
---|
109 | t, & ! Temperature (K) |
---|
110 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
111 | ps, & ! Surface Pressure (Pa) |
---|
112 | rho, & ! density (kg m^-3) |
---|
113 | q, & ! Specific Humidity (kg/kg) |
---|
114 | q1, & ! Tracer q1 (kg/kg) |
---|
115 | q2, & ! Tracer q2 (kg/kg) |
---|
116 | q3, & ! Tracer q3 (kg/kg) |
---|
117 | q4 ! Tracer q4 (kg/kg) |
---|
118 | |
---|
119 | ! if zcoords = 1, then we use z and output p |
---|
120 | ! if zcoords = 0, then we use p |
---|
121 | |
---|
122 | !----------------------------------------------------------------------- |
---|
123 | ! test case parameters |
---|
124 | !----------------------------------------------------------------------- |
---|
125 | real(rstd), parameter :: tau = 12.d0 * 86400.d0, & ! period of motion 12 days |
---|
126 | u0 = (2.d0*pi*a)/tau, & ! 2 pi a / 12 days |
---|
127 | k0 = (10.d0*a)/tau, & ! Velocity Magnitude |
---|
128 | omega0 = (23000.d0*pi)/tau, & ! Velocity Magnitude |
---|
129 | T0 = 300.d0, & ! temperature |
---|
130 | H = Rd * T0 / g, & ! scale height |
---|
131 | RR = 1.d0/2.d0, & ! horizontal half width divided by 'a' |
---|
132 | ZZ = 1000.d0, & ! vertical half width |
---|
133 | z0 = 5000.d0, & ! center point in z |
---|
134 | lambda0 = 5.d0*pi/6.d0, & ! center point in longitudes |
---|
135 | lambda1 = 7.d0*pi/6.d0, & ! center point in longitudes |
---|
136 | phi0 = 0.d0, & ! center point in latitudes |
---|
137 | phi1 = 0.d0 |
---|
138 | |
---|
139 | real(rstd) :: height ! The height of the model levels |
---|
140 | real(rstd) :: ptop ! Model top in p |
---|
141 | real(rstd) :: sin_tmp, cos_tmp, sin_tmp2, cos_tmp2 ! Calculate great circle distances |
---|
142 | real(rstd) :: d1, d2, r, r2, d3, d4 ! For tracer calculations |
---|
143 | real(rstd) :: s, bs ! Shape function, and parameter |
---|
144 | real(rstd) :: lonp ! Translational longitude, depends on time |
---|
145 | real(rstd) :: ud ! Divergent part of u |
---|
146 | real(rstd) :: time ! Initially set to zero seconds, needs |
---|
147 | ! to be modified when used in dycore |
---|
148 | |
---|
149 | !----------------------------------------------------------------------- |
---|
150 | ! HEIGHT AND PRESSURE |
---|
151 | !----------------------------------------------------------------------- |
---|
152 | |
---|
153 | ! Height and pressure are aligned (p = p0 exp(-z/H)) |
---|
154 | |
---|
155 | if (zcoords .eq. 1) then |
---|
156 | |
---|
157 | height = z |
---|
158 | p = p0 * exp(-z/H) |
---|
159 | |
---|
160 | else |
---|
161 | |
---|
162 | height = H * log(p0/p) |
---|
163 | |
---|
164 | endif |
---|
165 | |
---|
166 | ! Model top in p |
---|
167 | |
---|
168 | ptop = p0*exp(-12000.d0/H) |
---|
169 | |
---|
170 | !----------------------------------------------------------------------- |
---|
171 | ! THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED |
---|
172 | ! IN THE DYNAMICAL CORE |
---|
173 | !----------------------------------------------------------------------- |
---|
174 | |
---|
175 | ! These are initial conditions hence time = 0 |
---|
176 | |
---|
177 | time = 0.d0 |
---|
178 | |
---|
179 | ! Translational longitude = longitude when time = 0 |
---|
180 | |
---|
181 | lonp = lon - 2.d0*pi*time/tau |
---|
182 | |
---|
183 | ! Shape function |
---|
184 | !******** |
---|
185 | ! change in version 5: shape function |
---|
186 | !******** |
---|
187 | bs = 0.2 |
---|
188 | s = 1.0 + exp( (ptop-p0)/(bs*ptop) ) - exp( (p-p0)/(bs*ptop)) - exp( (ptop-p)/(bs*ptop)) |
---|
189 | |
---|
190 | ! Zonal Velocity |
---|
191 | !******** |
---|
192 | ! change in version 5: ud |
---|
193 | !******** |
---|
194 | |
---|
195 | ud = (omega0*a)/(bs*ptop) * cos(lonp) * (cos(lat)**2.0) * cos(2.0*pi*time/tau) * & |
---|
196 | ( - exp( (p-p0)/(bs*ptop)) + exp( (ptop-p)/(bs*ptop)) ) |
---|
197 | |
---|
198 | u = k0*sin(lonp)*sin(lonp)*sin(2.d0*lat)*cos(pi*time/tau) + u0*cos(lat) + ud |
---|
199 | |
---|
200 | ! Meridional Velocity |
---|
201 | |
---|
202 | v = k0*sin(2.d0*lonp)*cos(lat)*cos(pi*time/tau) |
---|
203 | |
---|
204 | ! Vertical Velocity - can be changed to vertical pressure velocity by |
---|
205 | ! omega = -(g*p)/(Rd*T0)*w |
---|
206 | ! |
---|
207 | !******** |
---|
208 | ! change in version 4: cos(2.0*pi*time/tau) is now used instead of cos(pi*time/tau) |
---|
209 | !******** |
---|
210 | !******** |
---|
211 | ! change in version 5: shape function in w |
---|
212 | !******** |
---|
213 | w = -((Rd*T0)/(g*p))*omega0*sin(lonp)*cos(lat)*cos(2.0*pi*time/tau)*s |
---|
214 | |
---|
215 | !----------------------------------------------------------------------- |
---|
216 | ! TEMPERATURE IS CONSTANT 300 K |
---|
217 | !----------------------------------------------------------------------- |
---|
218 | |
---|
219 | t = T0 |
---|
220 | |
---|
221 | !----------------------------------------------------------------------- |
---|
222 | ! PHIS (surface geopotential) |
---|
223 | !----------------------------------------------------------------------- |
---|
224 | |
---|
225 | phis = 0.d0 |
---|
226 | |
---|
227 | !----------------------------------------------------------------------- |
---|
228 | ! PS (surface pressure) |
---|
229 | !----------------------------------------------------------------------- |
---|
230 | |
---|
231 | ps = p0 |
---|
232 | |
---|
233 | !----------------------------------------------------------------------- |
---|
234 | ! RHO (density) |
---|
235 | !----------------------------------------------------------------------- |
---|
236 | |
---|
237 | rho = p/(Rd*t) |
---|
238 | |
---|
239 | !----------------------------------------------------------------------- |
---|
240 | ! initialize Q, set to zero |
---|
241 | !----------------------------------------------------------------------- |
---|
242 | |
---|
243 | q = 0.d0 |
---|
244 | |
---|
245 | !----------------------------------------------------------------------- |
---|
246 | ! initialize tracers |
---|
247 | !----------------------------------------------------------------------- |
---|
248 | |
---|
249 | ! Tracer 1 - Cosine Bells |
---|
250 | |
---|
251 | ! To calculate great circle distance |
---|
252 | |
---|
253 | sin_tmp = sin(lat) * sin(phi0) |
---|
254 | cos_tmp = cos(lat) * cos(phi0) |
---|
255 | sin_tmp2 = sin(lat) * sin(phi1) |
---|
256 | cos_tmp2 = cos(lat) * cos(phi1) |
---|
257 | |
---|
258 | ! great circle distance without 'a' |
---|
259 | |
---|
260 | r = ACOS (sin_tmp + cos_tmp*cos(lon-lambda0)) |
---|
261 | r2 = ACOS (sin_tmp2 + cos_tmp2*cos(lon-lambda1)) |
---|
262 | d1 = min( 1., (r/RR)**2 + ((height-z0)/ZZ)**2 ) |
---|
263 | d2 = min( 1., (r2/RR)**2 + ((height-z0)/ZZ)**2 ) |
---|
264 | |
---|
265 | q1 = 0.5d0 * (1.d0 + cos(pi*d1)) + 0.5d0 * (1.d0 + cos(pi*d2)) |
---|
266 | |
---|
267 | ! Tracer 2 - Correlated Cosine Bells |
---|
268 | |
---|
269 | q2 = 0.9d0 - 0.8d0*q1**2 |
---|
270 | |
---|
271 | ! Tracer 3 - Slotted Ellipse |
---|
272 | |
---|
273 | ! Make the ellipse |
---|
274 | |
---|
275 | if (d1 .le. RR) then |
---|
276 | q3 = 1.d0 |
---|
277 | elseif (d2 .le. RR) then |
---|
278 | q3 = 1.d0 |
---|
279 | else |
---|
280 | q3 = 0.1d0 |
---|
281 | endif |
---|
282 | |
---|
283 | ! Put in the slot |
---|
284 | |
---|
285 | if (height .gt. z0 .and. abs(lat) .lt. 0.125d0) then |
---|
286 | |
---|
287 | q3 = 0.1d0 |
---|
288 | |
---|
289 | endif |
---|
290 | |
---|
291 | ! Tracer 4: q4 is chosen so that, in combination with the other three tracer |
---|
292 | ! fields with weight (3/10), the sum is equal to one |
---|
293 | |
---|
294 | q4 = 1.d0 - 0.3d0*(q1+q2+q3) |
---|
295 | |
---|
296 | !************ |
---|
297 | ! change in version 4: added cutoff altitudes, tracers are set to zero outside this region |
---|
298 | ! prevents tracers from being trapped near the bottom and top of the domain |
---|
299 | !************ |
---|
300 | ! use a cutoff altitude for |
---|
301 | ! tracer 2 3 and 4 |
---|
302 | ! Set them to zero outside `buffer zone' |
---|
303 | !************ |
---|
304 | ! change in version 5: change from v4 reversed, no cutoff altitudes due to continuity equation being satisfied |
---|
305 | ! commented out below |
---|
306 | !************ |
---|
307 | |
---|
308 | !if (height .gt. (z0+1.25*ZZ) .or. height .lt. (z0-1.25*ZZ)) then |
---|
309 | |
---|
310 | ! q2 = 0.0 |
---|
311 | ! q3 = 0.0 |
---|
312 | ! q4 = 0.0 |
---|
313 | |
---|
314 | !endif |
---|
315 | |
---|
316 | |
---|
317 | |
---|
318 | |
---|
319 | END SUBROUTINE test1_advection_deformation |
---|
320 | |
---|
321 | |
---|
322 | |
---|
323 | |
---|
324 | |
---|
325 | !========================================================================================== |
---|
326 | ! TEST CASE 12 - PURE ADVECTION - 3D HADLEY-LIKE FLOW |
---|
327 | !========================================================================================== |
---|
328 | |
---|
329 | SUBROUTINE test1_advection_hadley (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1) |
---|
330 | |
---|
331 | IMPLICIT NONE |
---|
332 | !----------------------------------------------------------------------- |
---|
333 | ! input/output params parameters at given location |
---|
334 | !----------------------------------------------------------------------- |
---|
335 | |
---|
336 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
337 | lat, & ! Latitude (radians) |
---|
338 | z ! Height (m) |
---|
339 | |
---|
340 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
341 | |
---|
342 | integer, intent(in) :: zcoords ! 0 or 1 see below |
---|
343 | |
---|
344 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
345 | v, & ! Meridional wind (m s^-1) |
---|
346 | w, & ! Vertical Velocity (m s^-1) |
---|
347 | t, & ! Temperature (K) |
---|
348 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
349 | ps, & ! Surface Pressure (Pa) |
---|
350 | rho, & ! density (kg m^-3) |
---|
351 | q, & ! Specific Humidity (kg/kg) |
---|
352 | q1 ! Tracer q1 (kg/kg) |
---|
353 | |
---|
354 | ! if zcoords = 1, then we use z and output p |
---|
355 | ! if zcoords = 0, then we use p |
---|
356 | |
---|
357 | !----------------------------------------------------------------------- |
---|
358 | ! test case parameters |
---|
359 | !----------------------------------------------------------------------- |
---|
360 | real(rstd), parameter :: tau = 1.d0 * 86400.d0, & ! period of motion 1 day (in s) |
---|
361 | u0 = 40.d0, & ! Zonal velocity magnitude (m/s) |
---|
362 | w0 = 0.15d0, & ! Vertical velocity magnitude (m/s), changed in v5 |
---|
363 | T0 = 300.d0, & ! temperature (K) |
---|
364 | H = Rd * T0 / g, & ! scale height |
---|
365 | K = 5.d0, & ! number of Hadley-like cells |
---|
366 | z1 = 2000.d0, & ! position of lower tracer bound (m), changed in v5 |
---|
367 | z2 = 5000.d0, & ! position of upper tracer bound (m), changed in v5 |
---|
368 | z0 = 0.5d0*(z1+z2), & ! midpoint (m) |
---|
369 | ztop = 12000.d0 ! model top (m) |
---|
370 | |
---|
371 | real(rstd) :: rho0 ! reference density at z=0 m |
---|
372 | real(rstd) :: height ! Model level heights |
---|
373 | real(rstd) :: time ! Initially set to zero seconds, needs |
---|
374 | ! to be modified when used in dycore |
---|
375 | |
---|
376 | !----------------------------------------------------------------------- |
---|
377 | ! HEIGHT AND PRESSURE |
---|
378 | !----------------------------------------------------------------------- |
---|
379 | |
---|
380 | ! Height and pressure are aligned (p = p0 exp(-z/H)) |
---|
381 | |
---|
382 | if (zcoords .eq. 1) then |
---|
383 | |
---|
384 | height = z |
---|
385 | p = p0 * exp(-z/H) |
---|
386 | |
---|
387 | else |
---|
388 | |
---|
389 | height = H * log(p0/p) |
---|
390 | |
---|
391 | endif |
---|
392 | |
---|
393 | !----------------------------------------------------------------------- |
---|
394 | ! TEMPERATURE IS CONSTANT 300 K |
---|
395 | !----------------------------------------------------------------------- |
---|
396 | |
---|
397 | t = T0 |
---|
398 | |
---|
399 | !----------------------------------------------------------------------- |
---|
400 | ! PHIS (surface geopotential) |
---|
401 | !----------------------------------------------------------------------- |
---|
402 | |
---|
403 | phis = 0.d0 |
---|
404 | |
---|
405 | !----------------------------------------------------------------------- |
---|
406 | ! PS (surface pressure) |
---|
407 | !----------------------------------------------------------------------- |
---|
408 | |
---|
409 | ps = p0 |
---|
410 | |
---|
411 | !----------------------------------------------------------------------- |
---|
412 | ! RHO (density) |
---|
413 | !----------------------------------------------------------------------- |
---|
414 | |
---|
415 | rho = p/(Rd*t) |
---|
416 | rho0 = p0/(Rd*t) |
---|
417 | |
---|
418 | !----------------------------------------------------------------------- |
---|
419 | ! THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED |
---|
420 | ! IN THE DYNAMICAL CORE |
---|
421 | !----------------------------------------------------------------------- |
---|
422 | |
---|
423 | time = 0.d0 |
---|
424 | |
---|
425 | ! Zonal Velocity |
---|
426 | |
---|
427 | u = u0*cos(lat) |
---|
428 | |
---|
429 | ! Meridional Velocity |
---|
430 | |
---|
431 | !************ |
---|
432 | ! change in version 5: multiply v and w by rho0/rho |
---|
433 | !************ |
---|
434 | |
---|
435 | v = -(rho0/rho) * (a*w0*pi)/(K*ztop) *cos(lat)*sin(K*lat)*cos(pi*height/ztop)*cos(pi*time/tau) |
---|
436 | |
---|
437 | ! Vertical Velocity - can be changed to vertical pressure velocity by |
---|
438 | ! omega = -g*rho*w |
---|
439 | |
---|
440 | w = (rho0/rho) *(w0/K)*(-2.d0*sin(K*lat)*sin(lat) + K*cos(lat)*cos(K*lat)) & |
---|
441 | *sin(pi*height/ztop)*cos(pi*time/tau) |
---|
442 | |
---|
443 | |
---|
444 | !----------------------------------------------------------------------- |
---|
445 | ! initialize Q, set to zero |
---|
446 | !----------------------------------------------------------------------- |
---|
447 | |
---|
448 | q = 0.d0 |
---|
449 | |
---|
450 | !----------------------------------------------------------------------- |
---|
451 | ! initialize tracers |
---|
452 | !----------------------------------------------------------------------- |
---|
453 | |
---|
454 | ! Tracer 1 - Layer |
---|
455 | |
---|
456 | if (height .lt. z2 .and. height .gt. z1) then |
---|
457 | |
---|
458 | q1 = 0.5d0 * (1.d0 + cos( 2.d0*pi*(height-z0)/(z2-z1) ) ) |
---|
459 | |
---|
460 | else |
---|
461 | |
---|
462 | q1 = 0.d0 |
---|
463 | |
---|
464 | endif |
---|
465 | |
---|
466 | END SUBROUTINE test1_advection_hadley |
---|
467 | |
---|
468 | |
---|
469 | |
---|
470 | !========================================================================================== |
---|
471 | ! TEST CASE 13 - HORIZONTAL ADVECTION OF THIN CLOUD-LIKE TRACERS IN THE PRESENCE OF OROGRAPHY |
---|
472 | !========================================================================================== |
---|
473 | |
---|
474 | SUBROUTINE test1_advection_orography (lon,lat,p,z,zcoords,cfv,hybrid_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4) |
---|
475 | |
---|
476 | IMPLICIT NONE |
---|
477 | !----------------------------------------------------------------------- |
---|
478 | ! input/output params parameters at given location |
---|
479 | !----------------------------------------------------------------------- |
---|
480 | |
---|
481 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
482 | lat, & ! Latitude (radians) |
---|
483 | z, & ! Height (m) |
---|
484 | hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint |
---|
485 | hybm, & ! B coefficient for hybrid-eta coordinate, at model level midpoint |
---|
486 | gc ! bar{z} for Gal-Chen coordinate |
---|
487 | |
---|
488 | logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used |
---|
489 | ! if set to .true., then the pressure will be computed via the |
---|
490 | ! hybrid coefficients hyam and hybm, they need to be initialized |
---|
491 | ! if set to .false. (for pressure-based models): the pressure is already pre-computed |
---|
492 | ! and is an input value for this routine |
---|
493 | ! for height-based models: pressure will always be computed based on the height and |
---|
494 | ! hybrid_eta is not used |
---|
495 | |
---|
496 | ! Note that we only use hyam and hybm for the hybrid-eta coordiantes, and we only use |
---|
497 | ! gc for the Gal-Chen coordinates. If not required then they become dummy variables |
---|
498 | |
---|
499 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
500 | |
---|
501 | integer, intent(in) :: zcoords ! 0 or 1 see below |
---|
502 | integer, intent(in) :: cfv ! 0, 1 or 2 see below |
---|
503 | |
---|
504 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
505 | v, & ! Meridional wind (m s^-1) |
---|
506 | w, & ! Vertical Velocity (m s^-1) |
---|
507 | t, & ! Temperature (K) |
---|
508 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
509 | ps, & ! Surface Pressure (Pa) |
---|
510 | rho, & ! density (kg m^-3) |
---|
511 | q, & ! Specific Humidity (kg/kg) |
---|
512 | q1, & ! Tracer q1 (kg/kg) |
---|
513 | q2, & ! Tracer q2 (kg/kg) |
---|
514 | q3, & ! Tracer q3 (kg/kg) |
---|
515 | q4 ! Tracer q4 (kg/kg) |
---|
516 | |
---|
517 | ! if zcoords = 1, then we use z and output p |
---|
518 | ! if zcoords = 0, then we use p |
---|
519 | |
---|
520 | ! if cfv = 0 we assume that our horizontal velocities are not coordinate following |
---|
521 | ! if cfv = 1 then our velocities follow hybrid eta coordinates and we need to specify w |
---|
522 | ! if cfv = 2 then our velocities follow Gal-Chen coordinates and we need to specify w |
---|
523 | |
---|
524 | ! In hybrid-eta coords: p = hyam p0 + hybm ps |
---|
525 | ! In Gal-Chen coords: z = zs + (gc/ztop)*(ztop - zs) |
---|
526 | |
---|
527 | ! if other orography-following coordinates are used, the w wind needs to be newly derived for them |
---|
528 | |
---|
529 | !----------------------------------------------------------------------- |
---|
530 | ! test case parameters |
---|
531 | !----------------------------------------------------------------------- |
---|
532 | real(rstd), parameter :: tau = 12.d0 * 86400.d0, & ! period of motion 12 days (s) |
---|
533 | u0 = 2.d0*pi*a/tau, & ! Velocity Magnitude (m/s) |
---|
534 | T0 = 300.d0, & ! temperature (K) |
---|
535 | H = Rd * T0 / g, & ! scale height (m) |
---|
536 | alpha = pi/6.d0, & ! rotation angle (radians), 30 degrees |
---|
537 | lambdam = 3.d0*pi/2.d0, & ! mountain longitude center point (radians) |
---|
538 | phim = 0.d0, & ! mountain latitude center point (radians) |
---|
539 | h0 = 2000.d0, & ! peak height of the mountain range (m) |
---|
540 | Rm = 3.d0*pi/4.d0, & ! mountain radius (radians) |
---|
541 | zetam = pi/16.d0, & ! mountain oscillation half-width (radians) |
---|
542 | lambdap = pi/2.d0, & ! cloud-like tracer longitude center point (radians) |
---|
543 | phip = 0.d0, & ! cloud-like tracer latitude center point (radians) |
---|
544 | Rp = pi/4.d0, & ! cloud-like tracer radius (radians) |
---|
545 | zp1 = 3050.d0, & ! midpoint of first (lowermost) tracer (m) |
---|
546 | zp2 = 5050.d0, & ! midpoint of second tracer (m) |
---|
547 | zp3 = 8200.d0, & ! midpoint of third (topmost) tracer (m) |
---|
548 | dzp1 = 1000.d0, & ! thickness of first (lowermost) tracer (m) |
---|
549 | dzp2 = 1000.d0, & ! thickness of second tracer (m) |
---|
550 | dzp3 = 400.d0, & ! thickness of third (topmost) tracer (m) |
---|
551 | ztop = 12000.d0 ! model top (m) |
---|
552 | |
---|
553 | real(rstd) :: height ! Model level heights (m) |
---|
554 | real(rstd) :: r ! Great circle distance (radians) |
---|
555 | real(rstd) :: rz ! height differences |
---|
556 | real(rstd) :: zs ! Surface elevation (m) |
---|
557 | |
---|
558 | |
---|
559 | |
---|
560 | !----------------------------------------------------------------------- |
---|
561 | ! PHIS (surface geopotential) |
---|
562 | !----------------------------------------------------------------------- |
---|
563 | |
---|
564 | r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) |
---|
565 | |
---|
566 | if (r .lt. Rm) then |
---|
567 | |
---|
568 | zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0 |
---|
569 | |
---|
570 | else |
---|
571 | |
---|
572 | zs = 0.d0 |
---|
573 | |
---|
574 | endif |
---|
575 | |
---|
576 | phis = g*zs |
---|
577 | |
---|
578 | !----------------------------------------------------------------------- |
---|
579 | ! PS (surface pressure) |
---|
580 | !----------------------------------------------------------------------- |
---|
581 | |
---|
582 | ps = p0 * exp(-zs/H) |
---|
583 | |
---|
584 | |
---|
585 | !----------------------------------------------------------------------- |
---|
586 | ! HEIGHT AND PRESSURE |
---|
587 | !----------------------------------------------------------------------- |
---|
588 | |
---|
589 | ! Height and pressure are aligned (p = p0 exp(-z/H)) |
---|
590 | |
---|
591 | if (zcoords .eq. 1) then |
---|
592 | |
---|
593 | height = z |
---|
594 | p = p0 * exp(-z/H) |
---|
595 | |
---|
596 | else |
---|
597 | |
---|
598 | if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients |
---|
599 | height = H * log(p0/p) |
---|
600 | |
---|
601 | endif |
---|
602 | |
---|
603 | !----------------------------------------------------------------------- |
---|
604 | ! THE VELOCITIES ARE TIME INDEPENDENT |
---|
605 | !----------------------------------------------------------------------- |
---|
606 | |
---|
607 | ! Zonal Velocity |
---|
608 | |
---|
609 | u = u0*(cos(lat)*cos(alpha)+sin(lat)*cos(lon)*sin(alpha)) |
---|
610 | |
---|
611 | ! Meridional Velocity |
---|
612 | |
---|
613 | v = -u0*(sin(lon)*sin(alpha)) |
---|
614 | |
---|
615 | !----------------------------------------------------------------------- |
---|
616 | ! TEMPERATURE IS CONSTANT 300 K |
---|
617 | !----------------------------------------------------------------------- |
---|
618 | |
---|
619 | t = T0 |
---|
620 | |
---|
621 | !----------------------------------------------------------------------- |
---|
622 | ! RHO (density) |
---|
623 | !----------------------------------------------------------------------- |
---|
624 | |
---|
625 | rho = p/(Rd*t) |
---|
626 | |
---|
627 | !----------------------------------------------------------------------- |
---|
628 | ! initialize Q, set to zero |
---|
629 | !----------------------------------------------------------------------- |
---|
630 | |
---|
631 | q = 0.d0 |
---|
632 | |
---|
633 | !----------------------------------------------------------------------- |
---|
634 | ! VERTICAL VELOCITY IS TIME INDEPENDENT |
---|
635 | !----------------------------------------------------------------------- |
---|
636 | |
---|
637 | ! Vertical Velocity - can be changed to vertical pressure velocity by |
---|
638 | ! omega = -(g*p)/(Rd*T0)*w |
---|
639 | |
---|
640 | ! NOTE that if orography-following coordinates are used then the vertical |
---|
641 | ! velocity needs to be translated into the new coordinate system due to |
---|
642 | ! the variation of the height along coordinate surfaces |
---|
643 | ! See section 1.3 and the appendix of the test case document |
---|
644 | |
---|
645 | if (cfv .eq. 0) then |
---|
646 | |
---|
647 | ! if the horizontal velocities do not follow the vertical coordinate |
---|
648 | |
---|
649 | w = 0.d0 |
---|
650 | |
---|
651 | elseif (cfv .eq. 1) then |
---|
652 | |
---|
653 | ! if the horizontal velocities follow hybrid eta coordinates then |
---|
654 | ! the perceived vertical velocity is |
---|
655 | |
---|
656 | call test1_advection_orograph_hybrid_eta_velocity(w) |
---|
657 | |
---|
658 | elseif (cfv .eq. 2) then |
---|
659 | |
---|
660 | ! if the horizontal velocities follow Gal Chen coordinates then |
---|
661 | ! the perceived vertical velocity is |
---|
662 | |
---|
663 | call test1_advection_orograph_Gal_Chen_velocity(w) |
---|
664 | |
---|
665 | ! else |
---|
666 | ! compute your own vertical velocity if other orography-following |
---|
667 | ! vertical coordinate is used |
---|
668 | ! |
---|
669 | endif |
---|
670 | |
---|
671 | |
---|
672 | |
---|
673 | !----------------------------------------------------------------------- |
---|
674 | ! initialize tracers |
---|
675 | !----------------------------------------------------------------------- |
---|
676 | |
---|
677 | ! Tracer 1 - Cloud Layer |
---|
678 | |
---|
679 | r = acos( sin(phip)*sin(lat) + cos(phip)*cos(lat)*cos(lon - lambdap) ) |
---|
680 | |
---|
681 | rz = abs(height - zp1) |
---|
682 | |
---|
683 | if (rz .lt. 0.5d0*dzp1 .and. r .lt. Rp) then |
---|
684 | |
---|
685 | q1 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp1))*(1.d0+cos(pi*r/Rp)) |
---|
686 | |
---|
687 | else |
---|
688 | |
---|
689 | q1 = 0.d0 |
---|
690 | |
---|
691 | endif |
---|
692 | |
---|
693 | rz = abs(height - zp2) |
---|
694 | |
---|
695 | if (rz .lt. 0.5d0*dzp2 .and. r .lt. Rp) then |
---|
696 | |
---|
697 | q2 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp2))*(1.d0+cos(pi*r/Rp)) |
---|
698 | |
---|
699 | else |
---|
700 | |
---|
701 | q2 = 0.d0 |
---|
702 | |
---|
703 | endif |
---|
704 | |
---|
705 | rz = abs(height - zp3) |
---|
706 | |
---|
707 | if (rz .lt. 0.5d0*dzp3 .and. r .lt. Rp) then |
---|
708 | |
---|
709 | q3 = 1.d0 |
---|
710 | |
---|
711 | else |
---|
712 | |
---|
713 | q3 = 0.d0 |
---|
714 | |
---|
715 | endif |
---|
716 | |
---|
717 | q4 = q1 + q2 + q3 |
---|
718 | |
---|
719 | |
---|
720 | CONTAINS |
---|
721 | |
---|
722 | !----------------------------------------------------------------------- |
---|
723 | ! SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY |
---|
724 | ! UNDER HYBRID-ETA COORDINATES |
---|
725 | !----------------------------------------------------------------------- |
---|
726 | |
---|
727 | SUBROUTINE test1_advection_orograph_hybrid_eta_velocity(w) |
---|
728 | IMPLICIT NONE |
---|
729 | real(rstd), intent(out) :: w |
---|
730 | |
---|
731 | real(rstd) :: press, & ! hyam *p0 + hybm *ps |
---|
732 | r, & ! Great Circle Distance |
---|
733 | dzsdx, & ! Part of surface height derivative |
---|
734 | dzsdlambda, & ! Derivative of zs w.r.t lambda |
---|
735 | dzsdphi, & ! Derivative of zs w.r.t phi |
---|
736 | dzdlambda, & ! Derivative of z w.r.t lambda |
---|
737 | dzdphi, & ! Derivative of z w.r.t phi |
---|
738 | dpsdlambda, & ! Derivative of ps w.r.t lambda |
---|
739 | dpsdphi ! Derivative of ps w.r.t phi |
---|
740 | |
---|
741 | ! Calculate pressure and great circle distance to mountain center |
---|
742 | |
---|
743 | press = hyam*p0 + hybm*ps |
---|
744 | |
---|
745 | r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) |
---|
746 | |
---|
747 | ! Derivatives of surface height |
---|
748 | |
---|
749 | if (r .lt. Rm) then |
---|
750 | dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - & |
---|
751 | (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam) |
---|
752 | else |
---|
753 | dzsdx = 0.d0 |
---|
754 | endif |
---|
755 | |
---|
756 | ! Prevent division by zero |
---|
757 | |
---|
758 | if (1.d0-cos(r)**2 .gt. 0.d0) then |
---|
759 | dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) & |
---|
760 | /sqrt(1.d0-cos(r)**2) |
---|
761 | dzsdphi = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) & |
---|
762 | /sqrt(1.d0-cos(r)**2) |
---|
763 | else |
---|
764 | dzsdlambda = 0.d0 |
---|
765 | dzsdphi = 0.d0 |
---|
766 | endif |
---|
767 | |
---|
768 | ! Derivatives of surface pressure |
---|
769 | |
---|
770 | dpsdlambda = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdlambda |
---|
771 | dpsdphi = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdphi |
---|
772 | |
---|
773 | ! Derivatives of coordinate height |
---|
774 | |
---|
775 | dzdlambda = -(Rd*T0/(g*press))*hybm*dpsdlambda |
---|
776 | dzdphi = -(Rd*T0/(g*press))*hybm*dpsdphi |
---|
777 | |
---|
778 | ! Prevent division by zero |
---|
779 | |
---|
780 | if (abs(lat) .lt. pi/2.d0) then |
---|
781 | w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi |
---|
782 | else |
---|
783 | w = 0.d0 |
---|
784 | endif |
---|
785 | |
---|
786 | END SUBROUTINE test1_advection_orograph_hybrid_eta_velocity |
---|
787 | |
---|
788 | |
---|
789 | |
---|
790 | !----------------------------------------------------------------------- |
---|
791 | ! SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY |
---|
792 | ! UNDER GAL-CHEN COORDINATES |
---|
793 | !----------------------------------------------------------------------- |
---|
794 | |
---|
795 | SUBROUTINE test1_advection_orograph_Gal_Chen_velocity(w) |
---|
796 | IMPLICIT NONE |
---|
797 | real(rstd), intent(out) :: w |
---|
798 | |
---|
799 | real(rstd) :: r, & ! Great Circle Distance |
---|
800 | dzsdx, & ! Part of surface height derivative |
---|
801 | dzsdlambda, & ! Derivative of zs w.r.t lambda |
---|
802 | dzsdphi, & ! Derivative of zs w.r.t phi |
---|
803 | dzdlambda, & ! Derivative of z w.r.t lambda |
---|
804 | dzdphi ! Derivative of z w.r.t phi |
---|
805 | |
---|
806 | ! Calculate great circle distance to mountain center |
---|
807 | |
---|
808 | r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) |
---|
809 | |
---|
810 | ! Derivatives of surface height |
---|
811 | |
---|
812 | if (r .lt. Rm) then |
---|
813 | dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - & |
---|
814 | (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam) |
---|
815 | else |
---|
816 | dzsdx = 0.d0 |
---|
817 | endif |
---|
818 | |
---|
819 | ! Prevent division by zero |
---|
820 | |
---|
821 | if (1.d0-cos(r)**2 .gt. 0.d0) then |
---|
822 | dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) & |
---|
823 | /sqrt(1.d0-cos(r)**2) |
---|
824 | dzsdphi = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) & |
---|
825 | /sqrt(1.d0-cos(r)**2) |
---|
826 | else |
---|
827 | dzsdlambda = 0.d0 |
---|
828 | dzsdphi = 0.d0 |
---|
829 | endif |
---|
830 | |
---|
831 | ! Derivatives of coordinate height |
---|
832 | |
---|
833 | dzdlambda = (1.d0-gc/ztop)*dzsdlambda |
---|
834 | dzdphi = (1.d0-gc/ztop)*dzsdphi |
---|
835 | |
---|
836 | ! Prevent division by zero |
---|
837 | |
---|
838 | if (abs(lat) .lt. pi/2.d0) then |
---|
839 | w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi |
---|
840 | else |
---|
841 | w = 0.d0 |
---|
842 | endif |
---|
843 | |
---|
844 | END SUBROUTINE test1_advection_orograph_Gal_Chen_velocity |
---|
845 | |
---|
846 | END SUBROUTINE test1_advection_orography |
---|
847 | |
---|
848 | |
---|
849 | |
---|
850 | |
---|
851 | !========================================================================================== |
---|
852 | ! TEST CASE 2X - IMPACT OF OROGRAPHY ON A NON-ROTATING PLANET |
---|
853 | !========================================================================================== |
---|
854 | ! The tests in section 2-x examine the impact of 3D Schaer-like circular mountain profiles on an |
---|
855 | ! atmosphere at rest (2-0), and on flow fields with wind shear (2-1) and without vertical wind shear (2-2). |
---|
856 | ! A non-rotating planet is used for all configurations. Test 2-0 is conducted on an unscaled regular-size |
---|
857 | ! planet and primarily examines the accuracy of the pressure gradient calculation in a steady-state |
---|
858 | ! hydrostatically-balanced atmosphere at rest. This test is especially appealing for models with |
---|
859 | ! orography-following vertical coordinates. It increases the complexity of test 1-3, that investigated |
---|
860 | ! the impact of the same Schaer-type orographic profile on the accuracy of purely-horizontal passive |
---|
861 | ! tracer advection. |
---|
862 | ! |
---|
863 | ! Tests 2-1 and 2-2 increase the complexity even further since non-zero flow fields are now prescribed |
---|
864 | ! with and without vertical wind shear. In order to trigger non-hydrostatic responses the two tests are |
---|
865 | ! conducted on a reduced-size planet with reduction factor $X=500$ which makes the horizontal and |
---|
866 | ! vertical grid spacing comparable. This test clearly discriminates between non-hydrostatic and hydrostatic |
---|
867 | ! models since the expected response is in the non-hydrostatic regime. Therefore, the flow response is |
---|
868 | ! captured differently by hydrostatic models. |
---|
869 | |
---|
870 | |
---|
871 | |
---|
872 | |
---|
873 | !========================================================================= |
---|
874 | ! Test 2-0: Steady-State Atmosphere at Rest in the Presence of Orography |
---|
875 | !========================================================================= |
---|
876 | SUBROUTINE test2_steady_state_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,u,v,w,t,phis,ps,rho,q) |
---|
877 | |
---|
878 | IMPLICIT NONE |
---|
879 | !----------------------------------------------------------------------- |
---|
880 | ! input/output params parameters at given location |
---|
881 | !----------------------------------------------------------------------- |
---|
882 | |
---|
883 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
884 | lat, & ! Latitude (radians) |
---|
885 | z, & ! Height (m) |
---|
886 | hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint |
---|
887 | hybm ! B coefficient for hybrid-eta coordinate, at model level midpoint |
---|
888 | |
---|
889 | logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used |
---|
890 | ! if set to .true., then the pressure will be computed via the |
---|
891 | ! hybrid coefficients hyam and hybm, they need to be initialized |
---|
892 | ! if set to .false. (for pressure-based models): the pressure is already pre-computed |
---|
893 | ! and is an input value for this routine |
---|
894 | ! for height-based models: pressure will always be computed based on the height and |
---|
895 | ! hybrid_eta is not used |
---|
896 | |
---|
897 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
898 | |
---|
899 | integer, intent(in) :: zcoords ! 0 or 1 see below |
---|
900 | |
---|
901 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
902 | v, & ! Meridional wind (m s^-1) |
---|
903 | w, & ! Vertical Velocity (m s^-1) |
---|
904 | t, & ! Temperature (K) |
---|
905 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
906 | ps, & ! Surface Pressure (Pa) |
---|
907 | rho, & ! density (kg m^-3) |
---|
908 | q ! Specific Humidity (kg/kg) |
---|
909 | |
---|
910 | ! if zcoords = 1, then we use z and output p |
---|
911 | ! if zcoords = 0, then we compute or use p |
---|
912 | ! |
---|
913 | ! In hybrid-eta coords: p = hyam p0 + hybm ps |
---|
914 | ! |
---|
915 | ! The grid-point based initial data are computed in this routine. |
---|
916 | |
---|
917 | !----------------------------------------------------------------------- |
---|
918 | ! test case parameters |
---|
919 | !----------------------------------------------------------------------- |
---|
920 | real(rstd), parameter :: T0 = 300.d0, & ! temperature (K) |
---|
921 | gamma = 0.0065d0, & ! temperature lapse rate (K/m) |
---|
922 | lambdam = 3.d0*pi/2.d0, & ! mountain longitude center point (radians) |
---|
923 | phim = 0.d0, & ! mountain latitude center point (radians) |
---|
924 | h0 = 2000.d0, & ! peak height of the mountain range (m) |
---|
925 | Rm = 3.d0*pi/4.d0, & ! mountain radius (radians) |
---|
926 | zetam = pi/16.d0, & ! mountain oscillation half-width (radians) |
---|
927 | ztop = 12000.d0 ! model top (m) |
---|
928 | |
---|
929 | real(rstd) :: height ! Model level heights (m) |
---|
930 | real(rstd) :: r ! Great circle distance (radians) |
---|
931 | real(rstd) :: zs ! Surface elevation (m) |
---|
932 | real(rstd) :: exponent ! exponent: g/(Rd * gamma) |
---|
933 | real(rstd) :: exponent_rev ! reversed exponent |
---|
934 | |
---|
935 | |
---|
936 | !----------------------------------------------------------------------- |
---|
937 | ! compute exponents |
---|
938 | !----------------------------------------------------------------------- |
---|
939 | exponent = g/(Rd*gamma) |
---|
940 | exponent_rev = 1.d0/exponent |
---|
941 | |
---|
942 | !----------------------------------------------------------------------- |
---|
943 | ! PHIS (surface geopotential) |
---|
944 | !----------------------------------------------------------------------- |
---|
945 | |
---|
946 | r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) |
---|
947 | |
---|
948 | if (r .lt. Rm) then |
---|
949 | |
---|
950 | zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0 ! mountain height |
---|
951 | |
---|
952 | else |
---|
953 | |
---|
954 | zs = 0.d0 |
---|
955 | |
---|
956 | endif |
---|
957 | |
---|
958 | phis = g*zs |
---|
959 | |
---|
960 | !----------------------------------------------------------------------- |
---|
961 | ! PS (surface pressure) |
---|
962 | !----------------------------------------------------------------------- |
---|
963 | |
---|
964 | ps = p0 * (1.d0 - gamma/T0*zs)**exponent |
---|
965 | |
---|
966 | |
---|
967 | !----------------------------------------------------------------------- |
---|
968 | ! HEIGHT AND PRESSURE |
---|
969 | !----------------------------------------------------------------------- |
---|
970 | |
---|
971 | ! Height and pressure are aligned (p = p0 * (1.d0 - gamma/T0*z)**exponent) |
---|
972 | |
---|
973 | if (zcoords .eq. 1) then |
---|
974 | |
---|
975 | height = z |
---|
976 | p = p0 * (1.d0 - gamma/T0*z)**exponent |
---|
977 | |
---|
978 | else |
---|
979 | |
---|
980 | if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients |
---|
981 | height = T0/gamma * (1.d0 - (p/p0)**exponent_rev) ! compute the height at this pressure |
---|
982 | |
---|
983 | endif |
---|
984 | |
---|
985 | !----------------------------------------------------------------------- |
---|
986 | ! THE VELOCITIES ARE ZERO (STATE AT REST) |
---|
987 | !----------------------------------------------------------------------- |
---|
988 | |
---|
989 | ! Zonal Velocity |
---|
990 | |
---|
991 | u = 0.d0 |
---|
992 | |
---|
993 | ! Meridional Velocity |
---|
994 | |
---|
995 | v = 0.d0 |
---|
996 | |
---|
997 | ! Vertical Velocity |
---|
998 | |
---|
999 | w = 0.d0 |
---|
1000 | |
---|
1001 | !----------------------------------------------------------------------- |
---|
1002 | ! TEMPERATURE WITH CONSTANT LAPSE RATE |
---|
1003 | !----------------------------------------------------------------------- |
---|
1004 | |
---|
1005 | t = T0 - gamma*height |
---|
1006 | |
---|
1007 | !----------------------------------------------------------------------- |
---|
1008 | ! RHO (density) |
---|
1009 | !----------------------------------------------------------------------- |
---|
1010 | |
---|
1011 | rho = p/(Rd*t) |
---|
1012 | |
---|
1013 | !----------------------------------------------------------------------- |
---|
1014 | ! initialize Q, set to zero |
---|
1015 | !----------------------------------------------------------------------- |
---|
1016 | |
---|
1017 | q = 0.d0 |
---|
1018 | |
---|
1019 | END SUBROUTINE test2_steady_state_mountain |
---|
1020 | |
---|
1021 | |
---|
1022 | |
---|
1023 | !===================================================================================== |
---|
1024 | ! Tests 2-1 and 2-2: Non-hydrostatic Mountain Waves over a Schaer-type Mountain |
---|
1025 | !===================================================================================== |
---|
1026 | |
---|
1027 | SUBROUTINE test2_schaer_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,shear,u,v,w,t,phis,ps,rho,q) |
---|
1028 | |
---|
1029 | IMPLICIT NONE |
---|
1030 | !----------------------------------------------------------------------- |
---|
1031 | ! input/output params parameters at given location |
---|
1032 | !----------------------------------------------------------------------- |
---|
1033 | |
---|
1034 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
1035 | lat, & ! Latitude (radians) |
---|
1036 | z, & ! Height (m) |
---|
1037 | hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint |
---|
1038 | hybm ! B coefficient for hybrid-eta coordinate, at model level midpoint |
---|
1039 | |
---|
1040 | logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used |
---|
1041 | ! if set to .true., then the pressure will be computed via the |
---|
1042 | ! hybrid coefficients hyam and hybm, they need to be initialized |
---|
1043 | ! if set to .false. (for pressure-based models): the pressure is already pre-computed |
---|
1044 | ! and is an input value for this routine |
---|
1045 | ! for height-based models: pressure will always be computed based on the height and |
---|
1046 | ! hybrid_eta is not used |
---|
1047 | |
---|
1048 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
1049 | |
---|
1050 | |
---|
1051 | integer, intent(in) :: zcoords, & ! 0 or 1 see below |
---|
1052 | shear ! 0 or 1 see below |
---|
1053 | |
---|
1054 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
1055 | v, & ! Meridional wind (m s^-1) |
---|
1056 | w, & ! Vertical Velocity (m s^-1) |
---|
1057 | t, & ! Temperature (K) |
---|
1058 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
1059 | ps, & ! Surface Pressure (Pa) |
---|
1060 | rho, & ! density (kg m^-3) |
---|
1061 | q ! Specific Humidity (kg/kg) |
---|
1062 | |
---|
1063 | ! if zcoords = 1, then we use z and output p |
---|
1064 | ! if zcoords = 0, then we either compute or use p |
---|
1065 | |
---|
1066 | ! if shear = 1, then we use shear flow |
---|
1067 | ! if shear = 0, then we use constant u |
---|
1068 | |
---|
1069 | !----------------------------------------------------------------------- |
---|
1070 | ! test case parameters |
---|
1071 | !----------------------------------------------------------------------- |
---|
1072 | real(rstd), parameter :: X = 500.d0, & ! Reduced Earth reduction factor |
---|
1073 | Om = 0.d0, & ! Rotation Rate of Earth |
---|
1074 | as = a/X, & ! New Radius of small Earth |
---|
1075 | ueq = 20.d0, & ! Reference Velocity |
---|
1076 | Teq = 300.d0, & ! Temperature at Equator |
---|
1077 | Peq = 100000.d0, & ! Reference PS at Equator |
---|
1078 | ztop = 30000.d0, & ! Model Top |
---|
1079 | lambdac = pi/4.d0, & ! Lon of Schar Mountain Center |
---|
1080 | phic = 0.d0, & ! Lat of Schar Mountain Center |
---|
1081 | h0 = 250.d0, & ! Height of Mountain |
---|
1082 | d = 5000.d0, & ! Mountain Half-Width |
---|
1083 | xi = 4000.d0, & ! Mountain Wavelength |
---|
1084 | cs = 0.00025d0 ! Wind Shear (shear=1) |
---|
1085 | |
---|
1086 | real(rstd) :: height ! Model level heights |
---|
1087 | real(rstd) :: sin_tmp, cos_tmp ! Calculation of great circle distance |
---|
1088 | real(rstd) :: r ! Great circle distance |
---|
1089 | real(rstd) :: zs ! Surface height |
---|
1090 | real(rstd) :: c ! Shear |
---|
1091 | |
---|
1092 | !----------------------------------------------------------------------- |
---|
1093 | ! PHIS (surface geopotential) |
---|
1094 | !----------------------------------------------------------------------- |
---|
1095 | |
---|
1096 | sin_tmp = sin(lat) * sin(phic) |
---|
1097 | cos_tmp = cos(lat) * cos(phic) |
---|
1098 | |
---|
1099 | ! great circle distance with 'a/X' |
---|
1100 | |
---|
1101 | r = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac)) |
---|
1102 | zs = h0 * exp(-(r**2)/(d**2))*(cos(pi*r/xi)**2) |
---|
1103 | phis = g*zs |
---|
1104 | |
---|
1105 | !----------------------------------------------------------------------- |
---|
1106 | ! SHEAR FLOW OR CONSTANT FLOW |
---|
1107 | !----------------------------------------------------------------------- |
---|
1108 | |
---|
1109 | if (shear .eq. 1) then |
---|
1110 | |
---|
1111 | c = cs |
---|
1112 | |
---|
1113 | else |
---|
1114 | |
---|
1115 | c = 0.d0 |
---|
1116 | |
---|
1117 | endif |
---|
1118 | |
---|
1119 | !----------------------------------------------------------------------- |
---|
1120 | ! TEMPERATURE |
---|
1121 | !----------------------------------------------------------------------- |
---|
1122 | |
---|
1123 | t = Teq *(1.d0 - (c*ueq*ueq/(g))*(sin(lat)**2) ) |
---|
1124 | |
---|
1125 | !----------------------------------------------------------------------- |
---|
1126 | ! PS (surface pressure) |
---|
1127 | !----------------------------------------------------------------------- |
---|
1128 | |
---|
1129 | ps = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - phis/(Rd*t) ) |
---|
1130 | |
---|
1131 | !----------------------------------------------------------------------- |
---|
1132 | ! HEIGHT AND PRESSURE |
---|
1133 | !----------------------------------------------------------------------- |
---|
1134 | |
---|
1135 | if (zcoords .eq. 1) then |
---|
1136 | |
---|
1137 | height = z |
---|
1138 | p = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - g*height/(Rd*t) ) |
---|
1139 | |
---|
1140 | else |
---|
1141 | if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients |
---|
1142 | height = (Rd*t/(g))*log(peq/p) - (t*ueq*ueq/(2.d0*Teq*g))*(sin(lat)**2) |
---|
1143 | |
---|
1144 | endif |
---|
1145 | |
---|
1146 | !----------------------------------------------------------------------- |
---|
1147 | ! THE VELOCITIES |
---|
1148 | !----------------------------------------------------------------------- |
---|
1149 | |
---|
1150 | ! Zonal Velocity |
---|
1151 | |
---|
1152 | u = ueq * cos(lat) * sqrt( (2.d0*Teq/(t))*c*height + t/(Teq) ) |
---|
1153 | |
---|
1154 | ! Meridional Velocity |
---|
1155 | |
---|
1156 | v = 0.d0 |
---|
1157 | |
---|
1158 | ! Vertical Velocity = Vertical Pressure Velocity = 0 |
---|
1159 | |
---|
1160 | w = 0.d0 |
---|
1161 | |
---|
1162 | !----------------------------------------------------------------------- |
---|
1163 | ! RHO (density) |
---|
1164 | !----------------------------------------------------------------------- |
---|
1165 | |
---|
1166 | rho = p/(Rd*t) |
---|
1167 | |
---|
1168 | !----------------------------------------------------------------------- |
---|
1169 | ! initialize Q, set to zero |
---|
1170 | !----------------------------------------------------------------------- |
---|
1171 | |
---|
1172 | q = 0.d0 |
---|
1173 | |
---|
1174 | END SUBROUTINE test2_schaer_mountain |
---|
1175 | |
---|
1176 | |
---|
1177 | |
---|
1178 | |
---|
1179 | |
---|
1180 | |
---|
1181 | |
---|
1182 | |
---|
1183 | |
---|
1184 | |
---|
1185 | !========================================================================================== |
---|
1186 | ! TEST CASE 3 - GRAVITY WAVES |
---|
1187 | !========================================================================================== |
---|
1188 | |
---|
1189 | ! The non-hydrostatic gravity wave test examines the response of models to short time-scale wavemotion |
---|
1190 | ! triggered by a localized perturbation. The formulation presented in this document is new, |
---|
1191 | ! but is based on previous approaches by Skamarock et al. (JAS 1994), Tomita and Satoh (FDR 2004), and |
---|
1192 | ! Jablonowski et al. (NCAR Tech Report 2008) |
---|
1193 | |
---|
1194 | |
---|
1195 | !========== |
---|
1196 | ! Test 3-1 |
---|
1197 | !========== |
---|
1198 | SUBROUTINE test3_gravity_wave (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) |
---|
1199 | |
---|
1200 | IMPLICIT NONE |
---|
1201 | !----------------------------------------------------------------------- |
---|
1202 | ! input/output params parameters at given location |
---|
1203 | !----------------------------------------------------------------------- |
---|
1204 | |
---|
1205 | real(rstd), intent(in) :: lon, & ! Longitude (radians) |
---|
1206 | lat, & ! Latitude (radians) |
---|
1207 | z ! Height (m) |
---|
1208 | |
---|
1209 | real(rstd), intent(inout) :: p ! Pressure (Pa) |
---|
1210 | |
---|
1211 | |
---|
1212 | integer, intent(in) :: zcoords ! 0 or 1 see below |
---|
1213 | |
---|
1214 | real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) |
---|
1215 | v, & ! Meridional wind (m s^-1) |
---|
1216 | w, & ! Vertical Velocity (m s^-1) |
---|
1217 | t, & ! Temperature (K) |
---|
1218 | phis, & ! Surface Geopotential (m^2 s^-2) |
---|
1219 | ps, & ! Surface Pressure (Pa) |
---|
1220 | rho, & ! density (kg m^-3) |
---|
1221 | q ! Specific Humidity (kg/kg) |
---|
1222 | |
---|
1223 | ! if zcoords = 1, then we use z and output z |
---|
1224 | ! if zcoords = 0, then we use p |
---|
1225 | |
---|
1226 | !----------------------------------------------------------------------- |
---|
1227 | ! test case parameters |
---|
1228 | !----------------------------------------------------------------------- |
---|
1229 | real(rstd), parameter :: X = 125.d0, & ! Reduced Earth reduction factor |
---|
1230 | Om = 0.d0, & ! Rotation Rate of Earth |
---|
1231 | as = a/X, & ! New Radius of small Earth |
---|
1232 | u0 = 20.d0, & ! Reference Velocity |
---|
1233 | Teq = 300.d0, & ! Temperature at Equator |
---|
1234 | Peq = 100000.d0, & ! Reference PS at Equator |
---|
1235 | ztop = 10000.d0, & ! Model Top |
---|
1236 | lambdac = 2.d0*pi/3.d0, & ! Lon of Pert Center |
---|
1237 | d = 5000.d0, & ! Width for Pert |
---|
1238 | phic = 0.d0, & ! Lat of Pert Center |
---|
1239 | delta_theta = 1.d0, & ! Max Amplitude of Pert |
---|
1240 | Lz = 20000.d0, & ! Vertical Wavelength of Pert |
---|
1241 | N = 0.01d0, & ! Brunt-Vaisala frequency |
---|
1242 | N2 = N*N, & ! Brunt-Vaisala frequency Squared |
---|
1243 | bigG = (g*g)/(N2*cp) ! Constant |
---|
1244 | |
---|
1245 | real(rstd) :: height ! Model level height |
---|
1246 | real(rstd) :: sin_tmp, cos_tmp ! Calculation of great circle distance |
---|
1247 | real(rstd) :: r, s ! Shape of perturbation |
---|
1248 | real(rstd) :: TS ! Surface temperature |
---|
1249 | real(rstd) :: t_mean, t_pert ! Mean and pert parts of temperature |
---|
1250 | real(rstd) :: theta_pert ! Pot-temp perturbation |
---|
1251 | |
---|
1252 | !----------------------------------------------------------------------- |
---|
1253 | ! THE VELOCITIES |
---|
1254 | !----------------------------------------------------------------------- |
---|
1255 | |
---|
1256 | ! Zonal Velocity |
---|
1257 | |
---|
1258 | u = u0 * cos(lat) |
---|
1259 | |
---|
1260 | ! Meridional Velocity |
---|
1261 | |
---|
1262 | v = 0.d0 |
---|
1263 | |
---|
1264 | ! Vertical Velocity = Vertical Pressure Velocity = 0 |
---|
1265 | |
---|
1266 | w = 0.d0 |
---|
1267 | |
---|
1268 | !----------------------------------------------------------------------- |
---|
1269 | ! PHIS (surface geopotential) |
---|
1270 | !----------------------------------------------------------------------- |
---|
1271 | |
---|
1272 | phis = 0.d0 |
---|
1273 | |
---|
1274 | !----------------------------------------------------------------------- |
---|
1275 | ! SURFACE TEMPERATURE |
---|
1276 | !----------------------------------------------------------------------- |
---|
1277 | |
---|
1278 | TS = bigG + (Teq-bigG)*exp( -(u0*N2/(4.d0*g*g))*(u0+2.d0*om*as)*(cos(2.d0*lat)-1.d0) ) |
---|
1279 | |
---|
1280 | !----------------------------------------------------------------------- |
---|
1281 | ! PS (surface pressure) |
---|
1282 | !----------------------------------------------------------------------- |
---|
1283 | |
---|
1284 | ps = peq*exp( (u0/(4.0*bigG*Rd))*(u0+2.0*Om*as)*(cos(2.0*lat)-1.0) ) & |
---|
1285 | * (TS/Teq)**(cp/Rd) |
---|
1286 | |
---|
1287 | !----------------------------------------------------------------------- |
---|
1288 | ! HEIGHT AND PRESSURE AND MEAN TEMPERATURE |
---|
1289 | !----------------------------------------------------------------------- |
---|
1290 | |
---|
1291 | if (zcoords .eq. 1) then |
---|
1292 | |
---|
1293 | height = z |
---|
1294 | p = ps*( (bigG/TS)*exp(-N2*height/g)+1.d0 - (bigG/TS) )**(cp/Rd) |
---|
1295 | |
---|
1296 | else |
---|
1297 | |
---|
1298 | height = (-g/N2)*log( (TS/bigG)*( (p/ps)**(Rd/cp) - 1.d0 ) + 1.d0 ) |
---|
1299 | |
---|
1300 | endif |
---|
1301 | |
---|
1302 | t_mean = bigG*(1.d0 - exp(N2*height/g))+ TS*exp(N2*height/g) |
---|
1303 | |
---|
1304 | !----------------------------------------------------------------------- |
---|
1305 | ! rho (density), unperturbed using the background temperature t_mean |
---|
1306 | !----------------------------------------------------------------------- |
---|
1307 | |
---|
1308 | !*********** |
---|
1309 | ! change in version 3: density is now initialized with unperturbed background temperature, |
---|
1310 | ! temperature perturbation is added afterwards |
---|
1311 | !*********** |
---|
1312 | rho = p/(Rd*t_mean) |
---|
1313 | |
---|
1314 | !----------------------------------------------------------------------- |
---|
1315 | ! POTENTIAL TEMPERATURE PERTURBATION, |
---|
1316 | ! here: converted to temperature and added to the temperature field |
---|
1317 | ! models with a prognostic potential temperature field can utilize |
---|
1318 | ! the potential temperature perturbation theta_pert directly and add it |
---|
1319 | ! to the background theta field (not included here) |
---|
1320 | !----------------------------------------------------------------------- |
---|
1321 | |
---|
1322 | sin_tmp = sin(lat) * sin(phic) |
---|
1323 | cos_tmp = cos(lat) * cos(phic) |
---|
1324 | |
---|
1325 | ! great circle distance with 'a/X' |
---|
1326 | |
---|
1327 | r = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac)) |
---|
1328 | |
---|
1329 | s = (d**2)/(d**2 + r**2) |
---|
1330 | |
---|
1331 | theta_pert = delta_theta*s*sin(2.d0*pi*height/Lz) |
---|
1332 | |
---|
1333 | t_pert = theta_pert*(p/p0)**(Rd/cp) |
---|
1334 | |
---|
1335 | t = t_mean + t_pert |
---|
1336 | |
---|
1337 | !----------------------------------------------------------------------- |
---|
1338 | ! initialize Q, set to zero |
---|
1339 | !----------------------------------------------------------------------- |
---|
1340 | |
---|
1341 | q = 0.d0 |
---|
1342 | |
---|
1343 | END SUBROUTINE test3_gravity_wave |
---|
1344 | |
---|
1345 | |
---|
1346 | END MODULE dcmip_initial_conditions_test_1_2_3 |
---|
1347 | |
---|