Context Navigation

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

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

implement virtual temperature for moist physics from dcmip

YM

File size: 27.5 KB

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

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

Download in other formats: