Context Navigation

source: codes/icosagcm/trunk/src/physics_dcmip.f90 @ 213

Last change on this file since 213 was 213, checked in by dubos, 10 years ago
New dyn/phys interface - halo points not passed to physics any more (cleanup follows)
File size: 27.8 KB

Line
1	MODULE physics_dcmip_mod
2	USE ICOSA
3	PRIVATE
4
5	INTEGER,SAVE :: testcase
6	!$OMP THREADPRIVATE(testcase)
7
8	TYPE(t_field),POINTER :: f_out_i(:)
9	REAL(rstd),POINTER :: out_i(:,:)
10
11	TYPE(t_field),POINTER :: f_precl(:)
12	REAL(rstd),ALLOCATABLE :: precl_packed(:)
13
14	PUBLIC :: compute_phys_wrap, init_physics, &
15	init_physics_new, full_physics, write_physics
16
17	CONTAINS
18
19	!-------------------------- New interface ----------------------
20
21	SUBROUTINE init_physics_new
22	USE physics_interface_mod
23	IMPLICIT NONE
24	INTEGER :: ngrid
25	testcase=1 ! OK for 4.2 (moist baroclinic instability)
26	CALL getin("dcmip_physics",testcase)
27	ngrid = physics_inout%ngrid
28	! Input
29	ALLOCATE(physics_inout%Ai(ngrid))
30	ALLOCATE(physics_inout%lon(ngrid))
31	ALLOCATE(physics_inout%lat(ngrid))
32	ALLOCATE(physics_inout%phis(ngrid))
33	ALLOCATE(physics_inout%p(ngrid,llm+1))
34	ALLOCATE(physics_inout%Temp(ngrid,llm))
35	ALLOCATE(physics_inout%ulon(ngrid,llm))
36	ALLOCATE(physics_inout%ulat(ngrid,llm))
37	ALLOCATE(physics_inout%q(ngrid,llm,nqtot))
38	! Output (tendencies)
39	ALLOCATE(physics_inout%dTemp(ngrid,llm))
40	ALLOCATE(physics_inout%dulon(ngrid,llm))
41	ALLOCATE(physics_inout%dulat(ngrid,llm))
42	ALLOCATE(physics_inout%dq(ngrid,llm,nqtot))
43	! Physics-specific data
44	ALLOCATE(precl_packed(ngrid))
45	CALL allocate_field(f_precl, field_t,type_real)
46
47	PRINT *, 'init_physics_new', SIZE(physics_inout%Ai)
48	END SUBROUTINE init_physics_new
49
50	SUBROUTINE full_physics
51	USE physics_interface_mod
52	CALL compute_physics(physics_inout%ngrid, physics_inout%dt_phys, &
53	physics_inout%lat, physics_inout%p, physics_inout%Temp, &
54	physics_inout%ulon, physics_inout%ulat, physics_inout%q(:,:,1), &
55	physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat, &
56	physics_inout%dq(:,:,1), precl_packed)
57	END SUBROUTINE full_physics
58
59	SUBROUTINE write_physics
60	USE output_field_mod
61	USE physics_interface_mod
62	CALL unpack_field(f_precl, precl_packed)
63	CALL output_field("precl",f_precl)
64	END SUBROUTINE write_physics
65
66	!------------------------ Draft interface -----------------------
67
68	SUBROUTINE init_physics
69	USE physics_interface_mod
70	IMPLICIT NONE
71	testcase=1 ! OK for 4.2 (moist baroclinic instability)
72	CALL getin("dcmip_physics",testcase)
73	nb_extra_physics_2D=1 ! precl
74	nb_extra_physics_3D=0
75	END SUBROUTINE init_physics
76
77	SUBROUTINE compute_phys_wrap(args)
78	USE physics_interface_mod
79	TYPE(t_physics_inout) :: args
80	CALL compute_physics(args%ngrid, args%dt_phys, args%lat, &
81	args%p, args%Temp, args%ulon, args%ulat, args%q(:,:,1), &
82	args%dTemp, args%dulon, args%dulat, args%dq(:,:,1), args%extra_2D(:,1))
83	END SUBROUTINE compute_phys_wrap
84
85	!------------------ Interface-independent wrapper ---------------------------
86
87	SUBROUTINE compute_physics(ngrid,dt_phys,lat, p,Temp,u,v,q, dTemp,du,dv,dq, precl)
88	USE icosa
89	IMPLICIT NONE
90	INTEGER :: ngrid
91	REAL(rstd) :: lat(ngrid)
92	REAL(rstd) :: ps(ngrid)
93	REAL(rstd) :: precl(ngrid)
94	! arguments with bottom-up indexing (DYNAMICO)
95	REAL(rstd) :: p(ngrid,llm+1)
96	REAL(rstd) :: Temp(ngrid,llm)
97	REAL(rstd) :: u(ngrid,llm)
98	REAL(rstd) :: v(ngrid,llm)
99	REAL(rstd) :: q(ngrid,llm)
100	REAL(rstd) :: dTemp(ngrid,llm)
101	REAL(rstd) :: du(ngrid,llm)
102	REAL(rstd) :: dv(ngrid,llm)
103	REAL(rstd) :: dq(ngrid,llm)
104	! local arrays with top-down vertical indexing (DCMIP)
105	REAL(rstd) :: pint(ngrid,llm+1)
106	REAL(rstd) :: pmid(ngrid,llm)
107	REAL(rstd) :: pdel(ngrid,llm)
108	REAL(rstd) :: Tfi(ngrid,llm)
109	REAL(rstd) :: ufi(ngrid,llm)
110	REAL(rstd) :: vfi(ngrid,llm)
111	REAL(rstd) :: qfi(ngrid,llm)
112
113	INTEGER :: l,ll,ij
114	REAL(rstd) :: dt_phys, inv_dt
115
116	! prepare input fields and mirror vertical index
117	ps(:) = p(:,1) ! surface pressure
118
119	DO l=1,llm+1
120	DO ij=1,ngrid
121	pint(ij,l)=p(ij,llm+2-l)
122	ENDDO
123	ENDDO
124
125	DO l=1,llm
126	ll=llm+1-l
127	DO ij=1,ngrid
128	pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers
129	pdel(ij,l)=pint(ij,l+1)-pint(ij,l) ! Pressure difference between two layers
130	ufi(ij,l)=u(ij,ll)
131	vfi(ij,l)=v(ij,ll)
132	qfi(ij,l)=q(ij,ll)
133	Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l))
134	ENDDO
135	ENDDO
136
137	CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, testcase)
138
139	! Mirror vertical index and compute tendencies
140	inv_dt = 1./dt_phys
141	DO l=1,llm
142	ll=llm+1-l
143	DO ij=1,ngrid
144	dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)(1+0.608qfi(ij,ll)) - Temp(ij,l) )
145	du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l))
146	dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l))
147	dq(ij,l) = inv_dt * (qfi(ij,ll) - q(ij,l))
148	ENDDO
149	ENDDO
150
151	END SUBROUTINE compute_physics
152
153	END MODULE physics_dcmip_mod
154
155	subroutine simple_physics (pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test)
156	!-----------------------------------------------------------------------
157	!
158	! Purpose: Simple Physics Package
159	!
160	! Author: K. A. Reed (University of Michigan, kareed@umich.edu)
161	! version 5
162	! July/8/2012
163	!
164	! Change log:
165	! v2: removal of some NCAR CAM-specific 'use' associations
166	! v3: corrected precl(i) computation, the precipitation rate is now computed via a vertical integral, the previous single-level computation in v2 was a bug
167	! v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing the temperature variable: '-t(i,1)'
168	! v4: modified and enhanced parameter list to make the routine truly standalone, the number of columns and vertical levels have been added: pcols, pver
169	! v4: 'ncol' has been removed, 'pcols' is used instead
170	! v5: the sea surface temperature (SST) field Tsurf is now an array, the SST now depends on the latitude
171	! v5: addition of the latitude array 'lat' and the flag 'test' in the parameter list
172	! if test = 0: constant SST is used, correct setting for the tropical cyclone test case 5-1
173	! if test = 1: newly added latitude-dependent SST is used, correct setting for the moist baroclinic wave test with simple-physics (test 4-3)
174	!
175	! Description: Includes large-scale precipitation, surface fluxes and
176	! boundary-leyer mixing. The processes are time-split
177	! in that order. A partially implicit formulation is
178	! used to foster numerical stability.
179	! The routine assumes that the model levels are ordered
180	! in a top-down approach, e.g. level 1 denotes the uppermost
181	! full model level.
182	!
183	! This routine is based on an implementation which was
184	! developed for the NCAR Community Atmosphere Model (CAM).
185	! Adjustments for other models will be necessary.
186	!
187	! The routine provides both updates of the state variables
188	! u, v, T, q (these are local copies of u,v,T,q within this physics
189	! routine) and also collects their time tendencies.
190	! The latter might be used to couple the physics and dynamics
191	! in a process-split way. For a time-split coupling, the final
192	! state should be given to the dynamical core for the next time step.
193	! Test: 0 = Reed and Jablonowski (2011) tropical cyclone test case (test 5-1)
194	! 1 = Moist baroclinic instability test (test 4-3)
195	!
196	!
197	! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone
198	! simulations of intermediate complexity: A test case for AGCMs,
199	! J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099
200	!-----------------------------------------------------------------------
201	! use physics_types , only: physics_dme_adjust ! This is for CESM/CAM
202	! use cam_diagnostics, only: diag_phys_writeout ! This is for CESM/CAM
203
204	implicit none
205
206	integer, parameter :: r8 = selected_real_kind(12)
207
208	!
209	! Input arguments - MODEL DEPENDENT
210	!
211	integer, intent(in) :: pcols ! Set number of atmospheric columns
212	integer, intent(in) :: pver ! Set number of model levels
213	real(r8), intent(in) :: dtime ! Set model physics timestep
214	real(r8), intent(in) :: lat(pcols) ! Latitude
215	integer, intent(in) :: test ! Test number
216
217	!
218	! Input/Output arguments
219	!
220	! pcols is the maximum number of vertical columns per 'chunk' of atmosphere
221	!
222	real(r8), intent(inout) :: t(pcols,pver) ! Temperature at full-model level (K)
223	real(r8), intent(inout) :: q(pcols,pver) ! Specific Humidity at full-model level (kg/kg)
224	real(r8), intent(inout) :: u(pcols,pver) ! Zonal wind at full-model level (m/s)
225	real(r8), intent(inout) :: v(pcols,pver) ! Meridional wind at full-model level (m/s)
226	real(r8), intent(inout) :: pmid(pcols,pver) ! Pressure is full-model level (Pa)
227	real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa)
228	real(r8), intent(inout) :: pdel(pcols,pver) ! Layer thickness (Pa)
229	real(r8), intent(in) :: rpdel(pcols,pver) ! Reciprocal of layer thickness (1/Pa)
230	real(r8), intent(inout) :: ps(pcols) ! Surface Pressue (Pa)
231
232	!
233	! Output arguments
234	!
235	real(r8), intent(out) :: precl(pcols) ! Precipitation rate (m_water / s)
236
237	!
238	!---------------------------Local workspace-----------------------------
239	!
240
241	! Integers for loops
242
243	integer i,k ! Longitude, level indices
244
245	! Physical Constants - Many of these may be model dependent
246
247	real(r8) gravit ! Gravity
248	real(r8) rair ! Gas constant for dry air
249	real(r8) cpair ! Specific heat of dry air
250	real(r8) latvap ! Latent heat of vaporization
251	real(r8) rh2o ! Gas constant for water vapor
252	real(r8) epsilo ! Ratio of gas constant for dry air to that for vapor
253	real(r8) zvir ! Constant for virtual temp. calc. =(rh2o/rair) - 1
254	real(r8) a ! Reference Earth's Radius (m)
255	real(r8) omega ! Reference rotation rate of the Earth (s^-1)
256	real(r8) pi ! pi
257
258	! Simple Physics Specific Constants
259
260	!++++++++
261	real(r8) Tsurf(pcols) ! Sea Surface Temperature (constant for tropical cyclone)
262	!++++++++ Tsurf needs to be dependent on latitude for the
263	! moist baroclinic wave test 4-3 with simple-physics, adjust
264
265	real(r8) SST_tc ! Sea Surface Temperature for tropical cyclone test
266	real(r8) T0 ! Control temp for calculation of qsat
267	real(r8) e0 ! Saturation vapor pressure at T0 for calculation of qsat
268	real(r8) rhow ! Density of Liquid Water
269	real(r8) p0 ! Constant for calculation of potential temperature
270	real(r8) Cd0 ! Constant for calculating Cd from Smith and Vogl 2008
271	real(r8) Cd1 ! Constant for calculating Cd from Smith and Vogl 2008
272	real(r8) Cm ! Constant for calculating Cd from Smith and Vogl 2008
273	real(r8) v20 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
274	real(r8) C ! Drag coefficient for sensible heat and evaporation
275	real(r8) T00 ! Horizontal mean T at surface for moist baro test
276	real(r8) u0 ! Zonal wind constant for moist baro test
277	real(r8) latw ! halfwidth for for baro test
278	real(r8) eta0 ! Center of jets (hybrid) for baro test
279	real(r8) etav ! Auxiliary variable for baro test
280	real(r8) q0 ! Maximum specific humidity for baro test
281
282	! Physics Tendency Arrays
283	real(r8) dtdt(pcols,pver) ! Temperature tendency
284	real(r8) dqdt(pcols,pver) ! Specific humidity tendency
285	real(r8) dudt(pcols,pver) ! Zonal wind tendency
286	real(r8) dvdt(pcols,pver) ! Meridional wind tendency
287
288	! Temporary variables for tendency calculations
289
290	real(r8) tmp ! Temporary
291	real(r8) qsat ! Saturation vapor pressure
292	real(r8) qsats ! Saturation vapor pressure of SST
293
294	! Variables for Boundary Layer Calculation
295
296	real(r8) wind(pcols) ! Magnitude of Wind
297	real(r8) Cd(pcols) ! Drag coefficient for momentum
298	real(r8) Km(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations
299	real(r8) Ke(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations
300	real(r8) rho ! Density at lower/upper interface
301	real(r8) za(pcols) ! Heights at midpoints of first model level
302	real(r8) dlnpint ! Used for calculation of heights
303	real(r8) pbltop ! Top of boundary layer
304	real(r8) pblconst ! Constant for the calculation of the decay of diffusivity
305	real(r8) CA(pcols,pver) ! Matrix Coefficents for PBL Scheme
306	real(r8) CC(pcols,pver) ! Matrix Coefficents for PBL Scheme
307	real(r8) CE(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
308	real(r8) CAm(pcols,pver) ! Matrix Coefficents for PBL Scheme
309	real(r8) CCm(pcols,pver) ! Matrix Coefficents for PBL Scheme
310	real(r8) CEm(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
311	real(r8) CFu(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
312	real(r8) CFv(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
313	real(r8) CFt(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
314	real(r8) CFq(pcols,pver+1) ! Matrix Coefficents for PBL Scheme
315
316
317	! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to
318	! conserve the mass of the dry air
319
320	real(r8) qini(pcols,pver) ! Initial specific humidity
321
322	!===============================================================================
323	!
324	! Physical Constants - MAY BE MODEL DEPENDENT
325	!
326	!===============================================================================
327	gravit = 9.80616_r8 ! Gravity (9.80616 m/s^2)
328	rair = 287.0_r8 ! Gas constant for dry air: 287 J/(kg K)
329	cpair = 1.0045e3_r8 ! Specific heat of dry air: here we use 1004.5 J/(kg K)
330	latvap = 2.5e6_r8 ! Latent heat of vaporization (J/kg)
331	rh2o = 461.5_r8 ! Gas constant for water vapor: 461.5 J/(kg K)
332	epsilo = rair/rh2o ! Ratio of gas constant for dry air to that for vapor
333	zvir = (rh2o/rair) - 1._r8 ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608
334	a = 6371220.0_r8 ! Reference Earth's Radius (m)
335	omega = 7.29212d-5 ! Reference rotation rate of the Earth (s^-1)
336	pi = 4._r8*atan(1._r8) ! pi
337
338	!===============================================================================
339	!
340	! Local Constants for Simple Physics
341	!
342	!===============================================================================
343	C = 0.0011_r8 ! From Smith and Vogl 2008
344	SST_tc = 302.15_r8 ! Constant Value for SST for tropical cyclone test
345	T0 = 273.16_r8 ! control temp for calculation of qsat
346	e0 = 610.78_r8 ! saturation vapor pressure at T0 for calculation of qsat
347	rhow = 1000.0_r8 ! Density of Liquid Water
348	Cd0 = 0.0007_r8 ! Constant for Cd calc. Smith and Vogl 2008
349	Cd1 = 0.000065_r8 ! Constant for Cd calc. Smith and Vogl 2008
350	Cm = 0.002_r8 ! Constant for Cd calc. Smith and Vogl 2008
351	v20 = 20.0_r8 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
352	p0 = 100000.0_r8 ! Constant for potential temp calculation
353	pbltop = 85000._r8 ! Top of boundary layer
354	pblconst = 10000._r8 ! Constant for the calculation of the decay of diffusivity
355	T00 = 288.0_r8 ! Horizontal mean T at surface for moist baro test
356	u0 = 35.0_r8 ! Zonal wind constant for moist baro test
357	latw = 2.0_r8*pi/9.0_r8 ! Halfwidth for for baro test
358	eta0 = 0.252_r8 ! Center of jets (hybrid) for baro test
359	etav = (1._r8-eta0)0.5_r8pi ! Auxiliary variable for baro test
360	q0 = 0.021 ! Maximum specific humidity for baro test
361
362	!===============================================================================
363	!
364	! Definition of local arrays
365	!
366	!===============================================================================
367	!
368	! Calculate hydrostatic height za of the lowest model level
369	!
370	do i=1,pcols
371	dlnpint = log(ps(i)) - log(pint(i,pver)) ! ps(i) is identical to pint(i,pver+1), note: this is the correct sign (corrects typo in JAMES paper)
372	za(i) = rair/gravitt(i,pver)(1._r8+zvirq(i,pver))0.5_r8*dlnpint
373	end do
374	!
375	! Set Initial Specific Humidity - For dry mass adjustment at the end
376	!
377	qini(:pcols,:pver) = q(:pcols,:pver)
378	!--------------------------------------------------------------
379	! Set Sea Surface Temperature (constant for tropical cyclone)
380	! Tsurf needs to be dependent on latitude for the
381	! moist baroclinic wave test 4-3 with simple-physics
382	!--------------------------------------------------------------
383	if (test .eq. 1) then ! moist baroclinic wave with simple-physics
384	do i=1,pcols
385	Tsurf(i) = (T00 + piu0/rair 1.5_r8 * sin(etav) * (cos(etav))*0.5_r8 &
386	((-2._r8(sin(lat(i)))6 ((cos(lat(i)))*2 + 1._r8/3._r8) + 10._r8/63._r8) &
387	u0 * (cos(etav))**1.5_r8 + &
388	(8._r8/5._r8(cos(lat(i)))3 ((sin(lat(i)))*2 + 2._r8/3._r8) - pi/4._r8)aomega0.5_r8 ))/ &
389	(1._r8+zvirq0exp(-(lat(i)/latw)**4))
390
391	end do
392	else
393	do i=1,pcols ! constant SST for the tropical cyclone test case
394	Tsurf(i) = SST_tc
395	end do
396	end if
397
398	!===============================================================================
399	!
400	! Set initial physics time tendencies and precipitation field to zero
401	!
402	!===============================================================================
403	dtdt(:pcols,:pver) = 0._r8 ! initialize temperature tendency with zero
404	dqdt(:pcols,:pver) = 0._r8 ! initialize specific humidity tendency with zero
405	dudt(:pcols,:pver) = 0._r8 ! initialize zonal wind tendency with zero
406	dvdt(:pcols,:pver) = 0._r8 ! initialize meridional wind tendency with zero
407	precl(:pcols) = 0._r8 ! initialize precipitation rate with zero
408
409	!===============================================================================
410	!
411	! Large-Scale Condensation and Precipitation Rate
412	!
413	!===============================================================================
414	!
415	! Calculate Tendencies
416	!
417	do k=1,pver
418	do i=1,pcols
419	qsat = epsiloe0/pmid(i,k)exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0)) ! saturation specific humidity
420	! out_i(i,llm+1-k)=q(i,k)-qsat
421	if (q(i,k) > qsat) then ! saturated?
422	tmp = 1._r8/dtime(q(i,k)-qsat)/(1._r8+(latvap/cpair)(epsilolatvapqsat/(rairt(i,k)*2)))
423	dtdt(i,k) = dtdt(i,k)+latvap/cpair*tmp
424	dqdt(i,k) = dqdt(i,k)-tmp
425	precl(i) = precl(i) + tmppdel(i,k)/(gravitrhow) ! precipitation rate, computed via a vertical integral
426	! corrected in version 1.3
427	end if
428	end do
429	end do
430	!
431	! Update moisture and temperature fields from Large-Scale Precipitation Scheme
432	!
433	do k=1,pver
434	do i=1,pcols
435	t(i,k) = t(i,k) + dtdt(i,k)*dtime ! update the state variables T and q
436	q(i,k) = q(i,k) + dqdt(i,k)*dtime
437	end do
438	end do
439
440	!===============================================================================
441	! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC
442	!
443	! note: The variables, as done in many GCMs, are written to the history file
444	! after the moist physics process. This ensures that the moisture fields
445	! are somewhat in equilibrium.
446	!===============================================================================
447	! call diag_phys_writeout(state) ! This is for CESM/CAM
448
449	!===============================================================================
450	!
451	! Turbulent mixing coefficients for the PBL mixing of horizontal momentum,
452	! sensible heat and latent heat
453	!
454	! We are using Simplified Ekman theory to compute the diffusion coefficients
455	! Kx for the boundary-layer mixing. The Kx values are calculated at each time step
456	! and in each column.
457	!
458	!===============================================================================
459	!
460	! Compute magnitude of the wind and drag coeffcients for turbulence scheme:
461	! they depend on the conditions at the lowest model level and stay constant
462	! up to the 850 hPa level. Above this level the coefficients are decreased
463	! and tapered to zero. At the 700 hPa level the strength of the K coefficients
464	! is about 10% of the maximum strength.
465	!
466	do i=1,pcols
467	wind(i) = sqrt(u(i,pver)2+v(i,pver)2) ! wind magnitude at the lowest level
468	end do
469	do i=1,pcols
470	Ke(i,pver+1) = Cwind(i)za(i)
471	if( wind(i) .lt. v20) then
472	Cd(i) = Cd0+Cd1*wind(i)
473	Km(i,pver+1) = Cd(i)wind(i)za(i)
474	else
475	Cd(i) = Cm
476	Km(i,pver+1) = Cmwind(i)za(i)
477	endif
478	end do
479
480	do k=1,pver
481	do i=1,pcols
482	if( pint(i,k) .ge. pbltop) then
483	Km(i,k) = Km(i,pver+1) ! constant Km below 850 hPa level
484	Ke(i,k) = Ke(i,pver+1) ! constant Ke below 850 hPa level
485	else
486	Km(i,k) = Km(i,pver+1)exp(-(pbltop-pint(i,k))2/(pblconst)*2) ! Km tapered to 0
487	Ke(i,k) = Ke(i,pver+1)exp(-(pbltop-pint(i,k))2/(pblconst)*2) ! Ke tapered to 0
488	end if
489	end do
490	end do
491
492
493	!===============================================================================
494	! Update the state variables u, v, t, q with the surface fluxes at the
495	! lowest model level, this is done with an implicit approach
496	! see Reed and Jablonowski (JAMES, 2012)
497	!
498	! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1
499	! Tsurf needs to be dependent on latitude for the
500	! moist baroclinic wave test 4-3 with simple-physics, adjust
501	!===============================================================================
502
503	do i=1,pcols
504	qsats = epsiloe0/ps(i)exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0)) ! saturation specific humidity at the surface
505	dudt(i,pver) = dudt(i,pver) + (u(i,pver) &
506	/(1._r8+Cd(i)wind(i)dtime/za(i))-u(i,pver))/dtime
507	dvdt(i,pver) = dvdt(i,pver) + (v(i,pver) &
508	/(1._r8+Cd(i)wind(i)dtime/za(i))-v(i,pver))/dtime
509	u(i,pver) = u(i,pver)/(1._r8+Cd(i)wind(i)dtime/za(i))
510	v(i,pver) = v(i,pver)/(1._r8+Cd(i)wind(i)dtime/za(i))
511	dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+Cwind(i)Tsurf(i)*dtime/za(i)) &
512	/(1._r8+Cwind(i)dtime/za(i))-t(i,pver))/dtime
513	t(i,pver) = (t(i,pver)+Cwind(i)Tsurf(i)*dtime/za(i)) &
514	/(1._r8+Cwind(i)dtime/za(i))
515	dqdt(i,pver) = dqdt(i,pver) +((q(i,pver)+Cwind(i)qsats*dtime/za(i)) &
516	/(1._r8+Cwind(i)dtime/za(i))-q(i,pver))/dtime
517	q(i,pver) = (q(i,pver)+Cwind(i)qsatsdtime/za(i))/(1._r8+Cwind(i)*dtime/za(i))
518	end do
519	!===============================================================================
520
521
522	!===============================================================================
523	! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012)
524	!===============================================================================
525	! Calculate Diagonal Variables for Implicit PBL Scheme
526	!
527	do k=1,pver-1
528	do i=1,pcols
529	rho = (pint(i,k+1)/(rair*(t(i,k+1)+t(i,k))/2.0_r8))
530	CAm(i,k) = rpdel(i,k)dtimegravitgravitKm(i,k+1)rhorho &
531	/(pmid(i,k+1)-pmid(i,k))
532	CCm(i,k+1) = rpdel(i,k+1)dtimegravitgravitKm(i,k+1)rhorho &
533	/(pmid(i,k+1)-pmid(i,k))
534	CA(i,k) = rpdel(i,k)dtimegravitgravitKe(i,k+1)rhorho &
535	/(pmid(i,k+1)-pmid(i,k))
536	CC(i,k+1) = rpdel(i,k+1)dtimegravitgravitKe(i,k+1)rhorho &
537	/(pmid(i,k+1)-pmid(i,k))
538	end do
539	end do
540	do i=1,pcols
541	CAm(i,pver) = 0._r8
542	CCm(i,1) = 0._r8
543	CEm(i,pver+1) = 0._r8
544	CA(i,pver) = 0._r8
545	CC(i,1) = 0._r8
546	CE(i,pver+1) = 0._r8
547	CFu(i,pver+1) = 0._r8
548	CFv(i,pver+1) = 0._r8
549	CFt(i,pver+1) = 0._r8
550	CFq(i,pver+1) = 0._r8
551	end do
552	do i=1,pcols
553	do k=pver,1,-1
554	CE(i,k) = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
555	CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
556	CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1)) &
557	/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
558	CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1)) &
559	/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
560	CFt(i,k) = ((p0/pmid(i,k))*(rair/cpair)t(i,k)+CA(i,k)*CFt(i,k+1)) &
561	/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
562	CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1)) &
563	/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
564	end do
565	end do
566
567	!
568	! Calculate the updated temperature, specific humidity and horizontal wind
569	!
570	! First we need to calculate the updates at the top model level
571	!
572	do i=1,pcols
573	dudt(i,1) = dudt(i,1)+(CFu(i,1)-u(i,1))/dtime
574	dvdt(i,1) = dvdt(i,1)+(CFv(i,1)-v(i,1))/dtime
575	u(i,1) = CFu(i,1)
576	v(i,1) = CFv(i,1)
577	dtdt(i,1) = dtdt(i,1)+(CFt(i,1)(pmid(i,1)/p0)*(rair/cpair)-t(i,1))/dtime ! corrected in version 1.3
578	t(i,1) = CFt(i,1)(pmid(i,1)/p0)*(rair/cpair)
579	dqdt(i,1) = dqdt(i,1)+(CFq(i,1)-q(i,1))/dtime
580	q(i,1) = CFq(i,1)
581	end do
582	!
583	! Loop over the remaining level
584	!
585	do i=1,pcols
586	do k=2,pver
587	dudt(i,k) = dudt(i,k)+(CEm(i,k)*u(i,k-1)+CFu(i,k)-u(i,k))/dtime
588	dvdt(i,k) = dvdt(i,k)+(CEm(i,k)*v(i,k-1)+CFv(i,k)-v(i,k))/dtime
589	u(i,k) = CEm(i,k)*u(i,k-1)+CFu(i,k)
590	v(i,k) = CEm(i,k)*v(i,k-1)+CFv(i,k)
591	dtdt(i,k) = dtdt(i,k)+((CE(i,k)*t(i,k-1) &
592	(p0/pmid(i,k-1))*(rair/cpair)+CFt(i,k)) &
593	(pmid(i,k)/p0)*(rair/cpair)-t(i,k))/dtime
594	t(i,k) = (CE(i,k)t(i,k-1)(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) &
595	(pmid(i,k)/p0)*(rair/cpair)
596	dqdt(i,k) = dqdt(i,k)+(CE(i,k)*q(i,k-1)+CFq(i,k)-q(i,k))/dtime
597	q(i,k) = CE(i,k)*q(i,k-1)+CFq(i,k)
598	end do
599	end do
600
601	!===============================================================================
602	!
603	! Dry Mass Adjustment - THIS PROCESS WILL BE MODEL SPECIFIC
604	!
605	! note: Care needs to be taken to ensure that the model conserves the total
606	! dry air mass. Add your own routine here.
607	!===============================================================================
608	! call physics_dme_adjust(state, tend, qini, dtime) ! This is for CESM/CAM
609
610	return
611	end subroutine simple_physics
612

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

Download in other formats: