Context Navigation

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

Last change on this file since 81 was 81, checked in by ymipsl, 12 years ago

Add simplify physics for dcmip testcase

YM

File size: 27.0 KB

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

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

Download in other formats: