| blob | 991da22d11d81144b18864c33bef7b0d65a5b538 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 MODULE module_mp_gsfcgce
6 USE module_wrf_error
8 !JJS 1/3/2008 vvvvv
10 ! common /bt/
11 REAL, PRIVATE :: rd1, rd2, al, cp
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
22 ! common /size/
23 REAL, PRIVATE :: tnw, tns, tng, &
24 roqs, roqg, roqr
26 ! common /bterv/
27 REAL, PRIVATE :: zrc, zgc, zsc, &
28 vrc, vgc, vsc
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
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
44 REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a
46 ! common /rsnw1/
47 REAL, PRIVATE :: rn10b, rn10c, rnn191, rnn192, rn30, &
48 rnn30a, rn33, rn331, rn332
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/
67 !JJS 1/3/2008 ^^^^^
69 CONTAINS
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 )
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
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
134 REAL, DIMENSION( ims:ime , jms:jme ), &
135 INTENT(INOUT) :: rainnc, &
136 rainncv, &
137 snownc, &
138 snowncv, &
139 sr, &
140 graupelnc, &
141 graupelncv
146 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
148 REAL, INTENT(IN ) :: dt_in, &
149 grav, &
150 rhowater, &
151 rhosnow
153 LOGICAL, INTENT(IN), OPTIONAL :: F_QG
155 ! LOCAL VAR
158 !jjs INTEGER :: min_q, max_q
160 !jjs REAL, DIMENSION( its:ite , jts:jte ) &
161 !jjs :: rain, snow, graupel,ice
163 !
164 ! INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id
165 INTEGER :: itaobraun, istatmin, new_ice_sat, id
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
175 LOGICAL :: flag_qg
177 !
178 !c ihail = 0 for graupel, for tropical region
179 !c ihail = 1 for hail, for mid-lat region
181 ! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
182 !c if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
183 !c if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
184 itaobraun = 1
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
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
217 !c new_ice_sat = 0, 1 or 2
218 new_ice_sat = 2
220 !c istatmin
221 istatmin = 180
223 !c id = 0 without in-line staticstics
224 !c id = 1 with in-line staticstics
225 id = 0
227 !c ibud = 0 no calculation of dth, dqv, dqrest and dqall
228 !c ibud = 1 yes
229 ibud = 0
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
238 ! calculte fallflux and precipiation in MKS system
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 !-----------------------------------------------------------------------
251 !c set up constants used internally in GCE
253 call consat_s (ihail, itaobraun)
256 !c Negative values correction
258 iskip = 1
259 ! i12h=int(86400./dt_in)
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
285 !c microphysics in GCE
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 )
300 END SUBROUTINE gsfcgce
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 !-----------------------------------------------------------------------
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
331 REAL, INTENT(IN ) :: dt, grav, rhowater, rhosnow
334 REAL, DIMENSION( ims:ime , jms:jme ), &
335 INTENT(IN ) :: topo
337 ! temperary vars
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
348 INTEGER :: k, i, j
349 !
351 REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vti
353 REAL :: dtb, pi, consta, constc, gambp4, &
354 gamdp4, gam4pt5, gam4bbar
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.
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
374 REAL , PARAMETER :: rhoe_s = 1.29
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
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
395 !-----------------------------------------------------------------------
396 ! This program calculates precipitation fluxes due to terminal velocities.
397 !-----------------------------------------------------------------------
399 dtb=dt
400 pi=acos(-1.)
401 consta=2115.0*0.01**(1-constb)
402 ! constc=152.93*0.01**(1-constd)
403 constc=78.63*0.01**(1-constd)
405 ! Gamma function
406 gambp4=ggamma(constb+4.)
407 gamdp4=ggamma(constd+4.)
408 gam4pt5=ggamma(4.5)
409 gam4bbar=ggamma(4.+bbar)
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
418 j_loop: do j = jts, jte
419 i_loop: do i = its, ite
421 pptrain = 0.
422 pptsnow = 0.
423 pptgraul = 0.
424 pptice = 0.
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
436 DO k = kts, kte
437 qiz(k)=qi(i,k,j)
438 ENDDO
440 DO k = kts, kte
441 qsz(k)=qs(i,k,j)
442 ENDDO
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
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(pi*rhowater*xnor/rhoz(k)/qrz(k))
471 tmp1=sqrt(tmp1)
472 vtr(k)=consta*gambp4*sqrhoz(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
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)
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
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
523 !
524 !-- snow
525 !
526 t_del_tv=0.
527 del_tv=dtb
528 notlast=.true.
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(pi*rhosnow*xnos/rhoz(k)/qsz(k))
540 tmp1=sqrt(tmp1)
541 vts(k)=constc*gamdp4*sqrhoz(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
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)
559 t_del_tv=t_del_tv+del_tv
561 if ( t_del_tv .ge. dtb ) then
562 notlast=.false.
563 del_tv=dtb+del_tv-t_del_tv
564 endif
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
591 ENDDO
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
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(pi*rhohail*xnoh/rhoz(k)/qgz(k))
623 tmp1=sqrt(tmp1)
624 term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag)
625 vtg(k)=gam4pt5*term0*sqrt(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(pi*rhograul*xnog/rhoz(k)/qgz(k))
638 tmp1=sqrt(tmp1)
639 tmp1=tmp1**bbar
640 tmp1=1./tmp1
641 term0=abar*gam4bbar/6.
642 vtg(k)=term0*tmp1*(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
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)
660 t_del_tv=t_del_tv+del_tv
662 if ( t_del_tv .ge. dtb ) then
663 notlast=.false.
664 del_tv=dtb+del_tv-t_del_tv
665 endif
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 !
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
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)
728 t_del_tv=t_del_tv+del_tv
730 if ( t_del_tv .ge. dtb ) then
731 notlast=.false.
732 del_tv=dtb+del_tv-t_del_tv
733 endif
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 !
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
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 !
769 ! write(6,*) 'i=',i,' j=',j,' ', pptrain, pptsnow, pptgraul, pptice
770 ! call flush(6)
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)
781 ENDDO i_loop
782 ENDDO j_loop
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
795 END SUBROUTINE fall_flux
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
807 DATA B/-.577191652,.988205891,-.897056937,.918206857, &
808 -.756704078,.482199394,-.193527818,.035868343/
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
825 END FUNCTION ggamma
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
840 !c Local variables
841 ! REAL, DIMENSION( kts:kte ) :: Y1, Y2
842 REAL :: A0, A1, A2
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
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
862 A0=0.0
864 if (A1.NE.0.0.and.A1.GT.A2) then
865 A0=(A1-A2)/A1
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
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
905 END SUBROUTINE negcor
907 SUBROUTINE consat_s (ihail,itaobraun)
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
938 ! itaobraun=0 ! see Tao and Simpson (1993)
939 ! itaobraun=1 ! see Tao et al. (2003)
941 integer :: itaobraun
942 real :: cn0
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 ^^^^^
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=(cpi*roqr*tnw)**0.25
1058 zsc=(cpi*roqs*tns)**0.25
1059 zgc=(cpi*roqg*tng)**0.25
1060 vrc=aw*ga4b/(6.*zrc**bw)
1061 vsc=as*ga4d/(6.*zsc**bs)
1062 vgc=ag*ga4g/(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=.25*cpi*tns*cc1*ga3d
1072 esw=1.
1073 rn4=.25*cpi*esw*tns*cc1*ga3d
1074 ! ERI=1.
1075 eri=.1 ! 6/17/02 tao's ice913=0 (not 1)
1076 rn5=.25*cpi*eri*tnw*ac1*ga3b
1077 ! AMI=1./(24.*4.19E-10)
1078 ami=1./(24.*6.e-9) ! 6/15/02 tao's
1079 rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami
1080 ! ESR=1. ! also if ice913=1 for tao's
1081 esr=.5 ! 6/15/02 for ice913=0 tao's
1082 rn7=cpi2*esr*tnw*tns*roqs
1083 esr=1.
1084 rn8=cpi2*esr*tnw*tns*roqr
1085 rn9=cpi2*tns*tng*roqs
1086 rn10=2.*cpi*tns
1087 rn101=.31*ga5dh*sqrt(cc1)
1088 rn10a=als*als/rw
1089 !JJS
1090 rn10b=alv/tca
1091 rn10c=ars/(dwv*amw)
1092 !JJS
1093 rn11=2.*cpi*tns/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=cpi*eiw*ui50*ri50**2
1106 do 10 k=1,31
1107 y1=1.-aa2(k)
1108 rn13(k)=aa1(k)*y1/(ami50**y1-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
1114 egw=1.
1115 rn14=.25*cpi*egw*tng*ga3g*ag
1116 egi=.1
1117 rn15=.25*cpi*egi*tng*ga3g*ag
1118 egi=1.
1119 rn15a=.25*cpi*egi*tng*ga3g*ag
1120 egr=1.
1121 rn16=cpi2*egr*tng*tnw*roqr
1122 rn17=2.*cpi*tng
1123 rn17a=.31*ga5gh*sqrt(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.*cpi2*bpri*tnw*roqr
1130 rn18a=apri
1131 rn19=2.*cpi*tng/alf
1132 rn19a=.31*ga5gh*sqrt(ag)
1133 rn19b=cw/alf
1134 !
1135 rnn191=.78
1136 rnn192=.31*ga5gh*sqrt(ac2/dva)
1137 !
1138 rn20=2.*cpi*tng
1139 rn20a=als*als/rw
1140 rn20b=.31*ga5gh*sqrt(ag)
1141 bnd3=2.e-3
1142 rn21=1.e3*1.569e-12/0.15
1143 erw=1.
1144 rn22=.25*cpi*erw*ac1*tnw*ga3b
1145 rn23=2.*cpi*tnw
1146 rn23a=.31*ga5bh*sqrt(ac1)
1147 rn23b=alv*alv/rw
1150 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1151 !cc
1152 !cc "c0" in routine "consat" (2d), "consatrh" (3d)
1153 !cc if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
1154 !cc if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
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
1169 rn25=cn0
1170 rn30a=alv*als*amw/(tca*ars)
1171 rn30b=alv/tca
1172 rn30c=ars/(dwv*amw)
1173 rn31=1.e-17
1175 rn32=4.*51.545e-4
1176 !
1177 rn30=2.*cpi*tng
1178 rnn30a=alv*alv*amw/(tca*ars)
1179 !
1180 rn33=4.*tns
1181 rn331=.65
1182 rn332=.44*sqrt(ac3/dva)*ga5dh
1183 !
1185 return
1186 END SUBROUTINE consat_s
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
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
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
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
1250 !JJS the following common blocks have been moved to the top of
1251 !JJS module_mp_gsfcgce_driver_instat.F
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
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 ^^^^
1284 !JJS
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 ^^^^
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
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)
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
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
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
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
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
1391 !JJS common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny)
1393 !JJS COMMON/BLS/
1394 real, dimension (its:ite,jts:jte) :: &
1395 y0, ts0, qss0
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
1401 !JJS COMMON/MICRO/
1402 real, dimension (its:ite,kts:kte,jts:jte) :: &
1403 physc, physe, physd, &
1404 physs, physm, physf
1406 !23456789012345678901234567890123456789012345678901234567890123456789012
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
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/
1429 !JJS 1/3/2008 ^^^^^
1431 !
1432 save
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
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
1477 if (ice2 .eq. 2) ihail = 0
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
1488 !
1489 ! ****** THREE CLASSES OF ICE-PHASE (LIN ET AL, 1983) *********
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.97*roqs**.25/tns**.75
1529 ucog=687.97*roqg**.25/tng**.75
1530 uwet=4.464**.95
1532 rijl2 = 1. / (ide-ids) / (jde-jds)
1534 !JJScap $doacross local(j,i)
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
1544 f2=rd1*d2t
1545 f3=rd2*d2t
1547 ft=dt/d2t
1548 rft=rijl2*ft
1549 a0=.5*istatmin*rijl2
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
1567 ! ami50 for use in PINT
1568 ami50=3.76e-8
1569 ami100=1.51e-7
1570 ami40=2.41e-8
1572 !C ******************************************************************
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.
1582 do 2000 j=jts,jte
1583 do 2000 i=its,ite
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.e5*zrc*rrq
1602 zsr=1.e5*zsc*rrq
1603 zgr=1.e5*zgc*rrq
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=rn11*rr0*d2t
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=rn19*rr0*d2t
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=rn25*rr0*d2t
1648 r31r=rn31*rr0
1649 r32rt=rn32*d2t*rrs
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
1675 ! *** COMPUTE ZR,ZS,ZG,VR,VS,VG *****************************
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(vrcf*dd(i,j)**bwq, 0.)
1682 endif
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(vscf*dd(i,j)**bsq, 0.)
1689 endif
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(vgcr*dd(i,j)**bgq, 0.)
1697 else
1698 vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.)
1699 endif
1700 endif
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
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)
1712 y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.)
1713 dwv(i,j)=dwvp*tair(i,j)**1.81
1714 tca(i,j)=c141*y1(i,j)
1715 scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333)
1716 !JJS 100 CONTINUE
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**
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
1734 if (tair(i,j).lt.t0) then
1735 esi(i,j)=exp(.025*tairc(i,j))
1736 psaut(i,j)=r2iceg*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0)
1737 psaci(i,j)=r2iceg*r3f*esi(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)=R4F*QC(I,J)*DD(I,J)
1741 psacw(i,j)=r2iceg*0.5*r4f*qc(i,j)*dd(i,j)
1742 !JJS 3/30/06
1743 praci(i,j)=r2iceg*r5f*qi(i,j)/zr(i,j)**bw3
1744 piacr(i,j)=r2iceg*r6f*qi(i,j)*(zr(i,j)**(-bw6))
1745 !JJS PIACR(I,J)=R6F*QI(I,J)/ZR(I,J)**BW6
1746 else
1747 qsacw(i,j)=r2iceg*r4f*qc(i,j)*dd(i,j)
1748 endif
1750 !* 21 * PRAUT AUTOCONVERSION OF QC TO QR **21**
1751 !* 22 * PRACW : ACCRETION OF QC BY QR **22**
1753 pracw(i,j)=r22f*qc(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)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
1758 endif
1760 !* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971) **12**
1761 !* 13 * PSFI : BERGERON PROCESSES FOR QS **13**
1763 psfw(i,j)=0.0
1764 psfi(i,j)=0.0
1765 pidep(i,j)=0.0
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(d2t*y1(i,j)*(y2(i,j)+r12r*qc(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-6*exp(-.46*tairc(i,j)),1.)
1781 ! R_NCI=min(1.e-8*EXP(-.6*TAIRC(I,J)),1.) ! use Tao's
1782 dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
1783 rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
1784 y3(i,j)=1./tair(i,j)
1785 dd(i,j)=y3(i,j)*(rn30a*y3(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_nci*qi(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)=r2iceg*a2*.5*qi(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
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**
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)=ee1*exp(ee2*tairc(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)=.08*y3(i,j)*y3(i,j)
1828 y5(i,j)=.05*y3(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)=r2ice*r2iceg*r9r*egs(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)=r2ice*r2iceg*r9r*dd(i,j)
1840 ! WGACS(I,J)=0. ! 6/15/02 tao's
1841 y1(i,j)=1./zg(i,j)**bg3
1843 if(ihail .eq. 1) then
1844 dgacw(i,j)=r2ice*max(r14r*qc(i,j)*y1(i,j), 0.0)
1845 else
1846 dgacw(i,j)=r2ice*max(r14f*qc(i,j)*y1(i,j), 0.0)
1847 endif
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)=.08*y3(i,j)*y3(i,j)
1854 y5(i,j)=.05*y3(i,j)*y4(i,j)
1855 dd(i,j)=r16r*y1(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)
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
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**
1875 dgaci(i,j)=0.0
1876 wgaci(i,j)=0.0
1877 pgwet(i,j)=0.0
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)=r2ice*r15r*y1(i,j)
1883 wgaci(i,j)=r2ice*r15ar*y1(i,j)
1884 ! WGACI(I,J)=0. ! 6/15/02 tao's
1885 else
1887 !JJS DGACI(I,J)=r2ice*R15F*Y1(I,J)
1888 dgaci(i,j)=0.
1889 wgaci(i,j)=r2ice*r15af*y1(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+r17aq*scv(i,j)/zg(i,j)**bgh5
1900 else
1901 y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5
1902 endif
1903 y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) &
1904 -tca(i,j)*tairc(i,j)
1905 dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) &
1906 +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j)))
1907 pgwet(i,j)=r2ice*max(dd(i,j), 0.0)
1908 endif
1909 endif
1910 !JJS 125 CONTINUE
1912 !******** HANDLING THE NEGATIVE CLOUD WATER (QC) ******************
1913 !******** HANDLING THE NEGATIVE CLOUD ICE (QI) ******************
1915 !JJS DO 150 J=3,JLES
1916 !JJS DO 150 I=3,ILES
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))
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)
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
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 !
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
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
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
2021 !JJS 150 CONTINUE
2023 !* 7 * PRACS : ACCRETION OF QS BY QR ***7**
2024 !* 8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT) ***8**
2026 !JJS DO 175 J=3,JLES
2027 !JJS DO 175 I=3,ILES
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)=.08*y3(i,j)*y3(i,j)
2033 y5(i,j)=.05*y3(i,j)*y4(i,j)
2034 pracs(i,j)=r2ice*r2iceg*r7r*y1(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)=r2iceg*r8r*y1(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)
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**
2050 pgaut(i,j)=0.0
2051 pgfr(i,j)=0.0
2053 if (tair(i,j) .lt. t0) then
2054 ! Y1(I,J)=EXP(.09*TAIRC(I,J))
2055 ! PGAUT(I,J)=r2iceg*max(RN2*Y1(I,J)*(QS(I,J)-BND2), 0.0)
2056 ! IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2*Y1(I,J)*(QS(I,J)-BND2),0.0)
2057 y2(i,j)=exp(rn18a*(t0-tair(i,j)))
2058 !JJS PGFR(I,J)=r2ice*max(R18R*(Y2(I,J)-1.)/ZR(I,J)**7., 0.0)
2059 ! pgfr(i,j)=r2ice*max(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 = temp*temp*temp*temp*temp*temp*temp
2064 pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* &
2065 temp, 0.0)
2066 endif
2068 !JJS 175 CONTINUE
2070 !******** HANDLING THE NEGATIVE RAIN WATER (QR) *******************
2071 !******** HANDLING THE NEGATIVE SNOW (QS) *******************
2073 !JJS DO 200 J=3,JLES
2074 !JJS DO 200 I=3,ILES
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)-del*wgacr(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)
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)
2134 !JJS 200 CONTINUE
2136 !* 11 * PSMLT : MELTING OF QS **11**
2137 !* 19 * PGMLT : MELTING OF QG TO QR **19**
2139 !JJS DO 225 J=3,JLES
2140 !JJS DO 225 I=3,ILES
2142 psmlt(i,j)=0.0
2143 pgmlt(i,j)=0.0
2144 tair(i,j)=(pt(i,j)+tb0)*pi0
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)-alvr*dwv(i,j) &
2149 *(rp0-(qv(i,j)+qb0))
2150 y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
2151 dd(i,j)=r11rt*y1(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)))
2155 if(ihail.eq.1) then
2156 y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5
2157 else
2158 y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5
2159 endif
2161 dd1(i,j)=r19rt*y1(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**
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
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
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)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j)
2198 pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j))
2199 endif
2201 y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j)
2202 pt(i,j)=pt(i,j)+afcp*y1(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)
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
2230 !ccif ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
2231 !ccif ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
2233 y1(i,j)=1./tair(i,j)
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)
2238 tairccri=tairc(i,j) ! in degree c
2239 if(tairccri.le.-30.) tairccri=-30.
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)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) &
2245 /esi(i,j)
2246 pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0)
2248 r_nci=min(cn0*exp(beta*tairc(i,j)),1.)
2249 ! r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
2251 dd(i,j)=max(1.e-9*r_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)=cs580*qsi(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.)
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))
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)=r31r*exp(beta*tairc(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)=cs580*qsi(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)*(rn10a*y1(i,j)-rn10b) &
2289 +rn10c*tair(i,j)/esi(i,j)
2290 pidep(i,j)=max(r32rt*ssi(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
2301 !JJS 225 CONTINUE
2303 !***** TAO ET AL (1989) SATURATION TECHNIQUE ***********************
2305 if (new_ice_sat .eq. 0) then
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)=rp0*exp(c172-c409*y1(i,j))
2315 qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
2316 dd(i,j)=cp409*y1(i,j)*y1(i,j)
2317 dd1(i,j)=cp580*y2(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
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
2332 y5(i,j)=avcp*cnd(i,j)+ascp*dep(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 ******
2346 cnd(i,j)=max(-qc(i,j),cnd(i,j))
2348 !c ****** deposition or sublimation of qi ******
2350 dep(i,j)=max(-qi(i,j),dep(i,j))
2352 pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(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
2359 if (new_ice_sat .eq. 1) then
2361 !JJS DO J=3,JLES
2362 !JJS DO I=3,ILES
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)=rp0*exp(c172-c409*y1(i,j))
2370 qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
2371 dd(i,j)=cp409*y1(i,j)*y1(i,j)
2372 dd1(i,j)=cp580*y2(i,j)*y2(i,j)
2373 y5(i,j)=avcp*cnd(i,j)+ascp*dep(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
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
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))
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.
2405 !C ****** CONDENSATION OR EVAPORATION OF QC ******
2407 cnd(i,j)=max(-qc(i,j),cnd(i,j))
2409 !C ****** DEPOSITION OR SUBLIMATION OF QI ******
2411 dep(i,j)=max(-qi(i,j),dep(i,j))
2413 pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(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
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)=rp0*exp(c172-c409*y1(i,j))
2432 dd(i,j)=cp409*y1(i,j)*y1(i,j)
2433 dm(i,j)=qv(i,j)+qb0-qsw(i,j)
2434 cnd(i,j)=dm(i,j)/(1.+avcp*dd(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)=rp0*exp(c218-c580*y2(i,j))
2445 dd1(i,j)=cp580*y2(i,j)*y2(i,j)
2446 dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(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
2457 !c
2458 !
2459 !* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS **10**
2460 !* 20 * PGSUB : SUBLIMATION OF QG **20**
2462 !JJS DO 275 J=3,JLES
2463 !JJS DO 275 I=3,ILES
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
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)=rp0*exp(c218-c580*rtair(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+r101f*scv(i,j)/zs(i,j)**bsh5
2480 psdep(i,j)=r10t*ssi(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.))
2485 if(ihail.eq.1) then
2486 y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5
2487 else
2488 y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5
2489 endif
2491 pgsub(i,j)=r2ice*r20t*ssi(i,j)*y2(i,j)/y1(i,j)
2492 dm(i,j)=qv(i,j)+qb0-qsi(i,j)
2493 rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
2495 ! ******** DEPOSITION OR SUBLIMATION OF QS **********************
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))
2502 ! ******** SUBLIMATION OF QG ***********************************
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.))
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
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)
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
2523 !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) **23**
2525 ern(i,j)=0.0
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)=rp0*exp(c172-c409*rtair(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)=cv409*qsw(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+r23af*scv(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)=r23t*ssw(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
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)
2560 !JJS 275 CONTINUE
2562 scc=0.
2563 see=0.
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.)
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
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)
2580 ! 110 CONTINUE
2582 ! SC(K)=SCC+SC(K)
2583 ! SE(K)=SEE+SE(K)
2585 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2586 !c henry: please take a look (start)
2587 !JJS modified by JJS on 5/1/2007 vvvvv
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
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)
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)
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
2630 !JJS modified by JJS on 5/1/2007 ^^^^^
2632 2000 continue
2634 1000 continue
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
2652 ! ****************************************************************
2654 return
2655 END SUBROUTINE saticel_s
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
2687 END MODULE module_mp_gsfcgce