Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

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

Context Navigation

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

Download in other formats: