Context Navigation

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

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

Bug fix for DCMIP simple physics

YM

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

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

Download in other formats: