Context Navigation

source: codes/icosagcm/trunk/src/dcmip/physics_dcmip.f90 @ 705

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

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

Download in other formats: