Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_kessler.F @ 2756

Last change on this file since 2756 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 8.9 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3
4	MODULE module_mp_kessler
5
6	CONTAINS
7	!----------------------------------------------------------------
8	SUBROUTINE kessler( t, qv, qc, qr, rho, pii &
9	,dt_in, z, xlv, cp &
10	,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater &
11	,dz8w &
12	,RAINNC, RAINNCV &
13	,ids,ide, jds,jde, kds,kde & ! domain dims
14	,ims,ime, jms,jme, kms,kme & ! memory dims
15	,its,ite, jts,jte, kts,kte & ! tile dims
16	)
17	!----------------------------------------------------------------
18	IMPLICIT NONE
19	!----------------------------------------------------------------
20	! taken from the COMMAS code - WCS 10 May 1999.
21	! converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
22	!----------------------------------------------------------------
23	REAL , PARAMETER :: c1 = .001
24	REAL , PARAMETER :: c2 = .001
25	REAL , PARAMETER :: c3 = 2.2
26	REAL , PARAMETER :: c4 = .875
27	REAL , PARAMETER :: fudge = 1.0
28	REAL , PARAMETER :: mxfall = 10.0
29	!----------------------------------------------------------------
30	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
31	ims,ime, jms,jme, kms,kme, &
32	its,ite, jts,jte, kts,kte
33	REAL , INTENT(IN ) :: xlv, cp
34	REAL , INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
35	REAL , INTENT(IN ) :: rhowater
36
37	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
38	INTENT(INOUT) :: &
39	t , &
40	qv, &
41	qc, &
42	qr
43
44	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
45	INTENT(IN ) :: &
46	rho, &
47	pii, &
48	dz8w
49
50	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
51	INTENT(IN ) :: z
52
53	REAL, INTENT(IN ) :: dt_in
54
55	REAL, DIMENSION( ims:ime , jms:jme ), &
56	INTENT(INOUT) :: RAINNC, &
57	RAINNCV
58
59
60
61	! local variables
62
63	REAL :: qrprod, ern, gam, rcgs, rcgsi
64	REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: prod
65	REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
66	INTEGER :: i,j,k
67	INTEGER :: nfall, n, nfall_new
68	REAL :: qrr, pressure, temp, es, qvs, dz, dt
69	REAL :: f5, dtfall, rdz, product
70	REAL :: max_heating, max_condense, max_rain, maxqrp
71	REAL :: vtmax, ernmax, crmax, factorn, time_sediment
72	REAL :: qcr, factorr, ppt
73	REAL, PARAMETER :: max_cr_sedimentation = 0.75
74	!----------------------------------------------------------------
75
76	INTEGER :: imax, kmax
77
78	dt = dt_in
79
80	! f5 = 237.3 * 17.27 * 2.5e6 / cp
81	f5 = svp2(svpt0-svp3)xlv/cp
82	ernmax = 0.
83	maxqrp = -100.
84
85	!------------------------------------------------------------------------------
86	! parameters for the time split terminal advection
87	!------------------------------------------------------------------------------
88
89	max_heating = 0.
90	max_condense = 0.
91	max_rain = 0.
92
93	!-----------------------------------------------------------------------------
94	! outer J loop for entire microphysics, outer i loop for sedimentation
95	!-----------------------------------------------------------------------------
96
97	microphysics_outer_j_loop: DO j = jts, jte
98
99	sedimentation_outer_i_loop: DO i = its,ite
100
101	! vtmax = 0.
102	crmax = 0.
103
104
105	!------------------------------------------------------------------------------
106	! Terminal velocity calculation and advection, set up coefficients and
107	! compute stable timestep
108	!------------------------------------------------------------------------------
109
110	DO k = 1, kte
111	prodk(k) = qr(i,k,j)
112	rhok(k) = rho(i,k,j)
113	qrr = amax1(0.,qr(i,k,j)0.001rhok(k))
114	vtden(k) = sqrt(rhok(1)/rhok(k))
115	vt(k) = 36.34(qrr0.1364) vtden(k)
116	! vtmax = amax1(vt(k), vtmax)
117	rdzw(k) = 1./dz8w(i,k,j)
118	crmax = amax1(vt(k)dtrdzw(k),crmax)
119	ENDDO
120	DO k = 1, kte-1
121	rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
122	ENDDO
123	rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))
124
125	nfall = max(1,nint(0.5+crmax/max_cr_sedimentation)) ! courant number for big timestep.
126	dtfall = dt / float(nfall) ! splitting so courant number for sedimentation
127	time_sediment = dt ! is stable
128
129	!------------------------------------------------------------------------------
130	! Terminal velocity calculation and advection
131	! Do a time split loop on this for stability.
132	!------------------------------------------------------------------------------
133
134	column_sedimentation: DO WHILE ( nfall > 0 )
135
136	time_sediment = time_sediment - dtfall
137	DO k = 1, kte-1
138	factor(k) = dtfall*rdzk(k)/rhok(k)
139	ENDDO
140	factor(kte) = dtfall*rdzk(kte)
141
142	ppt=0.
143
144	k = 1
145	ppt=rhok(k)prodk(k)vt(k)*dtfall/rhowater
146	RAINNCV(i,j)=ppt*1000.
147	RAINNC(i,j)=RAINNC(i,j)+ppt*1000. ! unit = mm
148
149	!------------------------------------------------------------------------------
150	! Time split loop, Fallout done with flux upstream
151	!------------------------------------------------------------------------------
152
153	DO k = kts, kte-1
154	prodk(k) = prodk(k) - factor(k) &
155	* (rhok(k)prodk(k)vt(k) &
156	-rhok(k+1)prodk(k+1)vt(k+1))
157	ENDDO
158
159	k = kte
160	prodk(k) = prodk(k) - factor(k)prodk(k)vt(k)
161
162	!------------------------------------------------------------------------------
163	! compute new sedimentation velocity, and check/recompute new
164	! sedimentation timestep if this isn't the last split step.
165	!------------------------------------------------------------------------------
166
167	IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep
168
169	nfall = nfall - 1
170	crmax = 0.
171	DO k = kts, kte
172	qrr = amax1(0.,prodk(k)0.001rhok(k))
173	vt(k) = 36.34(qrr0.1364) vtden(k)
174	! vtmax = amax1(vt(k), vtmax)
175	crmax = amax1(vt(k)time_sedimentrdzw(k),crmax)
176	ENDDO
177
178	nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
179	if (nfall_new /= nfall ) then
180	nfall = nfall_new
181	dtfall = time_sediment/nfall
182	end if
183
184	ELSE ! this was the last timestep
185
186	DO k=kts,kte
187	prod(i,k,j) = prodk(k)
188	ENDDO
189	nfall = 0 ! exit condition for sedimentation loop
190
191	END IF
192
193	ENDDO column_sedimentation
194
195	ENDDO sedimentation_outer_i_loop
196
197	!------------------------------------------------------------------------------
198	! Production of rain and deletion of qc
199	! Production of qc from supersaturation
200	! Evaporation of QR
201	!------------------------------------------------------------------------------
202
203	DO k = kts, kte
204	DO i = its, ite
205	factorn = 1.0 / (1.+c3dtamax1(0.,qr(i,k,j))**c4)
206	qrprod = qc(i,k,j) * (1.0 - factorn) &
207	+ factornc1dt*amax1(qc(i,k,j)-c2,0.)
208	rcgs = 0.001*rho(i,k,j)
209
210	qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
211	qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
212	qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)
213
214	temp = pii(i,k,j)*t(i,k,j)
215	pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
216	gam = 2.5e+06/(1004.*pii(i,k,j))
217	! qvs = 380.exp(17.27(temp-273.)/(temp- 36.))/pressure
218	es = 1000.svp1exp(svp2*(temp-svpt0)/(temp-svp3))
219	qvs = ep2*es/(pressure-es)
220	! prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvsf5/(temp-36.)*2)
221	prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)qvsf5/(temp-svp3)**2)
222	ern = amin1(dt(((1.6+124.9(rcgsqr(i,k,j))*.2046) &
223	(rcgsqr(i,k,j))*.525)/(2.55e8/(pressureqvs) &
224	+5.4e5))(dim(qvs,qv(i,k,j))/(rcgsqvs)), &
225	amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))
226
227	! Update all variables
228
229	product = amax1(prod(i,k,j),-qc(i,k,j))
230	t (i,k,j) = t(i,k,j) + gam*(product - ern)
231	qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
232	qc(i,k,j) = qc(i,k,j) + product
233	qr(i,k,j) = qr(i,k,j) - ern
234
235	ENDDO
236	ENDDO
237
238	ENDDO microphysics_outer_j_loop
239
240	RETURN
241
242	END SUBROUTINE kessler
243
244	END MODULE module_mp_kessler

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

Download in other formats:

Original Format