Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_gsfcgce.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 96.5 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3
4	MODULE module_mp_gsfcgce
5
6	USE module_wrf_error
7
8	!JJS 1/3/2008 vvvvv
9
10	! common /bt/
11	REAL, PRIVATE :: rd1, rd2, al, cp
12
13	! common /cont/
14	REAL, PRIVATE :: c38, c358, c610, c149, &
15	c879, c172, c409, c76, &
16	c218, c580, c141
17	! common /b3cs/
18	REAL, PRIVATE :: ag, bg, as, bs, &
19	aw, bw, bgh, bgq, &
20	bsh, bsq, bwh, bwq
21
22	! common /size/
23	REAL, PRIVATE :: tnw, tns, tng, &
24	roqs, roqg, roqr
25
26	! common /bterv/
27	REAL, PRIVATE :: zrc, zgc, zsc, &
28	vrc, vgc, vsc
29
30	! common /bsnw/
31	REAL, PRIVATE :: alv, alf, als, t0, t00, &
32	avc, afc, asc, rn1, bnd1, &
33	rn2, bnd2, rn3, rn4, rn5, &
34	rn6, rn7, rn8, rn9, rn10, &
35	rn101,rn10a, rn11,rn11a, rn12
36
37	REAL, PRIVATE :: rn14, rn15,rn15a, rn16, rn17, &
38	rn17a,rn17b,rn17c, rn18, rn18a, &
39	rn19,rn19a,rn19b, rn20, rn20a, &
40	rn20b, bnd3, rn21, rn22, rn23, &
41	rn23a,rn23b, rn25,rn30a, rn30b, &
42	rn30c, rn31, beta, rn32
43
44	REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a
45
46	! common /rsnw1/
47	REAL, PRIVATE :: rn10b, rn10c, rnn191, rnn192, rn30, &
48	rnn30a, rn33, rn331, rn332
49
50	!
51	REAL, PRIVATE, DIMENSION( 31 ) :: aa1, aa2
52	DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, &
53	.3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, &
54	.1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, &
55	.1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, &
56	.4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, &
57	.7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, &
58	.4834e-6/
59	DATA aa2/.4006, .4831, .5320, .5307, .5319, &
60	.5249, .4888, .3894, .4047, .4318, &
61	.4771, .5183, .5463, .5651, .5813, &
62	.5655, .5478, .5203, .4906, .4447, &
63	.4126, .3960, .4149, .4320, .4506, &
64	.4483, .4460, .4433, .4413, .4382, &
65	.4361/
66
67	!JJS 1/3/2008 ^^^^^
68
69	CONTAINS
70
71	!-------------------------------------------------------------------
72	! NASA/GSFC GCE
73	! Tao et al, 2001, Meteo. & Atmos. Phy., 97-137
74	!-------------------------------------------------------------------
75	! SUBROUTINE gsfcgce( th, th_old, &
76	SUBROUTINE gsfcgce( th, &
77	qv, ql, &
78	qr, qi, &
79	qs, &
80	! qvold, qlold, &
81	! qrold, qiold, &
82	! qsold, &
83	rho, pii, p, dt_in, z, &
84	ht, dz8w, grav, &
85	rhowater, rhosnow, &
86	itimestep, &
87	ids,ide, jds,jde, kds,kde, & ! domain dims
88	ims,ime, jms,jme, kms,kme, & ! memory dims
89	its,ite, jts,jte, kts,kte, & ! tile dims
90	rainnc, rainncv, &
91	snownc, snowncv, sr, &
92	graupelnc, graupelncv, &
93	! f_qg, qg, pgold, &
94	f_qg, qg, &
95	ihail, ice2 &
96	)
97
98	!-------------------------------------------------------------------
99	IMPLICIT NONE
100	!-------------------------------------------------------------------
101	!
102	! JJS 2/15/2005
103	!
104	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
105	ims,ime, jms,jme, kms,kme , &
106	its,ite, jts,jte, kts,kte
107	INTEGER, INTENT(IN ) :: itimestep, ihail, ice2
108
109	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
110	INTENT(INOUT) :: &
111	th, &
112	qv, &
113	ql, &
114	qr, &
115	qi, &
116	qs, &
117	qg
118	!
119	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
120	INTENT(IN ) :: &
121	! th_old, &
122	! qvold, &
123	! qlold, &
124	! qrold, &
125	! qiold, &
126	! qsold, &
127	! qgold, &
128	rho, &
129	pii, &
130	p, &
131	dz8w, &
132	z
133
134	REAL, DIMENSION( ims:ime , jms:jme ), &
135	INTENT(INOUT) :: rainnc, &
136	rainncv, &
137	snownc, &
138	snowncv, &
139	sr, &
140	graupelnc, &
141	graupelncv
142
143	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
144
145	REAL, INTENT(IN ) :: dt_in, &
146	grav, &
147	rhowater, &
148	rhosnow
149
150	LOGICAL, INTENT(IN), OPTIONAL :: F_QG
151
152	! LOCAL VAR
153
154
155	!jjs INTEGER :: min_q, max_q
156
157	!jjs REAL, DIMENSION( its:ite , jts:jte ) &
158	!jjs :: rain, snow, graupel,ice
159
160	!
161	! INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id
162	INTEGER :: itaobraun, istatmin, new_ice_sat, id
163
164	INTEGER :: i, j, k
165	INTEGER :: iskip, ih, icount, ibud, i24h
166	REAL :: hour
167	REAL , PARAMETER :: cmin=1.e-20
168	REAL :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2
169	! REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: &
170	! th1, qv1, ql1, qr1, qi1, qs1, qg1
171
172	LOGICAL :: flag_qg
173
174	!
175	!c ihail = 0 for graupel, for tropical region
176	!c ihail = 1 for hail, for mid-lat region
177
178	! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
179	!c if ( itaobraun.eq.1 ) --> betah=0.5beta=-.460.5=-0.23; cn0=1.e-6
180	!c if ( itaobraun.eq.0 ) --> betah=0.5beta=-.60.5=-0.30; cn0=1.e-8
181	itaobraun = 1
182
183	!c ice2 = 0 for 3 ice --- ice, snow and graupel/hail
184	!c ice2 = 1 for 2 ice --- ice and snow only
185	!c ice2 = 2 for 2 ice --- ice and graupel only, use ihail = 0 only
186	!c ice2 = 3 for 0 ice --- no ice, warm only
187
188	! if (ice2 .eq. 2) ihail = 0
189
190	i24h=nint(86400./dt_in)
191	if (mod(itimestep,i24h).eq.1) then
192	write(6,*) 'ihail=',ihail,' ice2=',ice2
193	if (ice2.eq.0) then
194	write(6,*) 'Running 3-ice scheme in GSFCGCE with'
195	if (ihail.eq.0) then
196	write(6,*) ' ice, snow and graupel'
197	else if (ihail.eq.1) then
198	write(6,*) ' ice, snow and hail'
199	else
200	write(6,*) 'ihail has to be either 1 or 0'
201	stop
202	endif !ihail
203	else if (ice2.eq.1) then
204	write(6,*) 'Running 2-ice scheme in GSFCGCE with'
205	write(6,*) ' ice and snow'
206	else if (ice2.eq.2) then
207	write(6,*) 'Running 2-ice scheme in GSFCGCE with'
208	write(6,*) ' ice and graupel'
209	else if (ice2.eq.3) then
210	write(6,*) 'Running warm rain only scheme in GSFCGCE without any ice'
211	else
212	write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, 2, or 3'
213	stop
214	endif !ice2
215	endif !itimestep
216
217	!c new_ice_sat = 0, 1 or 2
218	new_ice_sat = 2
219
220	!c istatmin
221	istatmin = 180
222
223	!c id = 0 without in-line staticstics
224	!c id = 1 with in-line staticstics
225	id = 0
226
227	!c ibud = 0 no calculation of dth, dqv, dqrest and dqall
228	!c ibud = 1 yes
229	ibud = 0
230
231	!jjs dt=dt_in
232	!jjs rhoe_s=1.29
233	!
234	! IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN
235	! CALL wrf_error_fatal3 ( "module_mp_lin.b" , 130 , 'module_mp_lin: Improper use of Lin et al scheme; no ice phase. Please chose another one.')
236	! ENDIF
237
238	! calculte fallflux and precipiation in MKS system
239
240	call fall_flux(dt_in, qr, qi, qs, qg, p, &
241	rho, z, dz8w, ht, rainnc, &
242	rainncv, grav,itimestep, &
243	rhowater, rhosnow, &
244	snownc, snowncv, sr, &
245	graupelnc, graupelncv, &
246	ihail, ice2, &
247	ims,ime, jms,jme, kms,kme, & ! memory dims
248	its,ite, jts,jte, kts,kte ) ! tile dims
249	!-----------------------------------------------------------------------
250
251	!c set up constants used internally in GCE
252
253	call consat_s (ihail, itaobraun)
254
255
256	!c Negative values correction
257
258	iskip = 1
259
260	if (iskip.eq.0) then
261	call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, &
262	itimestep,1, &
263	its,ite,jts,jte,kts,kte)
264	call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, &
265	itimestep,2, &
266	its,ite,jts,jte,kts,kte)
267	call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, &
268	itimestep,3, &
269	its,ite,jts,jte,kts,kte)
270	call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, &
271	itimestep,4, &
272	its,ite,jts,jte,kts,kte)
273	call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, &
274	itimestep,5, &
275	its,ite,jts,jte,kts,kte)
276	call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, &
277	itimestep,6, &
278	its,ite,jts,jte,kts,kte)
279	! else if (mod(itimestep,i24h).eq.1) then
280	! print *,'no neg correction in mp at timestep=',itimestep
281	endif ! iskip
282
283	!c microphysics in GCE
284
285	call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin, &
286	new_ice_sat, id, &
287	! th, th_old, qv, ql, qr, &
288	th, qv, ql, qr, &
289	qi, qs, qg, &
290	! qvold, qlold, qrold, &
291	! qiold, qsold, qgold, &
292	rho, pii, p, itimestep, &
293	ids,ide, jds,jde, kds,kde, & ! domain dims
294	ims,ime, jms,jme, kms,kme, & ! memory dims
295	its,ite, jts,jte, kts,kte & ! tile dims
296	)
297
298	END SUBROUTINE gsfcgce
299
300	!-----------------------------------------------------------------------
301	SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p, &
302	rho, z, dz8w, topo, rainnc, &
303	rainncv, grav, itimestep, &
304	rhowater, rhosnow, &
305	snownc, snowncv, sr, &
306	graupelnc, graupelncv, &
307	ihail, ice2, &
308	ims,ime, jms,jme, kms,kme, & ! memory dims
309	its,ite, jts,jte, kts,kte ) ! tile dims
310	!-----------------------------------------------------------------------
311	! adopted from Jiun-Dar Chern's codes for Purdue Regional Model
312	! adopted by Jainn J. Shi, 6/10/2005
313	!-----------------------------------------------------------------------
314
315	IMPLICIT NONE
316	INTEGER, INTENT(IN ) :: ihail, ice2, &
317	ims,ime, jms,jme, kms,kme, &
318	its,ite, jts,jte, kts,kte
319	INTEGER, INTENT(IN ) :: itimestep
320	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
321	INTENT(INOUT) :: qr, qi, qs, qg
322	REAL, DIMENSION( ims:ime , jms:jme ), &
323	INTENT(INOUT) :: rainnc, rainncv, &
324	snownc, snowncv, sr, &
325	graupelnc, graupelncv
326	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
327	INTENT(IN ) :: rho, z, dz8w, p
328
329	REAL, INTENT(IN ) :: dt, grav, rhowater, rhosnow
330
331
332	REAL, DIMENSION( ims:ime , jms:jme ), &
333	INTENT(IN ) :: topo
334
335	! temperary vars
336
337	REAL, DIMENSION( kts:kte ) :: sqrhoz
338	REAL :: tmp1, term0
339	REAL :: pptrain, pptsnow, &
340	pptgraul, pptice
341	REAL, DIMENSION( kts:kte ) :: qrz, qiz, qsz, qgz, &
342	zz, dzw, prez, rhoz, &
343	orhoz
344
345
346	INTEGER :: k, i, j
347	!
348
349	REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vti
350
351	REAL :: dtb, pi, consta, constc, gambp4, &
352	gamdp4, gam4pt5, gam4bbar
353
354	! Lin
355	REAL , PARAMETER :: xnor = 8.0e6
356	! REAL , PARAMETER :: xnos = 3.0e6
357	REAL , PARAMETER :: xnos = 1.6e7 ! Tao's value
358	REAL , PARAMETER :: &
359	! constb = 0.8, constd = 0.25, o6 = 1./6., &
360	constb = 0.8, constd = 0.11, o6 = 1./6., &
361	cdrag = 0.6
362	! Lin
363	! REAL , PARAMETER :: xnoh = 4.0e4
364	REAL , PARAMETER :: xnoh = 2.0e5 ! Tao's value
365	REAL , PARAMETER :: rhohail = 917.
366
367	! Hobbs
368	REAL , PARAMETER :: xnog = 4.0e6
369	REAL , PARAMETER :: rhograul = 400.
370	REAL , PARAMETER :: abar = 19.3, bbar = 0.37, &
371	p0 = 1.0e5
372
373	REAL , PARAMETER :: rhoe_s = 1.29
374
375	! for terminal velocity flux
376	INTEGER :: min_q, max_q
377	REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
378	LOGICAL :: notlast
379
380	! if (itimestep.eq.1) then
381	! write(6, *) 'in fall_flux'
382	! write(6, *) 'ims=', ims, ' ime=', ime
383	! write(6, *) 'jms=', jms, ' jme=', jme
384	! write(6, *) 'kms=', kms, ' kme=', kme
385	! write(6, *) 'its=', its, ' ite=', ite
386	! write(6, *) 'jts=', jts, ' jte=', jte
387	! write(6, *) 'kts=', kts, ' kte=', kte
388	! write(6, *) 'dt=', dt
389	! write(6, *) 'ihail=', ihail
390	! write(6, *) 'ICE2=', ICE2
391	! write(6, *) 'dt=', dt
392	! endif
393
394	!-----------------------------------------------------------------------
395	! This program calculates precipitation fluxes due to terminal velocities.
396	!-----------------------------------------------------------------------
397
398	dtb=dt
399	pi=acos(-1.)
400	consta=2115.00.01*(1-constb)
401	! constc=152.930.01*(1-constd)
402	constc=78.630.01*(1-constd)
403
404	! Gamma function
405	gambp4=ggamma(constb+4.)
406	gamdp4=ggamma(constd+4.)
407	gam4pt5=ggamma(4.5)
408	gam4bbar=ggamma(4.+bbar)
409
410	!***********************************************************************
411	! Calculate precipitation fluxes due to terminal velocities.
412	!***********************************************************************
413	!
414	!- Calculate termianl velocity (vt?) of precipitation q?z
415	!- Find maximum vt? to determine the small delta t
416
417	j_loop: do j = jts, jte
418	i_loop: do i = its, ite
419
420	pptrain = 0.
421	pptsnow = 0.
422	pptgraul = 0.
423	pptice = 0.
424
425	do k = kts, kte
426	qrz(k)=qr(i,k,j)
427	rhoz(k)=rho(i,k,j)
428	orhoz(k)=1./rhoz(k)
429	prez(k)=p(i,k,j)
430	sqrhoz(k)=sqrt(rhoe_s/rhoz(k))
431	zz(k)=z(i,k,j)
432	dzw(k)=dz8w(i,k,j)
433	enddo !k
434
435	DO k = kts, kte
436	qiz(k)=qi(i,k,j)
437	ENDDO
438
439	DO k = kts, kte
440	qsz(k)=qs(i,k,j)
441	ENDDO
442
443	IF (ice2 .eq. 0) THEN
444	DO k = kts, kte
445	qgz(k)=qg(i,k,j)
446	ENDDO
447	ELSE
448	DO k = kts, kte
449	qgz(k)=0.
450	ENDDO
451	ENDIF
452
453
454	!
455	!-- rain
456	!
457	t_del_tv=0.
458	del_tv=dtb
459	notlast=.true.
460	DO while (notlast)
461	!
462	min_q=kte
463	max_q=kts-1
464	!
465	do k=kts,kte-1
466	if (qrz(k) .gt. 1.0e-8) then
467	min_q=min0(min_q,k)
468	max_q=max0(max_q,k)
469	tmp1=sqrt(pirhowaterxnor/rhoz(k)/qrz(k))
470	tmp1=sqrt(tmp1)
471	vtr(k)=constagambp4sqrhoz(k)/tmp1**constb
472	vtr(k)=vtr(k)/6.
473	if (k .eq. 1) then
474	del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k))
475	else
476	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k))
477	endif
478	else
479	vtr(k)=0.
480	endif
481	enddo
482
483	if (max_q .ge. min_q) then
484	!
485	!- Check if the summation of the small delta t >= big delta t
486	! (t_del_tv) (del_tv) (dtb)
487
488	t_del_tv=t_del_tv+del_tv
489	!
490	if ( t_del_tv .ge. dtb ) then
491	notlast=.false.
492	del_tv=dtb+del_tv-t_del_tv
493	endif
494
495	! use small delta t to calculate the qrz flux
496	! termi is the qrz flux pass in the grid box through the upper boundary
497	! termo is the qrz flux pass out the grid box through the lower boundary
498	!
499	fluxin=0.
500	do k=max_q,min_q,-1
501	fluxout=rhoz(k)vtr(k)qrz(k)
502	flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
503	! tmpqrz=qrz(k)
504	qrz(k)=qrz(k)+del_tv*flux
505	qrz(k)=amax1(0.,qrz(k))
506	qr(i,k,j)=qrz(k)
507	fluxin=fluxout
508	enddo
509	if (min_q .eq. 1) then
510	pptrain=pptrain+fluxin*del_tv
511	else
512	qrz(min_q-1)=qrz(min_q-1)+del_tv* &
513	fluxin/rhoz(min_q-1)/dzw(min_q-1)
514	qr(i,min_q-1,j)=qrz(min_q-1)
515	endif
516	!
517	else
518	notlast=.false.
519	endif
520	ENDDO
521
522	!
523	!-- snow
524	!
525	t_del_tv=0.
526	del_tv=dtb
527	notlast=.true.
528
529	DO while (notlast)
530	!
531	min_q=kte
532	max_q=kts-1
533	!
534	do k=kts,kte-1
535	if (qsz(k) .gt. 1.0e-8) then
536	min_q=min0(min_q,k)
537	max_q=max0(max_q,k)
538	tmp1=sqrt(pirhosnowxnos/rhoz(k)/qsz(k))
539	tmp1=sqrt(tmp1)
540	vts(k)=constcgamdp4sqrhoz(k)/tmp1**constd
541	vts(k)=vts(k)/6.
542	if (k .eq. 1) then
543	del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k))
544	else
545	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k))
546	endif
547	else
548	vts(k)=0.
549	endif
550	enddo
551
552	if (max_q .ge. min_q) then
553	!
554	!
555	!- Check if the summation of the small delta t >= big delta t
556	! (t_del_tv) (del_tv) (dtb)
557
558	t_del_tv=t_del_tv+del_tv
559
560	if ( t_del_tv .ge. dtb ) then
561	notlast=.false.
562	del_tv=dtb+del_tv-t_del_tv
563	endif
564
565	! use small delta t to calculate the qsz flux
566	! termi is the qsz flux pass in the grid box through the upper boundary
567	! termo is the qsz flux pass out the grid box through the lower boundary
568	!
569	fluxin=0.
570	do k=max_q,min_q,-1
571	fluxout=rhoz(k)vts(k)qsz(k)
572	flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
573	qsz(k)=qsz(k)+del_tv*flux
574	qsz(k)=amax1(0.,qsz(k))
575	qs(i,k,j)=qsz(k)
576	fluxin=fluxout
577	enddo
578	if (min_q .eq. 1) then
579	pptsnow=pptsnow+fluxin*del_tv
580	else
581	qsz(min_q-1)=qsz(min_q-1)+del_tv* &
582	fluxin/rhoz(min_q-1)/dzw(min_q-1)
583	qs(i,min_q-1,j)=qsz(min_q-1)
584	endif
585	!
586	else
587	notlast=.false.
588	endif
589
590	ENDDO
591
592	!
593	! ice2=0 --- with hail/graupel
594	! ice2=1 --- without hail/graupel
595	!
596	if (ice2.eq.0) then
597	!
598	!-- If IHAIL=1, use hail.
599	!-- If IHAIL=0, use graupel.
600	!
601	! if (ihail .eq. 1) then
602	! xnog = xnoh
603	! rhograul = rhohail
604	! endif
605
606	t_del_tv=0.
607	del_tv=dtb
608	notlast=.true.
609	!
610	DO while (notlast)
611	!
612	min_q=kte
613	max_q=kts-1
614	!
615	do k=kts,kte-1
616	if (qgz(k) .gt. 1.0e-8) then
617	if (ihail .eq. 1) then
618	! for hail, based on Lin et al (1983)
619	min_q=min0(min_q,k)
620	max_q=max0(max_q,k)
621	tmp1=sqrt(pirhohailxnoh/rhoz(k)/qgz(k))
622	tmp1=sqrt(tmp1)
623	term0=sqrt(4.gravrhohail/3./rhoz(k)/cdrag)
624	vtg(k)=gam4pt5term0sqrt(1./tmp1)
625	vtg(k)=vtg(k)/6.
626	if (k .eq. 1) then
627	del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
628	else
629	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
630	endif !k
631	else
632	! added by JJS
633	! for graupel, based on RH (1984)
634	min_q=min0(min_q,k)
635	max_q=max0(max_q,k)
636	tmp1=sqrt(pirhograulxnog/rhoz(k)/qgz(k))
637	tmp1=sqrt(tmp1)
638	tmp1=tmp1**bbar
639	tmp1=1./tmp1
640	term0=abar*gam4bbar/6.
641	vtg(k)=term0tmp1(p0/prez(k))**0.4
642	if (k .eq. 1) then
643	del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
644	else
645	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
646	endif !k
647	endif !ihail
648	else
649	vtg(k)=0.
650	endif !qgz
651	enddo !k
652
653	if (max_q .ge. min_q) then
654	!
655	!
656	!- Check if the summation of the small delta t >= big delta t
657	! (t_del_tv) (del_tv) (dtb)
658
659	t_del_tv=t_del_tv+del_tv
660
661	if ( t_del_tv .ge. dtb ) then
662	notlast=.false.
663	del_tv=dtb+del_tv-t_del_tv
664	endif
665
666	! use small delta t to calculate the qgz flux
667	! termi is the qgz flux pass in the grid box through the upper boundary
668	! termo is the qgz flux pass out the grid box through the lower boundary
669	!
670	fluxin=0.
671	do k=max_q,min_q,-1
672	fluxout=rhoz(k)vtg(k)qgz(k)
673	flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
674	qgz(k)=qgz(k)+del_tv*flux
675	qgz(k)=amax1(0.,qgz(k))
676	qg(i,k,j)=qgz(k)
677	fluxin=fluxout
678	enddo
679	if (min_q .eq. 1) then
680	pptgraul=pptgraul+fluxin*del_tv
681	else
682	qgz(min_q-1)=qgz(min_q-1)+del_tv* &
683	fluxin/rhoz(min_q-1)/dzw(min_q-1)
684	qg(i,min_q-1,j)=qgz(min_q-1)
685	endif
686	!
687	else
688	notlast=.false.
689	endif
690	!
691	ENDDO
692	ENDIF !ice2
693	!
694	!-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL
695	!
696
697	t_del_tv=0.
698	del_tv=dtb
699	notlast=.true.
700	!
701	DO while (notlast)
702	!
703	min_q=kte
704	max_q=kts-1
705	!
706	do k=kts,kte-1
707	if (qiz(k) .gt. 1.0e-8) then
708	min_q=min0(min_q,k)
709	max_q=max0(max_q,k)
710	vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16 ! Heymsfield and Donner
711	if (k .eq. 1) then
712	del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k))
713	else
714	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k))
715	endif
716	else
717	vti(k)=0.
718	endif
719	enddo
720
721	if (max_q .ge. min_q) then
722	!
723	!
724	!- Check if the summation of the small delta t >= big delta t
725	! (t_del_tv) (del_tv) (dtb)
726
727	t_del_tv=t_del_tv+del_tv
728
729	if ( t_del_tv .ge. dtb ) then
730	notlast=.false.
731	del_tv=dtb+del_tv-t_del_tv
732	endif
733
734	! use small delta t to calculate the qiz flux
735	! termi is the qiz flux pass in the grid box through the upper boundary
736	! termo is the qiz flux pass out the grid box through the lower boundary
737	!
738
739	fluxin=0.
740	do k=max_q,min_q,-1
741	fluxout=rhoz(k)vti(k)qiz(k)
742	flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
743	qiz(k)=qiz(k)+del_tv*flux
744	qiz(k)=amax1(0.,qiz(k))
745	qi(i,k,j)=qiz(k)
746	fluxin=fluxout
747	enddo
748	if (min_q .eq. 1) then
749	pptice=pptice+fluxin*del_tv
750	else
751	qiz(min_q-1)=qiz(min_q-1)+del_tv* &
752	fluxin/rhoz(min_q-1)/dzw(min_q-1)
753	qi(i,min_q-1,j)=qiz(min_q-1)
754	endif
755	!
756	else
757	notlast=.false.
758	endif
759	!
760	ENDDO !notlast
761
762	! prnc(i,j)=prnc(i,j)+pptrain
763	! psnowc(i,j)=psnowc(i,j)+pptsnow
764	! pgrauc(i,j)=pgrauc(i,j)+pptgraul
765	! picec(i,j)=picec(i,j)+pptice
766	!
767
768	! write(6,*) 'i=',i,' j=',j,' ', pptrain, pptsnow, pptgraul, pptice
769	! call flush(6)
770
771	snowncv(i,j) = pptsnow
772	snownc(i,j) = snownc(i,j) + pptsnow
773	graupelncv(i,j) = pptgraul
774	graupelnc(i,j) = graupelnc(i,j) + pptgraul
775	RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
776	RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
777	sr(i,j) = 0.
778	if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j)
779
780	ENDDO i_loop
781	ENDDO j_loop
782
783	! if (itimestep.eq.6480) then
784	! write(51,*) 'in the end of fallflux, itimestep=',itimestep
785	! do j=jts,jte
786	! do i=its,ite
787	! if (rainnc(i,j).gt.400.) then
788	! write(50,*) 'i=',i,' j=',j,' rainnc=',rainnc
789	! endif
790	! enddo
791	! enddo
792	! endif
793
794	END SUBROUTINE fall_flux
795
796	!----------------------------------------------------------------
797	REAL FUNCTION ggamma(X)
798	!----------------------------------------------------------------
799	IMPLICIT NONE
800	!----------------------------------------------------------------
801	REAL, INTENT(IN ) :: x
802	REAL, DIMENSION(8) :: B
803	INTEGER ::j, K1
804	REAL ::PF, G1TO2 ,TEMP
805
806	DATA B/-.577191652,.988205891,-.897056937,.918206857, &
807	-.756704078,.482199394,-.193527818,.035868343/
808
809	PF=1.
810	TEMP=X
811	DO 10 J=1,200
812	IF (TEMP .LE. 2) GO TO 20
813	TEMP=TEMP-1.
814	10 PF=PF*TEMP
815	100 FORMAT(//,5X,'module_gsfcgce: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
816	WRITE(wrf_err_message,100)X
817	CALL wrf_error_fatal(wrf_err_message)
818	20 G1TO2=1.
819	TEMP=TEMP - 1.
820	DO 30 K1=1,8
821	30 G1TO2=G1TO2 + B(K1)TEMP*K1
822	ggamma=PF*G1TO2
823
824	END FUNCTION ggamma
825
826	!-----------------------------------------------------------------------
827	!c Correction of negative values
828	SUBROUTINE negcor ( X, rho, dz8w, &
829	ims,ime, jms,jme, kms,kme, & ! memory dims
830	itimestep, ics, &
831	its,ite, jts,jte, kts,kte ) ! tile dims
832	!-----------------------------------------------------------------------
833	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
834	INTENT(INOUT) :: X
835	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
836	INTENT(IN ) :: rho, dz8w
837	integer, INTENT(IN ) :: itimestep, ics
838
839	!c Local variables
840	! REAL, DIMENSION( kts:kte ) :: Y1, Y2
841	REAL :: A0, A1, A2
842
843	A1=0.
844	A2=0.
845	do k=kts,kte
846	do j=jts,jte
847	do i=its,ite
848	A1=A1+max(X(i,k,j), 0.)rho(i,k,j)dz8w(i,k,j)
849	A2=A2+max(-X(i,k,j), 0.)rho(i,k,j)dz8w(i,k,j)
850	enddo
851	enddo
852	enddo
853
854	! A1=0.0
855	! A2=0.0
856	! do k=kts,kte
857	! A1=A1+Y1(k)
858	! A2=A2+Y2(k)
859	! enddo
860
861	A0=0.0
862
863	if (A1.NE.0.0.and.A1.GT.A2) then
864	A0=(A1-A2)/A1
865
866	if (mod(itimestep,540).eq.0) then
867	if (ics.eq.1) then
868	write(61,*) 'kms=',kms,' kme=',kme,' kts=',kts,' kte=',kte
869	write(61,*) 'jms=',jms,' jme=',jme,' jts=',jts,' jte=',jte
870	write(61,*) 'ims=',ims,' ime=',ime,' its=',its,' ite=',ite
871	endif
872	if (ics.eq.1) then
873	write(61,*) 'qv timestep=',itimestep
874	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
875	else if (ics.eq.2) then
876	write(61,*) 'ql timestep=',itimestep
877	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
878	else if (ics.eq.3) then
879	write(61,*) 'qr timestep=',itimestep
880	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
881	else if (ics.eq.4) then
882	write(61,*) 'qi timestep=',itimestep
883	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
884	else if (ics.eq.5) then
885	write(61,*) 'qs timestep=',itimestep
886	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
887	else if (ics.eq.6) then
888	write(61,*) 'qg timestep=',itimestep
889	write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
890	else
891	write(61,*) 'wrong cloud specieis number'
892	endif
893	endif
894
895	do k=kts,kte
896	do j=jts,jte
897	do i=its,ite
898	X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0)
899	enddo
900	enddo
901	enddo
902	endif
903
904	END SUBROUTINE negcor
905
906	SUBROUTINE consat_s (ihail,itaobraun)
907
908	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
909	! c
910	! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
911	! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
912	! c
913	! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
914	! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
915	! c
916	! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
917	! Model. Part I: Model description. Terrestrial, Atmospheric and c
918	! Oceanic Sciences, 4, 35-72. c
919	! c
920	! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
921	! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
922	! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
923	! radiation and surface processes in the Goddard Cumulus Ensemble c
924	! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
925	! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
926	! c
927	! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
928	! Rutledge, and J. Simpson, 2007: Improving simulations of c
929	! convective system from TRMM LBA: Easterly and Westerly regimes. c
930	! J. Atmos. Sci., 64, 1141-1164. c
931	! c
932	! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
933	! c
934	! Implemented into WRF by Roger Shi 2006/2007 c
935	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
936
937	! itaobraun=0 ! see Tao and Simpson (1993)
938	! itaobraun=1 ! see Tao et al. (2003)
939
940	integer :: itaobraun
941	real :: cn0
942
943	!JJS 1/3/2008 vvvvv
944	!JJS the following common blocks have been moved to the top of
945	!JJS module_mp_gsfcgce_driver_instat.F
946	!
947	! real, dimension (1:31) :: a1, a2
948	! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, &
949	! .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, &
950	! .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, &
951	! .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, &
952	! .6513e-6,.5956e-6,.5333e-6,.4834e-6/
953	! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, &
954	! .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, &
955	! .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, &
956	! .4382,.4361/
957	!JJS 1/3/2008 ^^^^^
958
959
960	! ******************************************************************
961	!JJS
962	al = 2.5e10
963	cp = 1.004e7
964	rd1 = 1.e-3
965	rd2 = 2.2
966	!JJS
967	cpi=4.*atan(1.)
968	cpi2=cpi*cpi
969	grvt=980.
970	cd1=6.e-1
971	cd2=4.grvt/(3.cd1)
972	tca=2.43e3
973	dwv=.226
974	dva=1.718e-4
975	amw=18.016
976	ars=8.314e7
977	scv=2.2904487
978	t0=273.16
979	t00=238.16
980	alv=2.5e10
981	alf=3.336e9
982	als=2.8336e10
983	avc=alv/cp
984	afc=alf/cp
985	asc=als/cp
986	rw=4.615e6
987	cw=4.187e7
988	ci=2.093e7
989	c76=7.66
990	c358=35.86
991	c172=17.26939
992	c409=4098.026
993	c218=21.87456
994	c580=5807.695
995	c610=6.1078e3
996	c149=1.496286e-5
997	c879=8.794142
998	c141=1.4144354e7
999	!*** DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY
1000	!*** DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION
1001	!******** HAIL OR GRAUPEL PARAMETERS ********
1002	if(ihail .eq. 1) then
1003	roqg=.9
1004	ag=sqrt(cd2*roqg)
1005	bg=.5
1006	tng=.002
1007	else
1008	roqg=.4
1009	ag=351.2
1010	! AG=372.3 ! if ice913=1 6/15/02 tao's
1011	bg=.37
1012	tng=.04
1013	endif
1014	!******** SNOW PARAMETERS ********
1015	!ccshie 6/15/02 tao's
1016	! TNS=1.
1017	! TNS=.08 ! if ice913=1, tao's
1018	tns=.16 ! if ice913=0, tao's
1019	roqs=.1
1020	! AS=152.93
1021	as=78.63
1022	! BS=.25
1023	bs=.11
1024	!******** RAIN PARAMETERS ********
1025	aw=2115.
1026	bw=.8
1027	roqr=1.
1028	tnw=.08
1029	!*****************************************************************
1030	bgh=.5*bg
1031	bsh=.5*bs
1032	bwh=.5*bw
1033	bgq=.25*bg
1034	bsq=.25*bs
1035	bwq=.25*bw
1036	!********GAMMA FUNCTION CALCULATIONS***********
1037	ga3b = gammagce(3.+bw)
1038	ga4b = gammagce(4.+bw)
1039	ga6b = gammagce(6.+bw)
1040	ga5bh = gammagce((5.+bw)/2.)
1041	ga3g = gammagce(3.+bg)
1042	ga4g = gammagce(4.+bg)
1043	ga5gh = gammagce((5.+bg)/2.)
1044	ga3d = gammagce(3.+bs)
1045	ga4d = gammagce(4.+bs)
1046	ga5dh = gammagce((5.+bs)/2.)
1047	!CCCCC LIN ET AL., 1983 OR LORD ET AL., 1984 CCCCCCCCCCCCCCCCC
1048	ac1=aw
1049	!JJS
1050	ac2=ag
1051	ac3=as
1052	!JJS
1053	bc1=bw
1054	cc1=as
1055	dc1=bs
1056	zrc=(cpiroqrtnw)**0.25
1057	zsc=(cpiroqstns)**0.25
1058	zgc=(cpiroqgtng)**0.25
1059	vrc=awga4b/(6.zrc**bw)
1060	vsc=asga4d/(6.zsc**bs)
1061	vgc=agga4g/(6.zgc**bg)
1062	! ****************************
1063	! RN1=1.E-3
1064	rn1=9.4e-15 ! 6/15/02 tao's
1065	bnd1=6.e-4
1066	rn2=1.e-3
1067	! BND2=1.25E-3
1068	! BND2=1.5E-3 ! if ice913=1 6/15/02 tao's
1069	bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's
1070	rn3=.25cpitnscc1ga3d
1071	esw=1.
1072	rn4=.25cpieswtnscc1*ga3d
1073	! ERI=1.
1074	eri=.1 ! 6/17/02 tao's ice913=0 (not 1)
1075	rn5=.25cpieritnwac1*ga3b
1076	! AMI=1./(24.*4.19E-10)
1077	ami=1./(24.*6.e-9) ! 6/15/02 tao's
1078	rn6=cpi2eritnwac1roqrga6bami
1079	! ESR=1. ! also if ice913=1 for tao's
1080	esr=.5 ! 6/15/02 for ice913=0 tao's
1081	rn7=cpi2esrtnwtnsroqs
1082	esr=1.
1083	rn8=cpi2esrtnwtnsroqr
1084	rn9=cpi2tnstng*roqs
1085	rn10=2.cpitns
1086	rn101=.31ga5dhsqrt(cc1)
1087	rn10a=als*als/rw
1088	!JJS
1089	rn10b=alv/tca
1090	rn10c=ars/(dwv*amw)
1091	!JJS
1092	rn11=2.cpitns/alf
1093	rn11a=cw/alf
1094	! AMI50=1.51e-7
1095	ami50=3.84e-6 ! 6/15/02 tao's
1096	! AMI40=2.41e-8
1097	ami40=3.08e-8 ! 6/15/02 tao's
1098	eiw=1.
1099	! UI50=20.
1100	ui50=100. ! 6/15/02 tao's
1101	ri50=2.*5.e-3
1102	cmn=1.05e-15
1103	rn12=cpieiwui50ri50*2
1104
1105	do 10 k=1,31
1106	y1=1.-aa2(k)
1107	rn13(k)=aa1(k)y1/(ami50y1-ami40*y1)
1108	rn12a(k)=rn13(k)/ami50
1109	rn12b(k)=aa1(k)ami50*aa2(k)
1110	rn25a(k)=aa1(k)cmn*aa2(k)
1111	10 continue
1112
1113	egw=1.
1114	rn14=.25cpiegwtngga3g*ag
1115	egi=.1
1116	rn15=.25cpiegitngga3g*ag
1117	egi=1.
1118	rn15a=.25cpiegitngga3g*ag
1119	egr=1.
1120	rn16=cpi2egrtngtnwroqr
1121	rn17=2.cpitng
1122	rn17a=.31ga5ghsqrt(ag)
1123	rn17b=cw-ci
1124	rn17c=cw
1125	apri=.66
1126	bpri=1.e-4
1127	bpri=0.5*bpri ! 6/17/02 tao's
1128	rn18=20.cpi2bpritnwroqr
1129	rn18a=apri
1130	rn19=2.cpitng/alf
1131	rn19a=.31ga5ghsqrt(ag)
1132	rn19b=cw/alf
1133	!
1134	rnn191=.78
1135	rnn192=.31ga5ghsqrt(ac2/dva)
1136	!
1137	rn20=2.cpitng
1138	rn20a=als*als/rw
1139	rn20b=.31ga5ghsqrt(ag)
1140	bnd3=2.e-3
1141	rn21=1.e3*1.569e-12/0.15
1142	erw=1.
1143	rn22=.25cpierwac1tnw*ga3b
1144	rn23=2.cpitnw
1145	rn23a=.31ga5bhsqrt(ac1)
1146	rn23b=alv*alv/rw
1147
1148
1149	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1150	!cc
1151	!cc "c0" in routine "consat" (2d), "consatrh" (3d)
1152	!cc if ( itaobraun.eq.1 ) --> betah=0.5beta=-.460.5=-0.23; cn0=1.e-6
1153	!cc if ( itaobraun.eq.0 ) --> betah=0.5beta=-.60.5=-0.30; cn0=1.e-8
1154
1155	if (itaobraun .eq. 0) then
1156	cn0=1.e-8
1157	beta=-.6
1158	elseif (itaobraun .eq. 1) then
1159	cn0=1.e-6
1160	beta=-.46
1161	endif
1162	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1163	! CN0=1.E-6
1164	! CN0=1.E-8 ! 6/15/02 tao's
1165	! BETA=-.46
1166	! BETA=-.6 ! 6/15/02 tao's
1167
1168	rn25=cn0
1169	rn30a=alvalsamw/(tca*ars)
1170	rn30b=alv/tca
1171	rn30c=ars/(dwv*amw)
1172	rn31=1.e-17
1173
1174	rn32=4.*51.545e-4
1175	!
1176	rn30=2.cpitng
1177	rnn30a=alvalvamw/(tca*ars)
1178	!
1179	rn33=4.*tns
1180	rn331=.65
1181	rn332=.44sqrt(ac3/dva)ga5dh
1182	!
1183
1184	return
1185	END SUBROUTINE consat_s
1186
1187	SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, &
1188	new_ice_sat, id, &
1189	ptwrf, qvwrf, qlwrf, qrwrf, &
1190	qiwrf, qswrf, qgwrf, &
1191	rho_mks, pi_mks, p0_mks,itimestep, &
1192	ids,ide, jds,jde, kds,kde, &
1193	ims,ime, jms,jme, kms,kme, &
1194	its,ite, jts,jte, kts,kte &
1195	)
1196	IMPLICIT NONE
1197	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1198	! c
1199	! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
1200	! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
1201	! c
1202	! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
1203	! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
1204	! c
1205	! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
1206	! Model. Part I: Model description. Terrestrial, Atmospheric and c
1207	! Oceanic Sciences, 4, 35-72. c
1208	! c
1209	! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
1210	! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
1211	! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
1212	! radiation and surface processes in the Goddard Cumulus Ensemble c
1213	! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
1214	! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
1215	! c
1216	! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
1217	! Rutledge, and J. Simpson, 2007: Improving simulations of c
1218	! convective system from TRMM LBA: Easterly and Westerly regimes. c
1219	! J. Atmos. Sci., 64, 1141-1164. c
1220	! c
1221	! Tao, W.-K., J. J. Shi, S. Lang, C. Peters-Lidard, A. Hou, S. c
1222	! Braun, and J. Simpson, 2007: New, improved bulk-microphysical c
1223	! schemes for studying precipitation processes in WRF. Part I: c
1224	! Comparisons with other schemes. to appear on Mon. Wea. Rev. C
1225	! c
1226	! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
1227	! c
1228	! Implemented into WRF by Roger Shi 2006/2007 c
1229	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1230	!
1231	! COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES
1232	!
1233	integer, parameter :: nt=2880, nt2=2*nt
1234
1235	!cc using scott braun's way for pint, pidep computations
1236	integer :: itaobraun,ice2,ihail,new_ice_sat,id,istatmin
1237	integer :: itimestep
1238	real :: tairccri, cn0, dt
1239	!cc
1240
1241	!JJS common/bxyz/ n,isec,nran,kt1,kt2
1242	!JJS common/option/ lipps,ijkadv,istatmin,iwater,itoga,imlifting,lin,
1243	!JJS 1 irf,iadvh,irfg,ismg,id
1244
1245	!JJS common/timestat/ ndt_stat,idq
1246	!JJS common/iice/ new_ice_sat
1247	!JJS common/bt/ dt,d2t,rijl2,dts,f5,rd1,rd2,bound,al,cp,ra,ck,ce,eps,
1248	!JJS 1 psfc,fcor,sec,aminut,rdt
1249
1250	!JJS the following common blocks have been moved to the top of
1251	!JJS module_mp_gsfcgce_driver_instat.F
1252
1253	! common/bt/ rd1,rd2,al,cp
1254	!
1255	!
1256	! common/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
1257	! common/size/ tnw,tns,tng,roqs,roqg,roqr
1258	! common/cont/ c38,c358,c610,c149,c879,c172,c409,c76,c218,c580,c141
1259	! common/b3cs/ ag,bg,as,bs,aw,bw,bgh,bgq,bsh,bsq,bwh,bwq
1260	! common/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2, &
1261	! rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a, &
1262	! rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17, &
1263	! rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b, &
1264	! bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b, &
1265	! rn30c,rn31,beta,rn32
1266	! common/rsnw1/ rn10b,rn10c,rnn191,rnn192,rn30,rnn30a,rn33,rn331, &
1267	! rn332
1268	!JJS
1269
1270	integer ids,ide,jds,jde,kds,kde
1271	integer ims,ime,jms,jme,kms,kme
1272	integer its,ite,jts,jte,kts,kte
1273	integer i,j,k, kp
1274
1275	real :: a0 ,a1 ,a2 ,afcp ,alvr ,ami100 ,ami40 ,ami50 ,ascp ,avcp ,betah &
1276	,bg3 ,bgh5 ,bs3 ,bs6 ,bsh5 ,bw3 ,bw6 ,bwh5 ,cmin ,cmin1 ,cmin2 ,cp409 &
1277	,cp580 ,cs580 ,cv409 ,d2t ,del ,dwvp ,ee1 ,ee2 ,f00 ,f2 ,f3 ,ft ,fv0 ,fvs &
1278	,pi0 ,pir ,pr0 ,qb0 ,r00 ,r0s ,r101f ,r10ar ,r10t ,r11at ,r11rt ,r12r ,r14f &
1279	,r14r ,r15af ,r15ar ,r15f ,r15r ,r16r ,r17aq ,r17as ,r17r ,r18r ,r19aq ,r19as &
1280	,r19bt ,r19rt ,r20bq ,r20bs ,r20t ,r22f ,r23af ,r23br ,r23t ,r25a ,r25rt ,r2ice &
1281	,r31r ,r32rt ,r3f ,r4f ,r5f ,r6f ,r7r ,r8r ,r9r ,r_nci ,rft ,rijl2 ,rp0 ,rr0 &
1282	,rrq ,rrs ,rt0 ,scc ,sccc ,sddd ,see ,seee ,sfff ,smmm ,ssss ,tb0 ,temp ,ucog &
1283	,ucor ,ucos ,uwet ,vgcf ,vgcr ,vrcf ,vscf ,zgr ,zrr ,zsr
1284
1285
1286	real, dimension (its:ite,jts:jte,kts:kte) :: fv
1287	real, dimension (its:ite,jts:jte,kts:kte) :: dpt, dqv
1288	real, dimension (its:ite,jts:jte,kts:kte) :: qcl, qrn, &
1289	qci, qcs, qcg
1290	!JJS 10/16/06 vvvv
1291	! real dpt1(ims:ime,jms:jme,kms:kme)
1292	! real dqv1(ims:ime,jms:jme,kms:kme),
1293	! 1 qcl1(ims:ime,jms:jme,kms:kme)
1294	! real qrn1(ims:ime,jms:jme,kms:kme),
1295	! 1 qci1(ims:ime,jms:jme,kms:kme)
1296	! real qcs1(ims:ime,jms:jme,kms:kme),
1297	! 1 qcg1(ims:ime,jms:jme,kms:kme)
1298	!JJS 10/16/06 ^^^^
1299
1300	!JJS
1301
1302	real, dimension (ims:ime, kms:kme, jms:jme) :: ptwrf, qvwrf
1303	real, dimension (ims:ime, kms:kme, jms:jme) :: qlwrf, qrwrf, &
1304	qiwrf, qswrf, qgwrf
1305	!JJS 10/16/06 vvvv
1306	! real ptwrfold(ims:ime, kms:kme, jms:jme)
1307	! real qvwrfold(ims:ime, kms:kme, jms:jme),
1308	! 1 qlwrfold(ims:ime, kms:kme, jms:jme)
1309	! real qrwrfold(ims:ime, kms:kme, jms:jme),
1310	! 1 qiwrfold(ims:ime, kms:kme, jms:jme)
1311	! real qswrfold(ims:ime, kms:kme, jms:jme),
1312	! 1 qgwrfold(ims:ime, kms:kme, jms:jme)
1313	!JJS 10/16/06 ^^^^
1314
1315	!JJS in MKS
1316	real, dimension (ims:ime, kms:kme, jms:jme) :: rho_mks
1317	real, dimension (ims:ime, kms:kme, jms:jme) :: pi_mks
1318	real, dimension (ims:ime, kms:kme, jms:jme) :: p0_mks
1319	!JJS
1320	! real, dimension (its:ite,jts:jte,kts:kte) :: ww1
1321	! real, dimension (its:ite,jts:jte,kts:kte) :: rsw
1322	! real, dimension (its:ite,jts:jte,kts:kte) :: rlw
1323
1324	!JJS COMMON /BADV/
1325	real, dimension (its:ite,jts:jte) :: &
1326	vg, zg, &
1327	ps, pg, &
1328	prn, psn, &
1329	pwacs, wgacr, &
1330	pidep, pint, &
1331	qsi, ssi, &
1332	esi, esw, &
1333	qsw, pr, &
1334	ssw, pihom, &
1335	pidw, pimlt, &
1336	psaut, qracs, &
1337	psaci, psacw, &
1338	qsacw, praci, &
1339	pmlts, pmltg, &
1340	asss, y1, y2
1341	!JJS Y2(its:ite,jts:jte), DDE(NB)
1342
1343	!JJS COMMON/BSAT/
1344	real, dimension (its:ite,jts:jte) :: &
1345	praut, pracw, &
1346	psfw, psfi, &
1347	dgacs, dgacw, &
1348	dgaci, dgacr, &
1349	pgacs, wgacs, &
1350	qgacw, wgaci, &
1351	qgacr, pgwet, &
1352	pgaut, pracs, &
1353	psacr, qsacr, &
1354	pgfr, psmlt, &
1355	pgmlt, psdep, &
1356	pgdep, piacr, &
1357	y5, scv, &
1358	tca, dwv, &
1359	egs, y3, &
1360	y4, ddb
1361
1362	!JJS COMMON/BSAT1/
1363	real, dimension (its:ite,jts:jte) :: &
1364	pt, qv, &
1365	qc, qr, &
1366	qi, qs, &
1367	qg, tair, &
1368	tairc, rtair, &
1369	dep, dd, &
1370	dd1, qvs, &
1371	dm, rq, &
1372	rsub1, col, &
1373	cnd, ern, &
1374	dlt1, dlt2, &
1375	dlt3, dlt4, &
1376	zr, vr, &
1377	zs, vs, &
1378	pssub, &
1379	pgsub, dda
1380
1381	!JJS COMMON/B5/
1382	real, dimension (its:ite,jts:jte,kts:kte) :: rho
1383	real, dimension (kts:kte) :: &
1384	tb, qb, rho1, &
1385	ta, qa, ta1, qa1, &
1386	coef, z1, z2, z3, &
1387	am, am1, ub, vb, &
1388	wb, ub1, vb1, rrho, &
1389	rrho1, wbx
1390
1391	!JJS COMMON/B6/
1392	real, dimension (its:ite,jts:jte,kts:kte) :: p0, pi, f0
1393	real, dimension (kts:kte) :: &
1394	fd, fe, &
1395	st, sv, &
1396	sq, sc, &
1397	se, sqa
1398
1399	!JJS COMMON/BRH1/
1400	real, dimension (kts:kte) :: &
1401	srro, qrro, sqc, sqr, &
1402	sqi, sqs, sqg, stqc, &
1403	stqr, stqi, stqs, stqg
1404	real, dimension (nt) :: &
1405	tqc, tqr, tqi, tqs, tqg
1406
1407	!JJS common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny)
1408
1409	!JJS COMMON/BLS/
1410	real, dimension (ims:ime,jms:jme) :: &
1411	y0, ts0, qss0
1412
1413	!JJS COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4)
1414	integer, dimension (its:ite,jts:jte) :: it
1415	integer, dimension (its:ite,jts:jte, 4) :: ics
1416
1417	integer :: i24h
1418	integer :: iwarm
1419	real :: r2is, r2ig
1420
1421
1422	!JJS COMMON/MICRO/
1423	! real, dimension (ims:ime,kms:kme,jms:jme) :: dbz
1424
1425	!23456789012345678901234567890123456789012345678901234567890123456789012
1426
1427	!
1428	!JJS 1/3/2008 vvvvv
1429	!JJS the following common blocks have been moved to the top of
1430	!JJS module_mp_gsfcgce_driver.F
1431
1432	! real, dimension (31) :: aa1, aa2
1433	! data aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, &
1434	! .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, &
1435	! .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, &
1436	! .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, &
1437	! .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, &
1438	! .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, &
1439	! .4834e-6/
1440	! data aa2/.4006, .4831, .5320, .5307, .5319, &
1441	! .5249, .4888, .3894, .4047, .4318, &
1442	! .4771, .5183, .5463, .5651, .5813, &
1443	! .5655, .5478, .5203, .4906, .4447, &
1444	! .4126, .3960, .4149, .4320, .4506, &
1445	! .4483, .4460, .4433, .4413, .4382, &
1446	! .4361/
1447
1448	!JJS 1/3/2008 ^^^^^
1449
1450	!
1451	!jm 20090220 save
1452
1453	! i24h=nint(86400./dt)
1454	! if (mod(itimestep,i24h).eq.1) then
1455	! write(6, *) 'ims=', ims, ' ime=', ime
1456	! write(6, *) 'jms=', jms, ' jme=', jme
1457	! write(6, *) 'kms=', kms, ' kme=', kme
1458	! write(6, *) 'its=', its, ' ite=', ite
1459	! write(6, *) 'jts=', jts, ' jte=', jte
1460	! write(6, *) 'kts=', kts, ' kte=', kte
1461	! write(6, *) ' ihail=', ihail
1462	! write(6, *) 'itaobraun=',itaobraun
1463	! write(6, *) ' ice2=', ICE2
1464	! write(6, *) 'istatmin=',istatmin
1465	! write(6, *) 'new_ice_sat=', new_ice_sat
1466	! write(6, *) 'id=', id
1467	! write(6, *) 'dt=', dt
1468	! endif
1469
1470	!JJS convert from mks to cgs, and move from WRF grid to GCE grid
1471	do k=kts,kte
1472	do j=jts,jte
1473	do i=its,ite
1474	rho(i,j,k)=rho_mks(i,k,j)*0.001
1475	p0(i,j,k)=p0_mks(i,k,j)*10.0
1476	pi(i,j,k)=pi_mks(i,k,j)
1477	dpt(i,j,k)=ptwrf(i,k,j)
1478	dqv(i,j,k)=qvwrf(i,k,j)
1479	qcl(i,j,k)=qlwrf(i,k,j)
1480	qrn(i,j,k)=qrwrf(i,k,j)
1481	qci(i,j,k)=qiwrf(i,k,j)
1482	qcs(i,j,k)=qswrf(i,k,j)
1483	qcg(i,j,k)=qgwrf(i,k,j)
1484	!JJS 10/16/06 vvvv
1485	! dpt1(i,j,k)=ptwrfold(i,k,j)
1486	! dqv1(i,j,k)=qvwrfold(i,k,j)
1487	! qcl1(i,j,k)=qlwrfold(i,k,j)
1488	! qrn1(i,j,k)=qrwrfold(i,k,j)
1489	! qci1(i,j,k)=qiwrfold(i,k,j)
1490	! qcs1(i,j,k)=qswrfold(i,k,j)
1491	! qcg1(i,j,k)=qgwrfold(i,k,j)
1492	!JJS 10/16/06 ^^^^
1493	enddo !i
1494	enddo !j
1495	enddo !k
1496
1497	do k=kts,kte
1498	do j=jts,jte
1499	do i=its,ite
1500	fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k))
1501	enddo !i
1502	enddo !j
1503	enddo !k
1504	!JJS
1505
1506	!
1507	! **** THREE CLASSES OF ICE-PHASE (LIN ET AL, 1983) *******
1508
1509	!JJS D22T=D2T
1510	!JJS IF(IJKADV .EQ. 0) THEN
1511	!JJS D2T=D2T
1512	!JJS ELSE
1513	d2t=dt
1514	!JJS ENDIF
1515	!
1516	! itaobraun=0 ! original pint and pidep & see Tao and Simpson 1993
1517	itaobraun=1 ! see Tao et al. (2003)
1518	!
1519	if ( itaobraun.eq.0 ) then
1520	cn0=1.e-8
1521	!c beta=-.6
1522	elseif ( itaobraun.eq.1 ) then
1523	cn0=1.e-6
1524	! cn0=1.e-8 ! special
1525	!c beta=-.46
1526	endif
1527	!C TAO 2007 START
1528	! ICE2=0 ! default, 3ice with loud ice, snow and graupel
1529	! r2is=1., r2ig=1.
1530	! ICE2=1 ! 2ice with cloud ice and snow (no graupel) - r2iceg=1, r2ice=0.
1531	! r2is=1., r2ig=0.
1532	! ICE2=2 ! 2ice with cloud ice and graupel (no snow) - r2ice=1, r2iceg=0.
1533	! r2is=0., r2ig=1.
1534	!c
1535	! r2ice=1.
1536	! r2iceg=1.
1537	r2ig=1.
1538	r2is=1.
1539	if (ice2 .eq. 1) then
1540	! r2ice=0.
1541	! r2iceg=1.
1542	r2ig=0.
1543	r2is=1.
1544	endif
1545	if (ice2 .eq. 2) then
1546	! r2ice=1.
1547	! r2iceg=0.
1548	r2ig=1.
1549	r2is=0.
1550	endif
1551	!C TAO 2007 END
1552
1553	!JJS 10/7/2008
1554	! ICE2=3 ! no ice, warm rain only
1555	iwarm = 0
1556	if (ice2 .eq. 3 ) iwarm = 1
1557
1558
1559
1560	cmin=1.e-19
1561	cmin1=1.e-20
1562	cmin2=1.e-12
1563	ucor=3071.29/tnw**0.75
1564	ucos=687.97roqs0.25/tns*0.75
1565	ucog=687.97roqg0.25/tng*0.75
1566	uwet=4.464**0.95
1567
1568	rijl2 = 1. / (ide-ids) / (jde-jds)
1569
1570	!JJScap $doacross local(j,i)
1571
1572	!JJS DO 1 J=1,JMAX
1573	!JJS DO 1 I=1,IMAX
1574	do j=jts,jte
1575	do i=its,ite
1576	it(i,j)=1
1577	enddo
1578	enddo
1579
1580	f2=rd1*d2t
1581	f3=rd2*d2t
1582
1583	ft=dt/d2t
1584	rft=rijl2*ft
1585	a0=.5istatminrijl2
1586	rt0=1./(t0-t00)
1587	bw3=bw+3.
1588	bs3=bs+3.
1589	bg3=bg+3.
1590	bsh5=2.5+bsh
1591	bgh5=2.5+bgh
1592	bwh5=2.5+bwh
1593	bw6=bw+6.
1594	bs6=bs+6.
1595	betah=.5*beta
1596	r10t=rn10*d2t
1597	r11at=rn11a*d2t
1598	r19bt=rn19b*d2t
1599	r20t=-rn20*d2t
1600	r23t=-rn23*d2t
1601	r25a=rn25
1602
1603	! ami50 for use in PINT
1604	ami50=3.76e-8
1605	ami100=1.51e-7
1606	ami40=2.41e-8
1607
1608	!C ******************************************************************
1609
1610	!JJS DO 1000 K=2,kles
1611	do 1000 k=kts,kte
1612	kp=k+1
1613	!JJS tb0=ta1(k)
1614	!JJS qb0=qa1(k)
1615	tb0=0.
1616	qb0=0.
1617
1618	do 2000 j=jts,jte
1619	do 2000 i=its,ite
1620
1621	rp0=3.799052e3/p0(i,j,k)
1622	pi0=pi(i,j,k)
1623	pir=1./(pi(i,j,k))
1624	pr0=1./p0(i,j,k)
1625	r00=rho(i,j,k)
1626	r0s=sqrt(rho(i,j,k))
1627	!JJS RR0=RRHO(K)
1628	rr0=1./rho(i,j,k)
1629	!JJS RRS=SRRO(K)
1630	rrs=sqrt(rr0)
1631	!JJS RRQ=QRRO(K)
1632	rrq=sqrt(rrs)
1633	f0(i,j,k)=al/cp/pi(i,j,k)
1634	f00=f0(i,j,k)
1635	fv0=fv(i,j,k)
1636	fvs=sqrt(fv(i,j,k))
1637	zrr=1.e5zrcrrq
1638	zsr=1.e5zscrrq
1639	zgr=1.e5zgcrrq
1640	cp409=c409*pi0
1641	cv409=c409*avc
1642	cp580=c580*pi0
1643	cs580=c580*asc
1644	alvr=r00*alv
1645	afcp=afc*pir
1646	avcp=avc*pir
1647	ascp=asc*pir
1648	vrcf=vrc*fv0
1649	vscf=vsc*fv0
1650	vgcf=vgc*fv0
1651	vgcr=vgc*rrs
1652	dwvp=c879*pr0
1653	r3f=rn3*fv0
1654	r4f=rn4*fv0
1655	r5f=rn5*fv0
1656	r6f=rn6*fv0
1657	r7r=rn7*rr0
1658	r8r=rn8*rr0
1659	r9r=rn9*rr0
1660	r101f=rn101*fvs
1661	r10ar=rn10a*r00
1662	r11rt=rn11rr0d2t
1663	r12r=rn12*r00
1664	r14r=rn14*rrs
1665	r14f=rn14*fv0
1666	r15r=rn15*rrs
1667	r15ar=rn15a*rrs
1668	r15f=rn15*fv0
1669	r15af=rn15a*fv0
1670	r16r=rn16*rr0
1671	r17r=rn17*rr0
1672	r17aq=rn17a*rrq
1673	r17as=rn17a*fvs
1674	r18r=rn18*rr0
1675	r19rt=rn19rr0d2t
1676	r19aq=rn19a*rrq
1677	r19as=rn19a*fvs
1678	r20bq=rn20b*rrq
1679	r20bs=rn20b*fvs
1680	r22f=rn22*fv0
1681	r23af=rn23a*fvs
1682	r23br=rn23b*r00
1683	r25rt=rn25rr0d2t
1684	r31r=rn31*rr0
1685	r32rt=rn32d2trrs
1686
1687	!JJS DO 100 J=3,JLES
1688	!JJS DO 100 I=3,ILES
1689	pt(i,j)=dpt(i,j,k)
1690	qv(i,j)=dqv(i,j,k)
1691	qc(i,j)=qcl(i,j,k)
1692	qr(i,j)=qrn(i,j,k)
1693	qi(i,j)=qci(i,j,k)
1694	qs(i,j)=qcs(i,j,k)
1695	qg(i,j)=qcg(i,j,k)
1696	! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
1697	if (qc(i,j) .le. cmin1) qc(i,j)=0.0
1698	if (qr(i,j) .le. cmin1) qr(i,j)=0.0
1699	if (qi(i,j) .le. cmin1) qi(i,j)=0.0
1700	if (qs(i,j) .le. cmin1) qs(i,j)=0.0
1701	if (qg(i,j) .le. cmin1) qg(i,j)=0.0
1702	tair(i,j)=(pt(i,j)+tb0)*pi0
1703	tairc(i,j)=tair(i,j)-t0
1704	zr(i,j)=zrr
1705	zs(i,j)=zsr
1706	zg(i,j)=zgr
1707	vr(i,j)=0.0
1708	vs(i,j)=0.0
1709	vg(i,j)=0.0
1710
1711	!JJS 10/7/2008 vvvvv
1712	IF (IWARM .EQ. 1) THEN
1713	!JJS for calculating processes related to warm rain only
1714	qi(i,j)=0.0
1715	qs(i,j)=0.0
1716	qg(i,j)=0.0
1717	dep(i,j)=0.
1718	pint(i,j)=0.
1719	psdep(i,j)=0.
1720	pgdep(i,j)=0.
1721	dd1(i,j)=0.
1722	pgsub(i,j)=0.
1723	psmlt(i,j)=0.
1724	pgmlt(i,j)=0.
1725	pimlt(i,j)=0.
1726	psacw(i,j)=0.
1727	piacr(i,j)=0.
1728	psfw(i,j)=0.
1729	pgfr(i,j)=0.
1730	dgacw(i,j)=0.
1731	dgacr(i,j)=0.
1732	psacr(i,j)=0.
1733	wgacr(i,j)=0.
1734	pihom(i,j)=0.
1735	pidw(i,j)=0.
1736
1737	if (qr(i,j) .gt. cmin1) then
1738	dd(i,j)=r00*qr(i,j)
1739	y1(i,j)=dd(i,j)**.25
1740	zr(i,j)=zrc/y1(i,j)
1741	vr(i,j)=max(vrcfdd(i,j)*bwq, 0.)
1742	endif
1743
1744	!* 21 * PRAUT AUTOCONVERSION OF QC TO QR 21
1745	!* 22 * PRACW : ACCRETION OF QC BY QR 22
1746	pracw(i,j)=0.
1747	praut(i,j)=0.0
1748	pracw(i,j)=r22fqc(i,j)/zr(i,j)*bw3
1749	y1(i,j)=qc(i,j)-bnd3
1750	if (y1(i,j).gt.0.0) then
1751	praut(i,j)=r00y1(i,j)y1(i,j)/(1.2e-4+rn21/y1(i,j))
1752	endif
1753
1754	!C****** HANDLING THE NEGATIVE CLOUD WATER (QC) ****************
1755	Y1(I,J)=QC(I,J)/D2T
1756	PRAUT(I,J)=MIN(Y1(I,J), PRAUT(I,J))
1757	PRACW(I,J)=MIN(Y1(I,J), PRACW(I,J))
1758	Y1(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
1759
1760	if (qc(i,j) .lt. y1(i,j) .and. y1(i,j) .ge. cmin2) then
1761	y2(i,j)=qc(i,j)/(y1(i,j)+cmin2)
1762	praut(i,j)=praut(i,j)*y2(i,j)
1763	pracw(i,j)=pracw(i,j)*y2(i,j)
1764	qc(i,j)=0.0
1765	else
1766	qc(i,j)=qc(i,j)-y1(i,j)
1767	endif
1768
1769	PR(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
1770	QR(I,J)=QR(I,J)+PR(I,J)
1771
1772	!*** TAO ET AL (1989) SATURATION TECHNIQUE *********************
1773
1774	cnd(i,j)=0.0
1775	tair(i,j)=(pt(i,j)+tb0)*pi0
1776	y1(i,j)=1./(tair(i,j)-c358)
1777	qsw(i,j)=rp0exp(c172-c409y1(i,j))
1778	dd(i,j)=cp409y1(i,j)y1(i,j)
1779	dm(i,j)=qv(i,j)+qb0-qsw(i,j)
1780	cnd(i,j)=dm(i,j)/(1.+avcpdd(i,j)qsw(i,j))
1781	!c **** condensation or evaporation of qc ****
1782	cnd(i,j)=max(-qc(i,j), cnd(i,j))
1783	pt(i,j)=pt(i,j)+avcp*cnd(i,j)
1784	qv(i,j)=qv(i,j)-cnd(i,j)
1785	qc(i,j)=qc(i,j)+cnd(i,j)
1786
1787	!C **** EVAPORATION ****
1788	!* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) 23
1789	ern(i,j)=0.0
1790
1791	if(qr(i,j).gt.0.0) then
1792	tair(i,j)=(pt(i,j)+tb0)*pi0
1793	rtair(i,j)=1./(tair(i,j)-c358)
1794	qsw(i,j)=rp0exp(c172-c409rtair(i,j))
1795	ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
1796	dm(i,j)=qv(i,j)+qb0-qsw(i,j)
1797	rsub1(i,j)=cv409qsw(i,j)rtair(i,j)*rtair(i,j)
1798	dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
1799	y1(i,j)=.78/zr(i,j)*2+r23afscv(i,j)/zr(i,j)**bwh5
1800	y2(i,j)=r23br/(tca(i,j)tair(i,j)*2)+1./(dwv(i,j) &
1801	*qsw(i,j))
1802	!cccc
1803	ern(i,j)=r23tssw(i,j)y1(i,j)/y2(i,j)
1804	ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
1805	pt(i,j)=pt(i,j)-avcp*ern(i,j)
1806	qv(i,j)=qv(i,j)+ern(i,j)
1807	qr(i,j)=qr(i,j)-ern(i,j)
1808	endif
1809
1810	ELSE ! part of if (iwarm.eq.1) then
1811	!JJS 10/7/2008 ^^^^^
1812
1813	!JJS for calculating processes related to both ice and warm rain
1814
1815	! * COMPUTE ZR,ZS,ZG,VR,VS,VG ***************************
1816
1817	if (qr(i,j) .gt. cmin1) then
1818	dd(i,j)=r00*qr(i,j)
1819	y1(i,j)=dd(i,j)**.25
1820	zr(i,j)=zrc/y1(i,j)
1821	vr(i,j)=max(vrcfdd(i,j)*bwq, 0.)
1822	endif
1823
1824	if (qs(i,j) .gt. cmin1) then
1825	dd(i,j)=r00*qs(i,j)
1826	y1(i,j)=dd(i,j)**.25
1827	zs(i,j)=zsc/y1(i,j)
1828	vs(i,j)=max(vscfdd(i,j)*bsq, 0.)
1829	endif
1830
1831	if (qg(i,j) .gt. cmin1) then
1832	dd(i,j)=r00*qg(i,j)
1833	y1(i,j)=dd(i,j)**.25
1834	zg(i,j)=zgc/y1(i,j)
1835	if(ihail .eq. 1) then
1836	vg(i,j)=max(vgcrdd(i,j)*bgq, 0.)
1837	else
1838	vg(i,j)=max(vgcfdd(i,j)*bgq, 0.)
1839	endif
1840	endif
1841
1842	if (qr(i,j) .le. cmin2) vr(i,j)=0.0
1843	if (qs(i,j) .le. cmin2) vs(i,j)=0.0
1844	if (qg(i,j) .le. cmin2) vg(i,j)=0.0
1845
1846	! ******************************************************************
1847	! *** Y1 : DYNAMIC VISCOSITY OF AIR (U)
1848	! *** DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
1849	! *** TCA : THERMAL CONDUCTIVITY OF AIR (KA)
1850	! *** Y2 : KINETIC VISCOSITY (V)
1851
1852	y1(i,j)=c149tair(i,j)*1.5/(tair(i,j)+120.)
1853	dwv(i,j)=dwvptair(i,j)*1.81
1854	tca(i,j)=c141*y1(i,j)
1855	scv(i,j)=1./((rr0y1(i,j)).1666667dwv(i,j)**.3333333)
1856	!JJS 100 CONTINUE
1857
1858	!* 1 * PSAUT : AUTOCONVERSION OF QI TO QS *1
1859	!* 3 * PSACI : ACCRETION OF QI TO QS *3
1860	!* 4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT) *4
1861	!* 5 * PRACI : ACCRETION OF QI BY QR *5
1862	!* 6 * PIACR : ACCRETION OF QR OR QG BY QI *6
1863
1864	!JJS DO 125 J=3,JLES
1865	!JJS DO 125 I=3,ILES
1866	psaut(i,j)=0.0
1867	psaci(i,j)=0.0
1868	praci(i,j)=0.0
1869	piacr(i,j)=0.0
1870	psacw(i,j)=0.0
1871	qsacw(i,j)=0.0
1872	dd(i,j)=1./zs(i,j)**bs3
1873
1874	if (tair(i,j).lt.t0) then
1875	esi(i,j)=exp(.025*tairc(i,j))
1876	psaut(i,j)=r2ismax(rn1esi(i,j)*(qi(i,j)-bnd1) ,0.0)
1877	psaci(i,j)=r2isr3fesi(i,j)qi(i,j)dd(i,j)
1878	!JJS 3/30/06
1879	! to cut water to snow accretion by half
1880	! PSACW(I,J)=R4FQC(I,J)DD(I,J)
1881	psacw(i,j)=r2is0.5r4fqc(i,j)dd(i,j)
1882	!JJS 3/30/06
1883	praci(i,j)=r2isr5fqi(i,j)/zr(i,j)**bw3
1884	piacr(i,j)=r2isr6fqi(i,j)(zr(i,j)*(-bw6))
1885	!JJS PIACR(I,J)=R6FQI(I,J)/ZR(I,J)*BW6
1886	else
1887	qsacw(i,j)=r2isr4fqc(i,j)*dd(i,j)
1888	endif
1889
1890	!* 21 * PRAUT AUTOCONVERSION OF QC TO QR 21
1891	!* 22 * PRACW : ACCRETION OF QC BY QR 22
1892
1893	pracw(i,j)=r22fqc(i,j)/zr(i,j)*bw3
1894	praut(i,j)=0.0
1895	y1(i,j)=qc(i,j)-bnd3
1896	if (y1(i,j).gt.0.0) then
1897	praut(i,j)=r00y1(i,j)y1(i,j)/(1.2e-4+rn21/y1(i,j))
1898	endif
1899
1900	!* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971) 12
1901	!* 13 * PSFI : BERGERON PROCESSES FOR QS 13
1902
1903	psfw(i,j)=0.0
1904	psfi(i,j)=0.0
1905	pidep(i,j)=0.0
1906
1907	if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then
1908	y1(i,j)=max( min(tairc(i,j), -1.), -31.)
1909	it(i,j)=int(abs(y1(i,j)))
1910	y1(i,j)=rn12a(it(i,j))
1911	y2(i,j)=rn12b(it(i,j))
1912	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1913	psfw(i,j)=r2is* &
1914	max(d2ty1(i,j)(y2(i,j)+r12rqc(i,j))qi(i,j),0.0)
1915	rtair(i,j)=1./(tair(i,j)-c76)
1916	y2(i,j)=exp(c218-c580*rtair(i,j))
1917	qsi(i,j)=rp0*y2(i,j)
1918	esi(i,j)=c610*y2(i,j)
1919	ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
1920	r_nci=min(1.e-6exp(-.46tairc(i,j)),1.)
1921	! R_NCI=min(1.e-8EXP(-.6TAIRC(I,J)),1.) ! use Tao's
1922	dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
1923	rsub1(i,j)=cs580qsi(i,j)rtair(i,j)*rtair(i,j)
1924	y3(i,j)=1./tair(i,j)
1925	dd(i,j)=y3(i,j)(rn30ay3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j)
1926	y1(i,j)=206.18*ssi(i,j)/dd(i,j)
1927	pidep(i,j)=y1(i,j)sqrt(r_nciqi(i,j)/r00)
1928	dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t
1929	if(dm(i,j).gt.cmin2) then
1930	a2=1.
1931	if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then
1932	a2=dep(i,j)/pidep(i,j)
1933	pidep(i,j)=dep(i,j)
1934	endif
1935	psfi(i,j)=r2isa2.5qi(i,j)y1(i,j)/(sqrt(ami100) &
1936	-sqrt(ami40))
1937	elseif(dm(i,j).lt.-cmin2) then
1938	!
1939	! SUBLIMATION TERMS USED ONLY WHEN SATURATION ADJUSTMENT FOR ICE
1940	! IS TURNED OFF
1941	!
1942	pidep(i,j)=0.
1943	psfi(i,j)=0.
1944	else
1945	pidep(i,j)=0.
1946	psfi(i,j)=0.
1947	endif
1948	endif
1949
1950	!TTT***** QG=QG+MIN(PGDRY,PGWET)
1951	!* 9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET) *9
1952	!* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT) 14
1953	!* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT) 16
1954
1955	if(qc(i,j)+qr(i,j).lt.1.e-4) then
1956	ee1=.01
1957	else
1958	ee1=1.
1959	endif
1960	ee2=0.09
1961	egs(i,j)=ee1exp(ee2tairc(i,j))
1962	! EGS(I,J)=0.1 ! 6/15/02 tao's
1963	if (tair(i,j).ge.t0) egs(i,j)=1.0
1964	y1(i,j)=abs(vg(i,j)-vs(i,j))
1965	y2(i,j)=zs(i,j)*zg(i,j)
1966	y3(i,j)=5./y2(i,j)
1967	y4(i,j)=.08y3(i,j)y3(i,j)
1968	y5(i,j)=.05y3(i,j)y4(i,j)
1969	dd(i,j)=y1(i,j)(y3(i,j)/zs(i,j)5+y4(i,j)/zs(i,j)*3 &
1970	+y5(i,j)/zs(i,j))
1971	pgacs(i,j)=r2igr2isr9regs(i,j)dd(i,j)
1972	!JJS 1/3/06 from Steve and Chunglin
1973	if (ihail.eq.1) then
1974	dgacs(i,j)=pgacs(i,j)
1975	else
1976	dgacs(i,j)=0.
1977	endif
1978	!JJS 1/3/06 from Steve and Chunglin
1979	wgacs(i,j)=r2igr2isr9r*dd(i,j)
1980	! WGACS(I,J)=0. ! 6/15/02 tao's
1981	y1(i,j)=1./zg(i,j)**bg3
1982
1983	if(ihail .eq. 1) then
1984	dgacw(i,j)=r2igmax(r14rqc(i,j)*y1(i,j), 0.0)
1985	else
1986	dgacw(i,j)=r2igmax(r14fqc(i,j)*y1(i,j), 0.0)
1987	endif
1988
1989	qgacw(i,j)=dgacw(i,j)
1990	y1(i,j)=abs(vg(i,j)-vr(i,j))
1991	y2(i,j)=zr(i,j)*zg(i,j)
1992	y3(i,j)=5./y2(i,j)
1993	y4(i,j)=.08y3(i,j)y3(i,j)
1994	y5(i,j)=.05y3(i,j)y4(i,j)
1995	dd(i,j)=r16ry1(i,j)(y3(i,j)/zr(i,j)5+y4(i,j)/zr(i,j)3 &
1996	+y5(i,j)/zr(i,j))
1997	dgacr(i,j)=r2ig*max(dd(i,j), 0.0)
1998	qgacr(i,j)=dgacr(i,j)
1999
2000	if (tair(i,j).ge.t0) then
2001	dgacs(i,j)=0.0
2002	wgacs(i,j)=0.0
2003	dgacw(i,j)=0.0
2004	dgacr(i,j)=0.0
2005	else
2006	pgacs(i,j)=0.0
2007	qgacw(i,j)=0.0
2008	qgacr(i,j)=0.0
2009	endif
2010
2011	!*****PGDRY : DGACW+DGACI+DGACR+DGACS ****
2012	!* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH) 15
2013	!* 17 * PGWET : WET GROWTH OF QG 17
2014
2015	dgaci(i,j)=0.0
2016	wgaci(i,j)=0.0
2017	pgwet(i,j)=0.0
2018
2019	if (tair(i,j).lt.t0) then
2020	y1(i,j)=qi(i,j)/zg(i,j)**bg3
2021	if (ihail.eq.1) then
2022	dgaci(i,j)=r2igr15ry1(i,j)
2023	wgaci(i,j)=r2igr15ary1(i,j)
2024	! WGACI(I,J)=0. ! 6/15/02 tao's
2025	else
2026
2027	!JJS DGACI(I,J)=r2igR15FY1(I,J)
2028	dgaci(i,j)=0.
2029	wgaci(i,j)=r2igr15afy1(i,j)
2030	! WGACI(I,J)=0. ! 6/15/02 tao's
2031	endif
2032	!
2033	if (tairc(i,j).ge.-50.) then
2034	if (alf+rn17c*tairc(i,j) .eq. 0.) then
2035	write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j)
2036	endif
2037	y1(i,j)=1./(alf+rn17c*tairc(i,j))
2038	if (ihail.eq.1) then
2039	y3(i,j)=.78/zg(i,j)*2+r17aqscv(i,j)/zg(i,j)**bgh5
2040	else
2041	y3(i,j)=.78/zg(i,j)*2+r17asscv(i,j)/zg(i,j)**bgh5
2042	endif
2043	y4(i,j)=alvrdwv(i,j)(rp0-(qv(i,j)+qb0)) &
2044	-tca(i,j)*tairc(i,j)
2045	dd(i,j)=y1(i,j)(r17ry4(i,j)*y3(i,j) &
2046	+(wgaci(i,j)+wgacs(i,j))(alf+rn17btairc(i,j)))
2047	pgwet(i,j)=r2ig*max(dd(i,j), 0.0)
2048	endif
2049	endif
2050	!JJS 125 CONTINUE
2051
2052	!****** HANDLING THE NEGATIVE CLOUD WATER (QC) ****************
2053	!****** HANDLING THE NEGATIVE CLOUD ICE (QI) ****************
2054
2055	!JJS DO 150 J=3,JLES
2056	!JJS DO 150 I=3,ILES
2057
2058	y1(i,j)=qc(i,j)/d2t
2059	psacw(i,j)=min(y1(i,j), psacw(i,j))
2060	praut(i,j)=min(y1(i,j), praut(i,j))
2061	pracw(i,j)=min(y1(i,j), pracw(i,j))
2062	psfw(i,j)= min(y1(i,j), psfw(i,j))
2063	dgacw(i,j)=min(y1(i,j), dgacw(i,j))
2064	qsacw(i,j)=min(y1(i,j), qsacw(i,j))
2065	qgacw(i,j)=min(y1(i,j), qgacw(i,j))
2066
2067	y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) &
2068	+dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t
2069	qc(i,j)=qc(i,j)-y1(i,j)
2070
2071	if (qc(i,j) .lt. 0.0) then
2072	a1=1.
2073	if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1.
2074	psacw(i,j)=psacw(i,j)*a1
2075	praut(i,j)=praut(i,j)*a1
2076	pracw(i,j)=pracw(i,j)*a1
2077	psfw(i,j)=psfw(i,j)*a1
2078	dgacw(i,j)=dgacw(i,j)*a1
2079	qsacw(i,j)=qsacw(i,j)*a1
2080	qgacw(i,j)=qgacw(i,j)*a1
2081	qc(i,j)=0.0
2082	endif
2083	!c
2084	!
2085	!******** SHED PROCESS (WGACR=PGWET-DGACW-WGACI-WGACS)
2086	!c
2087	wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j)
2088	y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j)
2089	if (pgwet(i,j).ge.y2(i,j)) then
2090	wgacr(i,j)=0.0
2091	wgaci(i,j)=0.0
2092	wgacs(i,j)=0.0
2093	else
2094	dgacr(i,j)=0.0
2095	dgaci(i,j)=0.0
2096	dgacs(i,j)=0.0
2097	endif
2098	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2099	!c
2100	y1(i,j)=qi(i,j)/d2t
2101	psaut(i,j)=min(y1(i,j), psaut(i,j))
2102	psaci(i,j)=min(y1(i,j), psaci(i,j))
2103	praci(i,j)=min(y1(i,j), praci(i,j))
2104	psfi(i,j)= min(y1(i,j), psfi(i,j))
2105	dgaci(i,j)=min(y1(i,j), dgaci(i,j))
2106	wgaci(i,j)=min(y1(i,j), wgaci(i,j))
2107	!
2108	y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) &
2109	+dgaci(i,j)+wgaci(i,j))*d2t
2110	qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t
2111
2112	if (qi(i,j).lt.0.0) then
2113	a2=1.
2114	if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1.
2115	psaut(i,j)=psaut(i,j)*a2
2116	psaci(i,j)=psaci(i,j)*a2
2117	praci(i,j)=praci(i,j)*a2
2118	psfi(i,j)=psfi(i,j)*a2
2119	dgaci(i,j)=dgaci(i,j)*a2
2120	wgaci(i,j)=wgaci(i,j)*a2
2121	qi(i,j)=0.0
2122	endif
2123	!
2124	dlt3(i,j)=0.0
2125	dlt2(i,j)=0.0
2126	!
2127
2128	! DLT4(I,J)=1.0
2129	! if(qc(i,j) .gt. 5.e-4) dlt4(i,j)=0.0
2130	! if(qs(i,j) .le. 1.e-4) dlt4(i,j)=1.0
2131	!
2132	! IF (TAIR(I,J).ge.T0) THEN
2133	! DLT4(I,J)=0.0
2134	! ENDIF
2135
2136	if (tair(i,j).lt.t0) then
2137	if (qr(i,j).lt.1.e-4) then
2138	dlt3(i,j)=1.0
2139	dlt2(i,j)=1.0
2140	endif
2141	if (qs(i,j).ge.1.e-4) then
2142	dlt2(i,j)=0.0
2143	endif
2144	endif
2145
2146	if (ice2 .eq. 1) then
2147	dlt3(i,j)=1.0
2148	dlt2(i,j)=1.0
2149	endif
2150	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2151	pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t
2152	ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) &
2153	+psfi(i,j)+dlt3(i,j)praci(i,j))d2t
2154	! PS(I,J)=(PSAUT(I,J)+PSACI(I,J)+dlt4(i,j)*PSACW(I,J)
2155	! 1 +PSFW(I,J)+PSFI(I,J)+DLT3(I,J)PRACI(I,J))D2T
2156	pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) &
2157	+dgacw(i,j))*d2t
2158	! PG(I,J)=((1.-DLT3(I,J))*PRACI(I,J)+DGACI(I,J)+WGACI(I,J)
2159	! 1 +DGACW(I,J)+(1.-dlt4(i,j))PSACW(I,J))D2T
2160
2161	!JJS 150 CONTINUE
2162
2163	!* 7 * PRACS : ACCRETION OF QS BY QR *7
2164	!* 8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT) *8
2165
2166	!JJS DO 175 J=3,JLES
2167	!JJS DO 175 I=3,ILES
2168
2169	y1(i,j)=abs(vr(i,j)-vs(i,j))
2170	y2(i,j)=zr(i,j)*zs(i,j)
2171	y3(i,j)=5./y2(i,j)
2172	y4(i,j)=.08y3(i,j)y3(i,j)
2173	y5(i,j)=.05y3(i,j)y4(i,j)
2174	pracs(i,j)=r2igr2isr7ry1(i,j)(y3(i,j)/zs(i,j)**5 &
2175	+y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j))
2176	psacr(i,j)=r2isr8ry1(i,j)(y3(i,j)/zr(i,j)*5 &
2177	+y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j))
2178	qsacr(i,j)=psacr(i,j)
2179
2180	if (tair(i,j).ge.t0) then
2181	pracs(i,j)=0.0
2182	psacr(i,j)=0.0
2183	else
2184	qsacr(i,j)=0.0
2185	endif
2186	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2187	!* 2 * PGAUT : AUTOCONVERSION OF QS TO QG *2
2188	!* 18 * PGFR : FREEZING OF QR TO QG 18
2189
2190	pgaut(i,j)=0.0
2191	pgfr(i,j)=0.0
2192
2193	if (tair(i,j) .lt. t0) then
2194	! Y1(I,J)=EXP(.09*TAIRC(I,J))
2195	! PGAUT(I,J)=r2ismax(RN2Y1(I,J)*(QS(I,J)-BND2), 0.0)
2196	! IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2Y1(I,J)(QS(I,J)-BND2),0.0)
2197	y2(i,j)=exp(rn18a*(t0-tair(i,j)))
2198	!JJS PGFR(I,J)=r2igmax(R18R(Y2(I,J)-1.)/ZR(I,J)**7., 0.0)
2199	! pgfr(i,j)=r2icemax(r18r(y2(i,j)-1.)* &
2200	! (zr(i,j)**(-7.)), 0.0)
2201	! modify to prevent underflow on some computers (JD)
2202	temp = 1./zr(i,j)
2203	temp = temptemptemptemptemptemptemp
2204	pgfr(i,j)=r2igmax(r18r(y2(i,j)-1.)* &
2205	temp, 0.0)
2206	endif
2207
2208	!JJS 175 CONTINUE
2209
2210	!****** HANDLING THE NEGATIVE RAIN WATER (QR) *****************
2211	!****** HANDLING THE NEGATIVE SNOW (QS) *****************
2212
2213	!JJS DO 200 J=3,JLES
2214	!JJS DO 200 I=3,ILES
2215
2216	y1(i,j)=qr(i,j)/d2t
2217	y2(i,j)=-qg(i,j)/d2t
2218	piacr(i,j)=min(y1(i,j), piacr(i,j))
2219	dgacr(i,j)=min(y1(i,j), dgacr(i,j))
2220	wgacr(i,j)=min(y1(i,j), wgacr(i,j))
2221	wgacr(i,j)=max(y2(i,j), wgacr(i,j))
2222	psacr(i,j)=min(y1(i,j), psacr(i,j))
2223	pgfr(i,j)= min(y1(i,j), pgfr(i,j))
2224	del=0.
2225	if(wgacr(i,j) .lt. 0.) del=1.
2226	y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) &
2227	+psacr(i,j)+pgfr(i,j))*d2t
2228	qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-delwgacr(i,j)d2t
2229	if (qr(i,j) .lt. 0.0) then
2230	a1=1.
2231	if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1.
2232	piacr(i,j)=piacr(i,j)*a1
2233	dgacr(i,j)=dgacr(i,j)*a1
2234	if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1
2235	pgfr(i,j)=pgfr(i,j)*a1
2236	psacr(i,j)=psacr(i,j)*a1
2237	qr(i,j)=0.0
2238	endif
2239	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2240	prn(i,j)=d2t((1.-dlt3(i,j))piacr(i,j)+dgacr(i,j) &
2241	+wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j))
2242	ps(i,j)=ps(i,j)+d2t(dlt3(i,j)piacr(i,j) &
2243	+dlt2(i,j)*psacr(i,j))
2244	pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j)
2245	y1(i,j)=qs(i,j)/d2t
2246	pgacs(i,j)=min(y1(i,j), pgacs(i,j))
2247	dgacs(i,j)=min(y1(i,j), dgacs(i,j))
2248	wgacs(i,j)=min(y1(i,j), wgacs(i,j))
2249	pgaut(i,j)=min(y1(i,j), pgaut(i,j))
2250	pracs(i,j)=min(y1(i,j), pracs(i,j))
2251	psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) &
2252	+pgaut(i,j)+pracs(i,j))
2253	qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j)
2254
2255	if (qs(i,j).lt.0.0) then
2256	a2=1.
2257	if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1.
2258	pgacs(i,j)=pgacs(i,j)*a2
2259	dgacs(i,j)=dgacs(i,j)*a2
2260	wgacs(i,j)=wgacs(i,j)*a2
2261	pgaut(i,j)=pgaut(i,j)*a2
2262	pracs(i,j)=pracs(i,j)*a2
2263	psn(i,j)=psn(i,j)*a2
2264	qs(i,j)=0.0
2265	endif
2266	!
2267	!C PSN(I,J)=D2T*(PGACS(I,J)+DGACS(I,J)+WGACS(I,J)
2268	!c +PGAUT(I,J)+PRACS(I,J))
2269	y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) &
2270	+dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j))
2271	pt(i,j)=pt(i,j)+afcp*y2(i,j)
2272	qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j)
2273
2274	!JJS 200 CONTINUE
2275
2276	!* 11 * PSMLT : MELTING OF QS 11
2277	!* 19 * PGMLT : MELTING OF QG TO QR 19
2278
2279	!JJS DO 225 J=3,JLES
2280	!JJS DO 225 I=3,ILES
2281
2282	psmlt(i,j)=0.0
2283	pgmlt(i,j)=0.0
2284	tair(i,j)=(pt(i,j)+tb0)*pi0
2285
2286	if (tair(i,j).ge.t0) then
2287	tairc(i,j)=tair(i,j)-t0
2288	y1(i,j)=tca(i,j)tairc(i,j)-alvrdwv(i,j) &
2289	*(rp0-(qv(i,j)+qb0))
2290	y2(i,j)=.78/zs(i,j)*2+r101fscv(i,j)/zs(i,j)**bsh5
2291	dd(i,j)=r11rty1(i,j)y2(i,j)+r11at*tairc(i,j) &
2292	*(qsacw(i,j)+qsacr(i,j))
2293	psmlt(i,j)=r2is*max(0.0, min(dd(i,j), qs(i,j)))
2294
2295	if(ihail.eq.1) then
2296	y3(i,j)=.78/zg(i,j)*2+r19aqscv(i,j)/zg(i,j)**bgh5
2297	else
2298	y3(i,j)=.78/zg(i,j)*2+r19asscv(i,j)/zg(i,j)**bgh5
2299	endif
2300
2301	dd1(i,j)=r19rty1(i,j)y3(i,j)+r19bt*tairc(i,j) &
2302	*(qgacw(i,j)+qgacr(i,j))
2303	pgmlt(i,j)=r2ig*max(0.0, min(dd1(i,j), qg(i,j)))
2304	pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j))
2305	qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j)
2306	qs(i,j)=qs(i,j)-psmlt(i,j)
2307	qg(i,j)=qg(i,j)-pgmlt(i,j)
2308	endif
2309	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2310	!* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00) 24
2311	!* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00) 25
2312	!* 26 * PIMLT : MELTING OF QI TO QC (T >= T0) 26
2313
2314	if (qc(i,j).le.cmin1) qc(i,j)=0.0
2315	if (qi(i,j).le.cmin1) qi(i,j)=0.0
2316	tair(i,j)=(pt(i,j)+tb0)*pi0
2317
2318	if(tair(i,j).le.t00) then
2319	pihom(i,j)=qc(i,j)
2320	else
2321	pihom(i,j)=0.0
2322	endif
2323	if(tair(i,j).ge.t0) then
2324	pimlt(i,j)=qi(i,j)
2325	else
2326	pimlt(i,j)=0.0
2327	endif
2328	pidw(i,j)=0.0
2329
2330	if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then
2331	tairc(i,j)=tair(i,j)-t0
2332	y1(i,j)=max( min(tairc(i,j), -1.), -31.)
2333	it(i,j)=int(abs(y1(i,j)))
2334	y2(i,j)=aa1(it(i,j))
2335	y3(i,j)=aa2(it(i,j))
2336	y4(i,j)=exp(abs(beta*tairc(i,j)))
2337	y5(i,j)=(r00qi(i,j)/(r25ay4(i,j)))**y3(i,j)
2338	pidw(i,j)=min(r25rty2(i,j)y4(i,j)*y5(i,j), qc(i,j))
2339	endif
2340
2341	y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j)
2342	pt(i,j)=pt(i,j)+afcpy1(i,j)+ascp(pidep(i,j))*d2t
2343	qv(i,j)=qv(i,j)-(pidep(i,j))*d2t
2344	qc(i,j)=qc(i,j)-y1(i,j)
2345	qi(i,j)=qi(i,j)+y1(i,j)
2346
2347	!* 31 * PINT : INITIATION OF QI 31
2348	!* 32 * PIDEP : DEPOSITION OF QI 32
2349	!
2350	! CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR
2351	! THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER
2352	! CONCENTRATION EXCEEDS THAT ALREADY PRESENT.
2353	!* 31 * pint : initiation of qi 31
2354	!* 32 * pidep : deposition of qi 32
2355	pint(i,j)=0.0
2356	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2357	if ( itaobraun.eq.1 ) then
2358	tair(i,j)=(pt(i,j)+tb0)*pi0
2359	if (tair(i,j) .lt. t0) then
2360	! if (qi(i,j) .le. cmin) qi(i,j)=0.
2361	if (qi(i,j) .le. cmin2) qi(i,j)=0.
2362	tairc(i,j)=tair(i,j)-t0
2363	rtair(i,j)=1./(tair(i,j)-c76)
2364	y2(i,j)=exp(c218-c580*rtair(i,j))
2365	qsi(i,j)=rp0*y2(i,j)
2366	esi(i,j)=c610*y2(i,j)
2367	ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2368	ami50=3.76e-8
2369
2370	!ccif ( itaobraun.eq.1 ) --> betah=0.5beta=-.460.5=-0.23; cn0=1.e-6
2371	!ccif ( itaobraun.eq.0 ) --> betah=0.5beta=-.60.5=-0.30; cn0=1.e-8
2372
2373	y1(i,j)=1./tair(i,j)
2374
2375	!cc insert a restriction on ice collection that ice collection
2376	!cc should be stopped at -30 c (with cn0=1.e-6, beta=-.46)
2377
2378	tairccri=tairc(i,j) ! in degree c
2379	if(tairccri.le.-30.) tairccri=-30.
2380
2381	! y2(i,j)=exp(betah*tairc(i,j))
2382	y2(i,j)=exp(betah*tairccri)
2383	y3(i,j)=sqrt(qi(i,j))
2384	dd(i,j)=y1(i,j)(rn10ay1(i,j)-rn10b)+rn10c*tair(i,j) &
2385	/esi(i,j)
2386	pidep(i,j)=max(r32rtssi(i,j)y2(i,j)*y3(i,j)/dd(i,j), 0.e0)
2387
2388	r_nci=min(cn0exp(betatairc(i,j)),1.)
2389	! r_nci=min(1.e-6exp(-.46tairc(i,j)),1.)
2390
2391	dd(i,j)=max(1.e-9r_nci/r00-qi(i,j)1.e-9/ami50, 0.)
2392	dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0)
2393	rsub1(i,j)=cs580qsi(i,j)rtair(i,j)*rtair(i,j)
2394	dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
2395	pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
2396
2397	! pint(i,j)=min(pint(i,j), dep(i,j))
2398	pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j))
2399
2400	! if (pint(i,j) .le. cmin) pint(i,j)=0.
2401	if (pint(i,j) .le. cmin2) pint(i,j)=0.
2402	pt(i,j)=pt(i,j)+ascp*pint(i,j)
2403	qv(i,j)=qv(i,j)-pint(i,j)
2404	qi(i,j)=qi(i,j)+pint(i,j)
2405	endif
2406	endif ! if ( itaobraun.eq.1 )
2407	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2408	!
2409	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2410	if ( itaobraun.eq.0 ) then
2411	tair(i,j)=(pt(i,j)+tb0)*pi0
2412	if (tair(i,j) .lt. t0) then
2413	if (qi(i,j) .le. cmin1) qi(i,j)=0.
2414	tairc(i,j)=tair(i,j)-t0
2415	dd(i,j)=r31rexp(betatairc(i,j))
2416	rtair(i,j)=1./(tair(i,j)-c76)
2417	y2(i,j)=exp(c218-c580*rtair(i,j))
2418	qsi(i,j)=rp0*y2(i,j)
2419	esi(i,j)=c610*y2(i,j)
2420	ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2421	dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
2422	rsub1(i,j)=cs580qsi(i,j)rtair(i,j)*rtair(i,j)
2423	dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
2424	pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
2425	y1(i,j)=1./tair(i,j)
2426	y2(i,j)=exp(betah*tairc(i,j))
2427	y3(i,j)=sqrt(qi(i,j))
2428	dd(i,j)=y1(i,j)(rn10ay1(i,j)-rn10b) &
2429	+rn10c*tair(i,j)/esi(i,j)
2430	pidep(i,j)=max(r32rtssi(i,j)y2(i,j)*y3(i,j)/dd(i,j), 0.)
2431	pint(i,j)=pint(i,j)+pidep(i,j)
2432	pint(i,j)=min(pint(i,j),dep(i,j))
2433	!c if (pint(i,j) .le. cmin2) pint(i,j)=0.
2434	pt(i,j)=pt(i,j)+ascp*pint(i,j)
2435	qv(i,j)=qv(i,j)-pint(i,j)
2436	qi(i,j)=qi(i,j)+pint(i,j)
2437	endif
2438	endif ! if ( itaobraun.eq.0 )
2439	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2440
2441	!JJS 225 CONTINUE
2442
2443	!*** TAO ET AL (1989) SATURATION TECHNIQUE *********************
2444
2445	if (new_ice_sat .eq. 0) then
2446
2447	!JJS DO 250 J=3,JLES
2448	!JJS DO 250 I=3,ILES
2449	tair(i,j)=(pt(i,j)+tb0)*pi0
2450	cnd(i,j)=rt0*(tair(i,j)-t00)
2451	dep(i,j)=rt0*(t0-tair(i,j))
2452	y1(i,j)=1./(tair(i,j)-c358)
2453	y2(i,j)=1./(tair(i,j)-c76)
2454	qsw(i,j)=rp0exp(c172-c409y1(i,j))
2455	qsi(i,j)=rp0exp(c218-c580y2(i,j))
2456	dd(i,j)=cp409y1(i,j)y1(i,j)
2457	dd1(i,j)=cp580y2(i,j)y2(i,j)
2458	if (qc(i,j).le.cmin) qc(i,j)=cmin
2459	if (qi(i,j).le.cmin) qi(i,j)=cmin
2460	if (tair(i,j).ge.t0) then
2461	dep(i,j)=0.0
2462	cnd(i,j)=1.
2463	qi(i,j)=0.0
2464	endif
2465
2466	if (tair(i,j).lt.t00) then
2467	cnd(i,j)=0.0
2468	dep(i,j)=1.
2469	qc(i,j)=0.0
2470	endif
2471
2472	y5(i,j)=avcpcnd(i,j)+ascpdep(i,j)
2473	! if (qc(i,j) .ge. cmin .or. qi(i,j) .ge. cmin) then
2474	y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j))
2475	y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j))
2476	y4(i,j)=dd(i,j)y1(i,j)+dd1(i,j)y2(i,j)
2477	qvs(i,j)=y1(i,j)+y2(i,j)
2478	rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
2479	cnd(i,j)=cnd(i,j)*rsub1(i,j)
2480	dep(i,j)=dep(i,j)*rsub1(i,j)
2481	if (qc(i,j).le.cmin) qc(i,j)=0.
2482	if (qi(i,j).le.cmin) qi(i,j)=0.
2483	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2484	!c **** condensation or evaporation of qc ****
2485
2486	cnd(i,j)=max(-qc(i,j),cnd(i,j))
2487
2488	!c **** deposition or sublimation of qi ****
2489
2490	dep(i,j)=max(-qi(i,j),dep(i,j))
2491
2492	pt(i,j)=pt(i,j)+avcpcnd(i,j)+ascpdep(i,j)
2493	qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
2494	qc(i,j)=qc(i,j)+cnd(i,j)
2495	qi(i,j)=qi(i,j)+dep(i,j)
2496	!JJS 250 continue
2497	endif
2498
2499	if (new_ice_sat .eq. 1) then
2500
2501	!JJS DO J=3,JLES
2502	!JJS DO I=3,ILES
2503
2504	tair(i,j)=(pt(i,j)+tb0)*pi0
2505	cnd(i,j)=rt0*(tair(i,j)-t00)
2506	dep(i,j)=rt0*(t0-tair(i,j))
2507	y1(i,j)=1./(tair(i,j)-c358)
2508	y2(i,j)=1./(tair(i,j)-c76)
2509	qsw(i,j)=rp0exp(c172-c409y1(i,j))
2510	qsi(i,j)=rp0exp(c218-c580y2(i,j))
2511	dd(i,j)=cp409y1(i,j)y1(i,j)
2512	dd1(i,j)=cp580y2(i,j)y2(i,j)
2513	y5(i,j)=avcpcnd(i,j)+ascpdep(i,j)
2514	y1(i,j)=rt0(tair(i,j)-t00)qsw(i,j)
2515	y2(i,j)=rt0(t0-tair(i,j))qsi(i,j)
2516	! IF (QC(I,J).LE.CMIN) QC(I,J)=CMIN
2517	! IF (QI(I,J).LE.CMIN) QI(I,J)=CMIN
2518
2519	if (tair(i,j).ge.t0) then
2520	! QI(I,J)=0.0
2521	dep(i,j)=0.0
2522	cnd(i,j)=1.
2523	y2(i,j)=0.
2524	y1(i,j)=qsw(i,j)
2525	endif
2526	if (tair(i,j).lt.t00) then
2527	cnd(i,j)=0.0
2528	dep(i,j)=1.
2529	y2(i,j)=qsi(i,j)
2530	y1(i,j)=0.
2531	! QC(I,J)=0.0
2532	endif
2533
2534	! Y1(I,J)=QC(I,J)*QSW(I,J)/(QC(I,J)+QI(I,J))
2535	! Y2(I,J)=QI(I,J)*QSI(I,J)/(QC(I,J)+QI(I,J))
2536
2537	y4(i,j)=dd(i,j)y1(i,j)+dd1(i,j)y2(i,j)
2538	qvs(i,j)=y1(i,j)+y2(i,j)
2539	rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
2540	cnd(i,j)=cnd(i,j)*rsub1(i,j)
2541	dep(i,j)=dep(i,j)*rsub1(i,j)
2542	! IF (QC(I,J).LE.CMIN) QC(I,J)=0.
2543	! IF (QI(I,J).LE.CMIN) QI(I,J)=0.
2544
2545	!C **** CONDENSATION OR EVAPORATION OF QC ****
2546
2547	cnd(i,j)=max(-qc(i,j),cnd(i,j))
2548
2549	!C **** DEPOSITION OR SUBLIMATION OF QI ****
2550
2551	dep(i,j)=max(-qi(i,j),dep(i,j))
2552
2553	pt(i,j)=pt(i,j)+avcpcnd(i,j)+ascpdep(i,j)
2554	qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
2555	qc(i,j)=qc(i,j)+cnd(i,j)
2556	qi(i,j)=qi(i,j)+dep(i,j)
2557	!JJS ENDDO
2558	!JJS ENDDO
2559	endif
2560
2561	!c
2562	!
2563	if (new_ice_sat .eq. 2) then
2564	!JJS do j=3,jles
2565	!JJS do i=3,iles
2566	dep(i,j)=0.0
2567	cnd(i,j)=0.0
2568	tair(i,j)=(pt(i,j)+tb0)*pi0
2569	if (tair(i,j) .ge. 253.16) then
2570	y1(i,j)=1./(tair(i,j)-c358)
2571	qsw(i,j)=rp0exp(c172-c409y1(i,j))
2572	dd(i,j)=cp409y1(i,j)y1(i,j)
2573	dm(i,j)=qv(i,j)+qb0-qsw(i,j)
2574	cnd(i,j)=dm(i,j)/(1.+avcpdd(i,j)qsw(i,j))
2575	!c **** condensation or evaporation of qc ****
2576	cnd(i,j)=max(-qc(i,j), cnd(i,j))
2577	pt(i,j)=pt(i,j)+avcp*cnd(i,j)
2578	qv(i,j)=qv(i,j)-cnd(i,j)
2579	qc(i,j)=qc(i,j)+cnd(i,j)
2580	endif
2581	if (tair(i,j) .le. 258.16) then
2582	!c cnd(i,j)=0.0
2583	y2(i,j)=1./(tair(i,j)-c76)
2584	qsi(i,j)=rp0exp(c218-c580y2(i,j))
2585	dd1(i,j)=cp580y2(i,j)y2(i,j)
2586	dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascpdd1(i,j)qsi(i,j))
2587	!c **** deposition or sublimation of qi ****
2588	dep(i,j)=max(-qi(i,j),dep(i,j))
2589	pt(i,j)=pt(i,j)+ascp*dep(i,j)
2590	qv(i,j)=qv(i,j)-dep(i,j)
2591	qi(i,j)=qi(i,j)+dep(i,j)
2592	endif
2593	!JJS enddo
2594	!JJS enddo
2595	endif
2596
2597	!c
2598	!
2599	!* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS 10
2600	!* 20 * PGSUB : SUBLIMATION OF QG 20
2601
2602	!JJS DO 275 J=3,JLES
2603	!JJS DO 275 I=3,ILES
2604
2605	psdep(i,j)=0.0
2606	pgdep(i,j)=0.0
2607	pssub(i,j)=0.0
2608	pgsub(i,j)=0.0
2609	tair(i,j)=(pt(i,j)+tb0)*pi0
2610
2611	if(tair(i,j).lt.t0) then
2612	if(qs(i,j).lt.cmin1) qs(i,j)=0.0
2613	if(qg(i,j).lt.cmin1) qg(i,j)=0.0
2614	rtair(i,j)=1./(tair(i,j)-c76)
2615	qsi(i,j)=rp0exp(c218-c580rtair(i,j))
2616	ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2617	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2618	y1(i,j)=r10ar/(tca(i,j)tair(i,j)*2)+1./(dwv(i,j) &
2619	*qsi(i,j))
2620	y2(i,j)=.78/zs(i,j)*2+r101fscv(i,j)/zs(i,j)**bsh5
2621	psdep(i,j)=r10tssi(i,j)y2(i,j)/y1(i,j)
2622	pssub(i,j)=psdep(i,j)
2623	psdep(i,j)=r2is*max(psdep(i,j), 0.)
2624	pssub(i,j)=r2is*max(-qs(i,j), min(pssub(i,j), 0.))
2625
2626	if(ihail.eq.1) then
2627	y2(i,j)=.78/zg(i,j)*2+r20bqscv(i,j)/zg(i,j)**bgh5
2628	else
2629	y2(i,j)=.78/zg(i,j)*2+r20bsscv(i,j)/zg(i,j)**bgh5
2630	endif
2631
2632	pgsub(i,j)=r2igr20tssi(i,j)*y2(i,j)/y1(i,j)
2633	dm(i,j)=qv(i,j)+qb0-qsi(i,j)
2634	rsub1(i,j)=cs580qsi(i,j)rtair(i,j)*rtair(i,j)
2635
2636	! ****** DEPOSITION OR SUBLIMATION OF QS ********************
2637
2638	y1(i,j)=dm(i,j)/(1.+rsub1(i,j))
2639	psdep(i,j)=r2is*min(psdep(i,j),max(y1(i,j),0.))
2640	y2(i,j)=min(y1(i,j),0.)
2641	pssub(i,j)=r2is*max(pssub(i,j),y2(i,j))
2642
2643	! ****** SUBLIMATION OF QG *********************************
2644
2645	dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.)
2646	pgsub(i,j)=r2ig*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.))
2647
2648	if(qc(i,j)+qi(i,j).gt.1.e-5) then
2649	dlt1(i,j)=1.
2650	else
2651	dlt1(i,j)=0.
2652	endif
2653
2654	psdep(i,j)=dlt1(i,j)*psdep(i,j)
2655	pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j)
2656	pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j)
2657
2658	pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j))
2659	qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j)
2660	qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j)
2661	qg(i,j)=qg(i,j)-pgsub(i,j)
2662	endif
2663
2664	!* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) 23
2665
2666	ern(i,j)=0.0
2667
2668	if(qr(i,j).gt.0.0) then
2669	tair(i,j)=(pt(i,j)+tb0)*pi0
2670	rtair(i,j)=1./(tair(i,j)-c358)
2671	qsw(i,j)=rp0exp(c172-c409rtair(i,j))
2672	ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
2673	dm(i,j)=qv(i,j)+qb0-qsw(i,j)
2674	rsub1(i,j)=cv409qsw(i,j)rtair(i,j)*rtair(i,j)
2675	dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
2676	y1(i,j)=.78/zr(i,j)*2+r23afscv(i,j)/zr(i,j)**bwh5
2677	y2(i,j)=r23br/(tca(i,j)tair(i,j)*2)+1./(dwv(i,j) &
2678	*qsw(i,j))
2679	!cccc
2680	ern(i,j)=r23tssw(i,j)y1(i,j)/y2(i,j)
2681	ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
2682	pt(i,j)=pt(i,j)-avcp*ern(i,j)
2683	qv(i,j)=qv(i,j)+ern(i,j)
2684	qr(i,j)=qr(i,j)-ern(i,j)
2685	endif
2686
2687	!JJS 10/7/2008 vvvvv
2688	ENDIF ! part of if (iwarm.eq.1) then
2689	!JJS 10/7/2008 ^^^^^
2690
2691	! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
2692	if (qc(i,j) .le. cmin1) qc(i,j)=0.
2693	if (qr(i,j) .le. cmin1) qr(i,j)=0.
2694	if (qi(i,j) .le. cmin1) qi(i,j)=0.
2695	if (qs(i,j) .le. cmin1) qs(i,j)=0.
2696	if (qg(i,j) .le. cmin1) qg(i,j)=0.
2697	dpt(i,j,k)=pt(i,j)
2698	dqv(i,j,k)=qv(i,j)
2699	qcl(i,j,k)=qc(i,j)
2700	qrn(i,j,k)=qr(i,j)
2701	qci(i,j,k)=qi(i,j)
2702	qcs(i,j,k)=qs(i,j)
2703	qcg(i,j,k)=qg(i,j)
2704
2705	!JJS 275 CONTINUE
2706
2707	scc=0.
2708	see=0.
2709
2710	! DO 110 J=3,JLES
2711	! DO 110 I=3,ILES
2712	! DD(I,J)=MAX(-CND(I,J), 0.)
2713	! CND(I,J)=MAX(CND(I,J), 0.)
2714	! DD1(I,J)=MAX(-DEP(I,J), 0.)
2715
2716	!ccshie 2/21/02 shie follow tao
2717	!CC for reference QI(I,J)=QI(I,J)-Y2(I,J)+PIDEP(I,J)*D2T
2718	!CC for reference QV(I,J)=QV(I,J)-(PIDEP(I,J))*D2T
2719
2720	!c DEP(I,J)=MAX(DEP(I,J), 0.)
2721	! DEP(I,J)=MAX(DEP(I,J), 0.)+PIDEP(I,J)*D2T
2722	! SCC=SCC+CND(I,J)
2723	! SEE=SEE+DD(I,J)+ERN(I,J)
2724
2725	! 110 CONTINUE
2726
2727	! SC(K)=SCC+SC(K)
2728	! SE(K)=SEE+SE(K)
2729
2730	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2731	!c henry: please take a look (start)
2732	!JJS modified by JJS on 5/1/2007 vvvvv
2733
2734	!JJS do 305 j=3,jles
2735	!JJS do 305 i=3,iles
2736	dd(i,j)=max(-cnd(i,j), 0.)
2737	cnd(i,j)=max(cnd(i,j), 0.)
2738	dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t
2739	dep(i,j)=max(dep(i,j), 0.)
2740	!JJS 305 continue
2741
2742	!JJS do 306 j=3,jles
2743	!JJS do 306 i=3,iles
2744	!JJS scc=scc+cnd(i,j)
2745	!JJS see=see+(dd(i,j)+ern(i,j))
2746	!
2747	!JJS sddd=sddd+(dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j))
2748	!JJS ssss=ssss+dd1(i,j)
2749	!JJS
2750	! shhh=shhh+rsw(i,j,k)*d2t
2751	! sccc=sccc+rlw(i,j,k)*d2t
2752	!jjs
2753	!JJS smmm=smmm+(psmlt(i,j)+pgmlt(i,j)+pimlt(i,j))
2754	!JJS sfff=sfff+d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j)
2755	!JJS 1 +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j))
2756	!JJS 2 -qracs(i,j)+pihom(i,j)+pidw(i,j)
2757
2758
2759	sccc=cnd(i,j)
2760	seee=dd(i,j)+ern(i,j)
2761	sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j)
2762	ssss=dd1(i,j) + pgsub(i,j)
2763	smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j)
2764	sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) &
2765	+pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) &
2766	+wgacr(i,j))+pihom(i,j)+pidw(i,j)
2767
2768	! physc(i,k,j) = avcp * sccc / d2t
2769	! physe(i,k,j) = avcp * seee / d2t
2770	! physd(i,k,j) = ascp * sddd / d2t
2771	! physs(i,k,j) = ascp * ssss / d2t
2772	! physf(i,k,j) = afcp * sfff / d2t
2773	! physm(i,k,j) = afcp * smmm / d2t
2774	! physc(i,k,j) = physc(i,k,j) + avcp * sccc
2775	! physe(i,k,j) = physc(i,k,j) + avcp * seee
2776	! physd(i,k,j) = physd(i,k,j) + ascp * sddd
2777	! physs(i,k,j) = physs(i,k,j) + ascp * ssss
2778	! physf(i,k,j) = physf(i,k,j) + afcp * sfff
2779	! physm(i,k,j) = physm(i,k,j) + afcp * smmm
2780
2781	!JJS modified by JJS on 5/1/2007 ^^^^^
2782
2783	2000 continue
2784
2785	1000 continue
2786
2787	!JJS ****************************************************************
2788	!JJS convert from GCE grid back to WRF grid
2789	do k=kts,kte
2790	do j=jts,jte
2791	do i=its,ite
2792	ptwrf(i,k,j) = dpt(i,j,k)
2793	qvwrf(i,k,j) = dqv(i,j,k)
2794	qlwrf(i,k,j) = qcl(i,j,k)
2795	qrwrf(i,k,j) = qrn(i,j,k)
2796	qiwrf(i,k,j) = qci(i,j,k)
2797	qswrf(i,k,j) = qcs(i,j,k)
2798	qgwrf(i,k,j) = qcg(i,j,k)
2799	enddo !i
2800	enddo !j
2801	enddo !k
2802
2803	! ****************************************************************
2804
2805	return
2806	END SUBROUTINE saticel_s
2807
2808	!JJS
2809	!JJS REAL FUNCTION GAMMA(X)
2810	!JJS Y=GAMMLN(X)
2811	!JJS GAMMA=EXP(Y)
2812	!JJS RETURN
2813	!JJS END
2814	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2815	!JJS real function GAMMLN (xx)
2816	real function gammagce (xx)
2817	!**********************************************************************
2818	real*8 cof(6),stp,half,one,fpf,x,tmp,ser
2819	data cof,stp / 76.18009173,-86.50532033,24.01409822, &
2820	-1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 /
2821	data half,one,fpf / .5, 1., 5.5 /
2822	!
2823	x=xx-one
2824	tmp=x+fpf
2825	tmp=(x+half)*log(tmp)-tmp
2826	ser=one
2827	do j=1,6
2828	x=x+one
2829	ser=ser+cof(j)/x
2830	enddo !j
2831	gammln=tmp+log(stp*ser)
2832	!JJS
2833	gammagce=exp(gammln)
2834	!JJS
2835	return
2836	END FUNCTION gammagce
2837
2838	END MODULE module_mp_gsfcgce

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

Download in other formats:

Original Format