Context Navigation

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

Last change on this file since 192 was 192, checked in by dubos, 10 years ago
Fixed DCMIP5 physics/etat0
File size: 29.5 KB

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

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

Download in other formats: