Context Navigation

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

Last change on this file since 187 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM

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

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

Download in other formats: