Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_gsfcgce.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 91.0 KB

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

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

Download in other formats:

Original Format