Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

dcmip_initial_conditions_test_1_2_3_v5.f90

Last change on this file was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC

File size: 47.3 KB

Line
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.d0pia)/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 ! 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.d0pitime/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)/(bsptop) ) - exp( (p-p0)/(bsptop)) - exp( (ptop-p)/(bs*ptop))
189
190	! Zonal Velocity
191	!********
192	! change in version 5: ud
193	!********
194
195	ud = (omega0a)/(bsptop) * cos(lonp) * (cos(lat)*2.0) cos(2.0pitime/tau) * &
196	( - exp( (p-p0)/(bsptop)) + exp( (ptop-p)/(bsptop)) )
197
198	u = k0sin(lonp)sin(lonp)sin(2.d0lat)cos(pitime/tau) + u0*cos(lat) + ud
199
200	! Meridional Velocity
201
202	v = k0sin(2.d0lonp)cos(lat)cos(pi*time/tau)
203
204	! Vertical Velocity - can be changed to vertical pressure velocity by
205	! omega = -(gp)/(RdT0)*w
206	!
207	!********
208	! change in version 4: cos(2.0pitime/tau) is now used instead of cos(pi*time/tau)
209	!********
210	!********
211	! change in version 5: shape function in w
212	!********
213	w = -((RdT0)/(gp))omega0sin(lonp)cos(lat)cos(2.0pitime/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(pid1)) + 0.5d0 (1.d0 + cos(pi*d2))
266
267	! Tracer 2 - Correlated Cosine Bells
268
269	q2 = 0.9d0 - 0.8d0q1*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.25ZZ) .or. height .lt. (z0-1.25ZZ)) 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) * (aw0pi)/(Kztop) cos(lat)sin(Klat)cos(piheight/ztop)cos(pitime/tau)
436
437	! Vertical Velocity - can be changed to vertical pressure velocity by
438	! omega = -grhow
439
440	w = (rho0/rho) (w0/K)(-2.d0sin(Klat)sin(lat) + Kcos(lat)cos(Klat)) &
441	sin(piheight/ztop)cos(pitime/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.d0pi(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.d0pia/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(pir/Rm))cos(pir/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 = hyamp0 + hybmps ! 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 = -(gp)/(RdT0)*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.d0pirz/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.d0pirz/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 = hyamp0 + hybmps
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 = -h0pi/(2.d0Rm)sin(pir/Rm)cos(pir/zetam)**2 - &
751	(h0pi/zetam)(1.d0+cos(pir/Rm))cos(pir/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 = -(gp0/(RdT0))exp(-gzs/(RdT0))dzsdlambda
771	dpsdphi = -(gp0/(RdT0))exp(-gzs/(RdT0))dzsdphi
772
773	! Derivatives of coordinate height
774
775	dzdlambda = -(RdT0/(gpress))hybmdpsdlambda
776	dzdphi = -(RdT0/(gpress))hybmdpsdphi
777
778	! Prevent division by zero
779
780	if (abs(lat) .lt. pi/2.d0) then
781	w = - (u/(acos(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 = -h0pi/(2.d0Rm)sin(pir/Rm)cos(pir/zetam)**2 - &
814	(h0pi/zetam)(1.d0+cos(pir/Rm))cos(pir/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/(acos(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(pir/Rm))cos(pir/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/T0zs)*exponent
965
966
967	!-----------------------------------------------------------------------
968	! HEIGHT AND PRESSURE
969	!-----------------------------------------------------------------------
970
971	! Height and pressure are aligned (p = p0 * (1.d0 - gamma/T0z)*exponent)
972
973	if (zcoords .eq. 1) then
974
975	height = z
976	p = p0 * (1.d0 - gamma/T0z)*exponent
977
978	else
979
980	if (hybrid_eta) p = hyamp0 + hybmps ! 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(-(r2)/(d2))(cos(pir/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 - (cuequeq/(g))(sin(lat)**2) )
1124
1125	!-----------------------------------------------------------------------
1126	! PS (surface pressure)
1127	!-----------------------------------------------------------------------
1128
1129	ps = peqexp( -(uequeq/(2.d0RdTeq))(sin(lat)2) - phis/(Rdt) )
1130
1131	!-----------------------------------------------------------------------
1132	! HEIGHT AND PRESSURE
1133	!-----------------------------------------------------------------------
1134
1135	if (zcoords .eq. 1) then
1136
1137	height = z
1138	p = peqexp( -(uequeq/(2.d0RdTeq))(sin(lat)2) - gheight/(Rd*t) )
1139
1140	else
1141	if (hybrid_eta) p = hyamp0 + hybmps ! compute the pressure based on the surface pressure and hybrid coefficients
1142	height = (Rdt/(g))log(peq/p) - (tuequeq/(2.d0Teqg))(sin(lat)*2)
1143
1144	endif
1145
1146	!-----------------------------------------------------------------------
1147	! THE VELOCITIES
1148	!-----------------------------------------------------------------------
1149
1150	! Zonal Velocity
1151
1152	u = ueq * cos(lat) * sqrt( (2.d0Teq/(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 (X,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	X ! Reduced Earth reduction factor (DCMIP value = 125)
1208
1209	real(rstd), intent(inout) :: p, & ! Pressure (Pa)
1210	z ! Height (m)
1211
1212
1213	integer, intent(in) :: zcoords ! 0 or 1 see below
1214
1215	real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1)
1216	v, & ! Meridional wind (m s^-1)
1217	w, & ! Vertical Velocity (m s^-1)
1218	t, & ! Temperature (K)
1219	phis, & ! Surface Geopotential (m^2 s^-2)
1220	ps, & ! Surface Pressure (Pa)
1221	rho, & ! density (kg m^-3)
1222	q ! Specific Humidity (kg/kg)
1223
1224	! if zcoords = 1, then we use z and output z
1225	! if zcoords = 0, then we use p
1226
1227	!-----------------------------------------------------------------------
1228	! test case parameters
1229	!-----------------------------------------------------------------------
1230	real(rstd), parameter :: &
1231	Om = 0.d0, & ! Rotation Rate of Earth
1232	u0 = 20.d0, & ! Reference Velocity
1233	! u0 = 0.d0, & ! FIXME : no zonal wind for NH tests
1234	Teq = 300.d0, & ! Temperature at Equator
1235	Peq = 100000.d0, & ! Reference PS at Equator
1236	ztop = 10000.d0, & ! Model Top
1237	lambdac = 2.d0*pi/3.d0, & ! Lon of Pert Center
1238	d = 5000.d0, & ! Width for Pert
1239	phic = 0.d0, & ! Lat of Pert Center
1240	delta_theta = 1.d0, & ! Max Amplitude of Pert
1241	Lz = 20000.d0, & ! Vertical Wavelength of Pert
1242	N = 0.01d0, & ! Brunt-Vaisala frequency
1243	N2 = N*N, & ! Brunt-Vaisala frequency Squared
1244	bigG = (gg)/(N2cp) ! Constant
1245
1246	real(rstd) :: as ! New Radius of small Earth
1247
1248	real(rstd) :: height ! Model level height
1249	real(rstd) :: sin_tmp, cos_tmp ! Calculation of great circle distance
1250	real(rstd) :: r, s ! Shape of perturbation
1251	real(rstd) :: TS ! Surface temperature
1252	real(rstd) :: t_mean, t_pert ! Mean and pert parts of temperature
1253	real(rstd) :: theta_pert ! Pot-temp perturbation
1254
1255	as = a/X
1256
1257	!-----------------------------------------------------------------------
1258	! THE VELOCITIES
1259	!-----------------------------------------------------------------------
1260
1261	! Zonal Velocity
1262
1263	u = u0 * cos(lat)
1264
1265	! Meridional Velocity
1266
1267	v = 0.d0
1268
1269	! Vertical Velocity = Vertical Pressure Velocity = 0
1270
1271	w = 0.d0
1272
1273	!-----------------------------------------------------------------------
1274	! PHIS (surface geopotential)
1275	!-----------------------------------------------------------------------
1276
1277	phis = 0.d0
1278
1279	!-----------------------------------------------------------------------
1280	! SURFACE TEMPERATURE
1281	!-----------------------------------------------------------------------
1282
1283	TS = bigG + (Teq-bigG)exp( -(u0N2/(4.d0gg))(u0+2.d0omas)(cos(2.d0*lat)-1.d0) )
1284
1285	!-----------------------------------------------------------------------
1286	! PS (surface pressure)
1287	!-----------------------------------------------------------------------
1288
1289	ps = peqexp( (u0/(4.0bigGRd))(u0+2.0Omas)(cos(2.0lat)-1.0) ) &
1290	* (TS/Teq)**(cp/Rd)
1291
1292	!-----------------------------------------------------------------------
1293	! HEIGHT AND PRESSURE AND MEAN TEMPERATURE
1294	!-----------------------------------------------------------------------
1295
1296	if (zcoords .eq. 1) then
1297
1298	height = z
1299	p = ps( (bigG/TS)exp(-N2height/g)+1.d0 - (bigG/TS) )*(cp/Rd)
1300
1301	else
1302
1303	height = (-g/N2)log( (TS/bigG)( (p/ps)**(Rd/cp) - 1.d0 ) + 1.d0 )
1304	z = height ! modified : initialize z when pressure-based
1305
1306	endif
1307
1308	t_mean = bigG(1.d0 - exp(N2height/g))+ TSexp(N2height/g)
1309
1310	!-----------------------------------------------------------------------
1311	! rho (density), unperturbed using the background temperature t_mean
1312	!-----------------------------------------------------------------------
1313
1314	!***********
1315	! change in version 3: density is now initialized with unperturbed background temperature,
1316	! temperature perturbation is added afterwards
1317	!***********
1318	rho = p/(Rd*t_mean)
1319
1320	!-----------------------------------------------------------------------
1321	! POTENTIAL TEMPERATURE PERTURBATION,
1322	! here: converted to temperature and added to the temperature field
1323	! models with a prognostic potential temperature field can utilize
1324	! the potential temperature perturbation theta_pert directly and add it
1325	! to the background theta field (not included here)
1326	!-----------------------------------------------------------------------
1327
1328	sin_tmp = sin(lat) * sin(phic)
1329	cos_tmp = cos(lat) * cos(phic)
1330
1331	! great circle distance with 'a/X'
1332
1333	r = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac))
1334
1335	s = (d2)/(d2 + r**2)
1336
1337	theta_pert = delta_thetassin(2.d0piheight/Lz)
1338
1339	t_pert = theta_pert(p/p0)*(Rd/cp)
1340
1341	t = t_mean + t_pert
1342
1343	!-----------------------------------------------------------------------
1344	! initialize Q, set to zero
1345	!-----------------------------------------------------------------------
1346
1347	q = 0.d0
1348
1349	END SUBROUTINE test3_gravity_wave
1350
1351
1352	END MODULE dcmip_initial_conditions_test_1_2_3
1353

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: codes/icosagcm/trunk/src/dcmip/dcmip_initial_conditions_test_1_2_3_v5.f90

Download in other formats: