File: C:\NOAA\NEMS_11731\src\atmos\phys\calpreciptype.f

1           SUBROUTINE CALPRECIPTYPE(kdt,nrcm,im,ix,lm,lp1,randomno,  &
2                                xlat,xlon, &
3                                gt0,gq0,prsl,prsi,PREC, & !input
4                                phii,num_p3d,TSKIN,SR,phy_f3d,  & !input
5     			   DOMR,DOMZR,DOMIP,DOMS)  !output
6     !      SUBROUTINE CALPRECIPTYPE(nrcm,randomno,im,lm,lp1,T,Q,PMID,PINT,PREC, & !input
7     !                           ZINT,num_p3d,TSKIN,SR,F_RimeF,  & !input
8     !			   DOMR,DOMZR,DOMIP,DOMS)  !output
9     !$$$  SUBPROGRAM DOCUMENTATION BLOCK
10     !                .      .    .     
11     ! SUBPROGRAM:    CALPRECIPTYPE      COMPUTE DOMINANT PRECIP TYPE
12     !   PRGRMMR: CHUANG         ORG: W/NP2      DATE: 2008-05-28
13     !          
14     !     
15     ! ABSTRACT:
16     !     THIS ROUTINE COMPUTES PRECIPITATION TYPE.
17     !   . It is adopted from post but was made into a column to used by GFS model    
18     !     
19     !
20     !      use vrbls3d   
21     !      use vrbls2d   
22     !      use soil
23     !      use masks
24     !      use params_mod
25     !      use ctlblk_mod
26     !      use rqstfld_mod
27           USE FUNCPHYS, ONLY : fpvs,FTDP,fpkap,ftlcl,stma,fthe
28           USE PHYSCONS
29     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30           implicit none
31     !
32     !      INCLUDE "mpif.h"
33     !     
34     !     IN NGM SUBROUTINE OUTPUT WE FIND THE FOLLOWING COMMENT.
35     !     "IF THE FOLLOWING THRESHOLD VALUES ARE CHANGED, CONTACT
36     !     TDL/SYNOPTIC-SCALE TECHNIQUES BRANCH (PAUL DALLAVALLE
37     !     AND JOHN JENSENIUS).  THEY MAY BE USING IT IN ONE OF 
38     !     THEIR PACKING CODES."  THE THRESHOLD VALUE IS 0.01 INCH
39     !     OR 2.54E-4 METER.  PRECIPITATION VALUES LESS THAN THIS
40     !     THRESHOLD ARE SET TO MINUS ONE TIMES THIS THRESHOLD.
41           real,PARAMETER :: PTHRESH = 0.0
42     !     
43     !     SET CELCIUS TO KELVIN AND SECOND TO HOUR CONVERSION.
44           integer,PARAMETER :: NALG    = 5
45     !     
46     !     DECLARE VARIABLES.
47     !     
48           integer,intent(in) :: kdt,nrcm,im,ix,lm,lp1,num_p3d
49           real,intent(in) :: xlat(im),xlon(im) 
50           real(kind=kind_phys),dimension(im),intent(in) :: PREC,SR,TSKIN
51           real,intent(in) :: randomno(ix,nrcm)
52           real(kind=kind_phys),dimension(ix,LM),intent(in) :: gt0,gq0,prsl,phy_f3d
53           real(kind=kind_phys),dimension(ix,lp1),intent(in) :: prsi,phii
54           
55     !      real(kind=kind_phys),dimension(im,LM),intent(in) :: T,Q,PMID,F_RimeF
56     !      real(kind=kind_phys),dimension(im,lp1),intent(in) :: pint,zint
57           real(kind=kind_phys),dimension(im),intent(out) :: DOMR,DOMZR,DOMIP,DOMS
58           
59           INTEGER IWX1,IWX4,IWX5
60           REAL IWX2,IWX3
61           REAL(kind=kind_phys) ES,QC,PV
62           REAL SLEET(NALG),RAIN(NALG),FREEZR(NALG),SNOW(NALG)
63           real(kind=kind_phys),dimension(LM) :: T,Q,PMID,F_RimeF
64           real(kind=kind_phys),dimension(lp1) :: pint,zint
65           REAL(kind=kind_phys), ALLOCATABLE :: TWET(:),RH(:),TD(:)
66     !     REAL DOMS(IM,JM),DOMR(IM,JM),DOMIP(IM,JM),DOMZR(IM,JM) 
67     !
68           integer I,J,L,IWX,ISNO,IIP,IZR,IRAIN,k,k1
69           real(kind=kind_phys) tdpd,pr,tr,pk,tlcl,thelcl,qwet,               &
70                time_vert,time_ncep,time_ramer,time_bourg,time_revised,time_dominant &
71           ,btim,timef
72     !      real RDTPHS,TLOW,TSFCK,QSAT,DTOP,DBOT,SNEQV,RRNUM,SFCPRS,SFCQ,     &
73     !           RC,SFCTMP,SNCOVR,FACTRS,SOLAR
74     !     
75     !     computes wet bulb here since two algorithms use it
76     !      lp1=lm+1
77     ! convert geopotential to height
78     !      do l=1,lp1
79     !        zint(l)=zint(l)/con_g
80     !      end do
81     ! DON'T FORGET TO FLIP 3D ARRAYS AROUND BECAUSE GFS COUNTS FROM BOTTOM UP      
82           	    
83           ALLOCATE ( TWET(LM),RH(LM),TD(LM) )
84     !      print*,'debug calpreciptype: ', &
85     !      im,lm,lp1,nrcm
86     
87           time_vert    = 0.
88           time_ncep    = 0.
89           time_ramer   = 0.
90           time_bourg   = 0.
91           time_revised = 0.
92     
93           do i=1,im
94     !       if(kdt>15. and. kdt<20) btim = timef()
95             do k=1,lm
96     	  k1          = lm-k+1
97     	  t(k1)       = gt0(i,k)
98     	  q(k1)       = gq0(i,k)
99     	  pmid(k1)    = prsl(i,k) * 1000.0
100     	  f_rimef(k1) = phy_f3d(i,k) 
101     !
102     ! Compute wet bulb
103     !        do l=1,lm
104     	   pv     = pmid(k1)*q(k1)/(con_eps-con_epsm1*q(k1))
105     	   td(k1) = ftdp(pv)
106                tdpd   = t(k1)-td(k1)
107                if(pmid(k1)>=50000.)then ! only compute twet below 500mb to save time
108                  if(tdpd.gt.0.) then
109                    pr     = pmid(k1)
110                    tr     = t(k1)
111                    pk     = fpkap(pr)
112                    tlcl   = ftlcl(tr,tdpd)
113                    thelcl = fthe(tlcl,pk*tlcl/tr)
114                    call stma(thelcl,pk,twet(k1),qwet)
115                  else
116                    twet(k1)=t(k1)
117                  endif
118     	   endif 
119                ES     = FPVS(T(k1))
120     	   ES     = MIN(ES,PMID(k1))
121     	   QC     = CON_EPS*ES/(PMID(k1)+CON_EPSM1*ES)
122                RH(k1) = MAX(con_epsq,Q(k1))/QC
123     	  
124     	   k1       = lp1-k+1
125     	   pint(k1) = prsi(i,k) * 1000.0
126     	   zint(k1) = phii(i,k)/con_g
127     
128     	 enddo
129     	 pint(1) = prsi(i,lp1) * 1000.0
130     	 zint(1) = phii(i,lp1)/con_g
131     
132     !	 if(kdt>15.and.kdt<20) time_vert = time_vert + (timef() - btim)
133     ! debug print statement
134     !	if (abs(xlon(i)*57.29578-114.0) .lt. 0.2  .and. &
135     !	   abs(xlat(i)*57.29578-40.0) .lt. 0.2)then
136     !         print*,'debug in calpreciptype: i,im,lm,lp1,xlon,xlat,prec,tskin,sr,nrcm,randomno,num_p3d ', &
137     !         i,im,lm,lp1,xlon(i)*57.29578,xlat(i)*57.29578,prec(i),tskin(i),sr(i),  &
138     !	 nrcm,randomno(i,1:nrcm),num_p3d
139     !         do l=1,lm
140     !          print*,'debug in calpreciptype: l,t,q,p,pint,z,twet', &
141     !	  l,t(l),q(l), &
142     !          pmid(l),pint(l),zint(l),twet(l)
143     !         end do
144     !	 print*,'debug in calpreciptype: lp1,pint,z ', lp1,pint(lp1),zint(lp1)
145     !        end if  
146     ! end debug print statement		
147     !        CALL WETBULB(lm,con_rocp,con_epsq,T,Q,PMID,TWET)       
148     !     INSTANTANEOUS PRECIPITATION TYPE.
149     !        if(kdt>10.and.kdt<20)btim = timef() 
150     
151             CALL CALWXT(lm,lp1,T(1),Q(1),PMID(1),PINT(1),PREC(i),  &
152                         PTHRESH,con_fvirt,con_rog,con_epsq,   &
153                         ZINT(1),IWX1,TWET(1))
154     !        if(kdt>10.and.kdt<20)time_ncep=time_ncep+(timef() - btim)
155             IWX       = IWX1
156             ISNO      = MOD(IWX,2)
157             IIP       = MOD(IWX,4)/2
158             IZR       = MOD(IWX,8)/4
159             IRAIN     = IWX/8
160             SNOW(1)   = ISNO*1.0
161             SLEET(1)  = IIP*1.0
162             FREEZR(1) = IZR*1.0
163             RAIN(1)   = IRAIN*1.0
164     
165     
166     !     DOMINANT PRECIPITATION TYPE
167     !GSM  IF DOMINANT PRECIP TYPE IS REQUESTED, 4 MORE ALGORITHMS
168     !GSM    WILL BE CALLED.  THE TALLIES ARE THEN SUMMED IN
169     !GSM    CALWXT_DOMINANT
170     
171     !  RAMER ALGORITHM
172     !        ALLOCATE ( RH(LM),TD(LM) )
173     !        DO L=1,LM
174     !HC: use RH and TD consistent with GFS ice physics
175     !          ES=FPVS(T(L))
176     !	  ES=MIN(ES,PMID(L))
177     !	  QC=CON_EPS*ES/(PMID(L)+CON_EPSM1*ES)
178     !          RH(L)=MAX(con_epsq,Q(L))/QC
179     !	  PV   = PMID(L)*Q(L)/(CON_EPS-CON_EPSM1*Q(L))
180     !	  TD(L)=FTDP(PV)
181     !        END DO	
182     !        if(kdt>10.and.kdt<20)btim = timef()
183     
184             CALL CALWXT_RAMER(lm,lp1,T(1),Q(1),PMID(1),RH(1),TD(1), &
185     	                  PINT(1),PREC(i),PTHRESH,IWX2)
186     
187     !	if(kdt>10.and.kdt<20)time_ramer=time_ramer+(timef() - btim)
188     !        deallocate(RH,TD)         
189     !
190     !     DECOMPOSE IWX2 ARRAY
191     !
192             IWX       = NINT(IWX2)
193             ISNO      = MOD(IWX,2)
194             IIP       = MOD(IWX,4)/2
195             IZR       = MOD(IWX,8)/4
196             IRAIN     = IWX/8
197             SNOW(2)   = ISNO*1.0
198             SLEET(2)  = IIP*1.0
199             FREEZR(2) = IZR*1.0
200             RAIN(2)   = IRAIN*1.0
201     
202     ! BOURGOUIN ALGORITHM
203     !      ISEED=44641*(INT(SDAT(1)-1)*24*31+INT(SDAT(2))*24+IHRST)+   &
204     !     &  MOD(IFHR*60+IFMIN,44641)+4357
205     !        if(kdt>10.and.kdt<20)btim = timef()
206             CALL CALWXT_BOURG(LM,LP1,randomno(i,1),con_g,PTHRESH,           &
207          &                        T(1),Q(1),PMID(1),PINT(1),PREC(i),ZINT(1),IWX3)
208     !          print *,'in SURFCE,me=',me,'IWX3=',IWX3(1:30,JSTA),'PTHRESH=',PTHRESH
209     !        if(kdt>10.and.kdt<20)time_bourg=time_bourg+(timef() - btim)
210     !
211     !     DECOMPOSE IWX3 ARRAY
212     !
213             IWX       = NINT(IWX3)
214             ISNO      = MOD(IWX,2)
215             IIP       = MOD(IWX,4)/2
216             IZR       = MOD(IWX,8)/4
217             IRAIN     = IWX/8
218             SNOW(3)   = ISNO*1.0
219             SLEET(3)  = IIP*1.0
220             FREEZR(3) = IZR*1.0
221             RAIN(3)   = IRAIN*1.0
222     
223     !
224     ! REVISED NCEP ALGORITHM
225     !
226     !        if(kdt>10.and.kdt<20)btim = timef()
227     
228             CALL CALWXT_REVISED(LM,LP1,T(1),Q(1),PMID(1),PINT(1),PREC(i),PTHRESH,  &
229                                 con_fvirt,con_rog,con_epsq,ZINT(1),TWET(1),IWX4)
230     
231     !	if(kdt>10.and.kdt<20)time_revised=time_revised+(timef() - btim)	
232     !          print *,'in SURFCE,me=',me,'IWX4=',IWX4(1:30,JSTA)
233     !
234     !     DECOMPOSE IWX2 ARRAY
235     !
236             IWX       = IWX4
237             ISNO      = MOD(IWX,2)
238             IIP       = MOD(IWX,4)/2
239             IZR       = MOD(IWX,8)/4
240             IRAIN     = IWX/8
241             SNOW(4)   = ISNO*1.0
242             SLEET(4)  = IIP*1.0
243             FREEZR(4) = IZR*1.0
244             RAIN(4)   = IRAIN*1.0
245                   
246     ! EXPLICIT ALGORITHM (UNDER 18 NOT ADMITTED WITHOUT PARENT 
247     !     OR GUARDIAN)
248      
249             IF(num_p3d == 3) then ! Ferrier's scheme
250               CALL CALWXT_EXPLICIT(LM,PTHRESH,                           &
251     	  TSKIN(i),PREC(i),SR(i),F_RimeF(1),IWX5)
252             else
253               IWX5 = 0
254             endif
255     !     DECOMPOSE IWX2 ARRAY
256     !
257             IWX       = IWX5
258             ISNO      = MOD(IWX,2)
259             IIP       = MOD(IWX,4)/2
260             IZR       = MOD(IWX,8)/4
261             IRAIN     = IWX/8
262             SNOW(5)   = ISNO*1.0
263             SLEET(5)  = IIP*1.0
264             FREEZR(5) = IZR*1.0
265             RAIN(5)   = IRAIN*1.0
266     !               
267     !        if(kdt>10.and.kdt<20)btim = timef()
268     
269             CALL CALWXT_DOMINANT(NALG,PREC(i),PTHRESH,RAIN(1),FREEZR(1),SLEET(1), &
270                                 SNOW(1),DOMR(i),DOMZR(i),DOMIP(i),DOMS(i))
271     
272     !	if(kdt>10.and.kdt<20)time_dominant=time_dominant+(timef() - btim)       
273     !debug print statement
274     !        if (abs(xlon(i)*57.29578-114.0) .lt. 0.2  .and. &
275     !	   abs(xlat(i)*57.29578-40.0) .lt. 0.2) &
276     !         print*,'debug in calpreciptype: DOMR,DOMZR,DOMIP,DOMS ', &
277     !	  DOMR(i),DOMZR(i),DOMIP(i),DOMS(i)
278     ! end of debug print statement      
279     
280           enddo ! end loop for i
281     
282     !      if(kdt>10.and.kdt<20)print*, &
283     !      'time_vert,time_ncep,time_ramer,time_bourg,time_revised,time_dominant='&
284     !      ,time_vert,time_ncep,time_ramer,time_bourg,time_revised,time_dominant
285     
286           DEALLOCATE (TWET,RH,TD)        
287           RETURN
288           END
289     !
290     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
291     !
292            SUBROUTINE CALWXT(lm,lp1,T,Q,PMID,PINT,PREC,  &
293                      PTHRESH,D608,ROG,EPSQ,    &
294     		 ZINT,IWX,TWET)
295     ! 
296     !     FILE: CALWXT.f
297     !     WRITTEN: 11 NOVEMBER 1993, MICHAEL BALDWIN
298     !     REVISIONS:
299     !               30 SEPT 1994-SETUP NEW DECISION TREE (M BALDWIN)
300     !               12 JUNE 1998-CONVERSION TO 2-D (T BLACK)
301     !     01-10-25  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
302     !     02-01-15  MIKE BALDWIN - WRF VERSION
303     !                              
304     !
305     !     ROUTINE TO COMPUTE PRECIPITATION TYPE USING A DECISION TREE
306     !     APPROACH THAT USES VARIABLES SUCH AS INTEGRATED WET BULB TEMP
307     !     BELOW FREEZING AND LOWEST LAYER TEMPERATURE
308     !
309     !     SEE BALDWIN AND CONTORNO PREPRINT FROM 13TH WEATHER ANALYSIS
310     !     AND FORECASTING CONFERENCE FOR MORE DETAILS
311     !     (OR BALDWIN ET AL, 10TH NWP CONFERENCE PREPRINT)
312     ! 
313     !      use params_mod
314     !      use ctlblk_mod
315     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
316           implicit none
317     !
318     !    INPUT:
319     !      T,Q,PMID,HTM,LMH,PREC,ZINT
320     !
321           integer,intent(in):: lm,lp1
322     !      real,intent(in):: pthresh
323           real,dimension(LM),intent(in) :: T,Q,PMID,TWET
324           real,dimension(LP1),intent(in) :: ZINT,PINT
325           integer,intent(out)  :: IWX
326           real,intent(in) :: PREC,PTHRESH,D608,ROG,EPSQ
327     !      real,intent(out)  :: ZWET
328     
329     
330     !    OUTPUT:
331     !      IWX - INSTANTANEOUS WEATHER TYPE.
332     !        ACTS LIKE A 4 BIT BINARY
333     !          1111 = RAIN/FREEZING RAIN/ICE PELLETS/SNOW
334     !          WHERE THE ONE'S DIGIT IS FOR SNOW
335     !                THE TWO'S DIGIT IS FOR ICE PELLETS
336     !                THE FOUR'S DIGIT IS FOR FREEZING RAIN
337     !            AND THE EIGHT'S DIGIT IS FOR RAIN
338     !
339     !    INTERNAL:
340     !
341     !      REAL, ALLOCATABLE :: TWET(:)
342           real, parameter :: D00=0.0 
343           integer KARR,LICEE
344           real TCOLD,TWARM
345     
346     !    SUBROUTINES CALLED:
347     !     WETBULB
348     !     
349     !
350     !     INITIALIZE WEATHER TYPE ARRAY TO ZERO (IE, OFF).
351     !     WE DO THIS SINCE WE WANT IWX TO REPRESENT THE
352     !     INSTANTANEOUS WEATHER TYPE ON RETURN.
353     !     
354     !
355     !     ALLOCATE LOCAL STORAGE
356     !
357     
358           integer I,J,L,LMHK,LICE,IFREL,IWRML,IFRZL
359           real PSFCK,TDCHK,A,TDKL,TDPRE,TLMHK,TWRMK,AREAS8,AREAP4,   &
360                SURFW,SURFC,DZKL,AREA1,PINTK1,PINTK2,PM150,PKL,TKL,QKL
361     
362     !      ALLOCATE ( TWET(LM) )
363     !
364     !!$omp  parallel do
365           IWX = 0
366     !      ZWET=SPVAL
367     !
368     !!$omp  parallel do
369     !!$omp& private(a,lmhk,pkl,psfck,qkl,tdchk,tdkl,tdpre,tkl)
370     
371     !
372     !   SKIP THIS POINT IF NO PRECIP THIS TIME STEP 
373     !
374           IF (PREC.LE.PTHRESH) GOTO 800
375     !
376     !   FIND COLDEST AND WARMEST TEMPS IN SATURATED LAYER BETWEEN
377     !   70 MB ABOVE GROUND AND 500 MB
378     !   ALSO FIND HIGHEST SATURATED LAYER IN THAT RANGE
379     !
380     !meb
381           PSFCK=PINT(LM+1)
382     !meb
383           TDCHK=2.0
384       760 TCOLD=T(LM)
385           TWARM=T(LM)
386           LICEE=LM
387     !
388           DO 775 L=1,LM
389             QKL=Q(L)
390             QKL=AMAX1(EPSQ,QKL)
391             TKL=T(L)
392             PKL=PMID(L)
393     !
394     !   SKIP PAST THIS IF THE LAYER IS NOT BETWEEN 70 MB ABOVE GROUND
395     !       AND 500 MB
396     !
397             IF (PKL.LT.50000.0.OR.PKL.GT.PSFCK-7000.0) GOTO 775
398             A=ALOG(QKL*PKL/(6.1078*(0.378*QKL+0.622)))
399             TDKL=(237.3*A)/(17.269-A)+273.15
400             TDPRE=TKL-TDKL
401             IF (TDPRE.LT.TDCHK.AND.TKL.LT.TCOLD) TCOLD=TKL
402             IF (TDPRE.LT.TDCHK.AND.TKL.GT.TWARM) TWARM=TKL
403             IF (TDPRE.LT.TDCHK.AND.L.LT.LICEE) LICEE=L
404       775 CONTINUE
405     !
406     !    IF NO SAT LAYER AT DEW POINT DEP=TDCHK, INCREASE TDCHK
407     !     AND START AGAIN (BUT DON'T MAKE TDCHK > 6)
408     !
409           IF (TCOLD==T(LM).AND.TDCHK<6.0) THEN
410             TDCHK=TDCHK+2.0
411             GOTO 760
412           ENDIF
413       800 CONTINUE
414     !
415     !    LOWEST LAYER T
416     !
417           KARR=0
418           IF (PREC.LE.PTHRESH) GOTO 850
419           TLMHK=T(LM)
420     !
421     !    DECISION TREE TIME
422     !
423           IF (TCOLD>269.15) THEN
424               IF (TLMHK.LE.273.15) THEN
425     !             TURN ON THE FLAG FOR
426     !             FREEZING RAIN = 4
427     !             IF ITS NOT ON ALREADY
428     !             IZR=MOD(IWX(I,J),8)/4
429     !             IF (IZR.LT.1) IWX(I,J)=IWX(I,J)+4
430                 IWX=IWX+4
431                 GOTO 850
432               ELSE
433     !             TURN ON THE FLAG FOR
434     !             RAIN = 8
435     !             IF ITS NOT ON ALREADY
436     !             IRAIN=IWX(I,J)/8
437     !             IF (IRAIN.LT.1) IWX(I,J)=IWX(I,J)+8
438                 IWX=IWX+8
439                 GOTO 850
440               ENDIF
441           ENDIF
442           KARR=1
443       850 CONTINUE
444     !
445     !   COMPUTE WET BULB ONLY AT POINTS THAT NEED IT
446     !
447     !      CALL WETBULB(lm,T,Q,PMID,KARR,TWET)
448     !      CALL WETFRZLVL(TWET,ZWET)
449     !
450     !!$omp  parallel do
451     !!$omp& private(area1,areap4,areas8,dzkl,ifrzl,iwrml,lice,
452     !!$omp&         lmhk,pintk1,pintk2,pm150,psfck,surfc,surfw,
453     !!$omp&         tlmhk,twrmk)
454     
455           IF(KARR.GT.0)THEN
456             LICE=LICEE
457     !meb
458             PSFCK=PINT(LM+1)
459     !meb
460             TLMHK=T(LM)
461             TWRMK=TWARM
462     !
463     !    TWET AREA VARIABLES
464     !     CALCULATE ONLY WHAT IS NEEDED
465     !      FROM GROUND TO 150 MB ABOVE SURFACE
466     !      FROM GROUND TO TCOLD LAYER
467     !      AND FROM GROUND TO 1ST LAYER WHERE WET BULB T < 0.0
468     !
469     !     PINTK1 IS THE PRESSURE AT THE BOTTOM OF THE LAYER
470     !     PINTK2 IS THE PRESSURE AT THE TOP OF THE LAYER
471     !
472     !     AREAP4 IS THE AREA OF TWET ABOVE -4 C BELOW HIGHEST SAT LYR 
473     !
474             AREAS8=D00
475             AREAP4=D00
476             SURFW =D00
477             SURFC =D00
478     !
479             DO 1945 L=LM,LICE,-1
480             DZKL=ZINT(L)-ZINT(L+1)
481             AREA1=(TWET(L)-269.15)*DZKL
482             IF (TWET(L).GE.269.15) AREAP4=AREAP4+AREA1
483      1945   CONTINUE
484     !
485             IF (AREAP4.LT.3000.0) THEN
486     !             TURN ON THE FLAG FOR
487     !             SNOW = 1
488     !             IF ITS NOT ON ALREADY
489     !             ISNO=MOD(IWX(I,J),2)
490     !             IF (ISNO.LT.1) IWX(I,J)=IWX(I,J)+1
491               IWX=IWX+1
492               GO TO 1900
493             ENDIF
494     !
495     !     AREAS8 IS THE NET AREA OF TWET W.R.T. FREEZING IN LOWEST 150MB
496     !
497             PINTK1=PSFCK
498             PM150=PSFCK-15000.
499     !
500             DO 1955 L=LM,1,-1
501             PINTK2=PINT(L)
502             IF(PINTK1.LT.PM150)GO TO 1950
503             DZKL=ZINT(L)-ZINT(L+1)
504     !
505     !    SUM PARTIAL LAYER IF IN 150 MB AGL LAYER
506     !
507             IF(PINTK2.LT.PM150)                                      &
508               DZKL=T(L)*(Q(L)*D608+1.0)*ROG*ALOG(PINTK1/PM150)
509             AREA1=(TWET(L)-273.15)*DZKL
510             AREAS8=AREAS8+AREA1
511      1950   PINTK1=PINTK2
512      1955   CONTINUE
513     !
514     !     SURFW IS THE AREA OF TWET ABOVE FREEZING BETWEEN THE GROUND
515     !       AND THE FIRST LAYER ABOVE GROUND BELOW FREEZING
516     !     SURFC IS THE AREA OF TWET BELOW FREEZING BETWEEN THE GROUND
517     !       AND THE WARMEST SAT LAYER
518     !
519             IFRZL=0
520             IWRML=0
521     !
522             DO 2050 L=LM,1,-1
523             IF (IFRZL.EQ.0.AND.T(L).LT.273.15) IFRZL=1
524             IF (IWRML.EQ.0.AND.T(L).GE.TWRMK) IWRML=1
525     !
526             IF (IWRML.EQ.0.OR.IFRZL.EQ.0) THEN
527     !	  if(pmid(l) < 50000.)print*,'need twet above 500mb'
528               DZKL=ZINT(L)-ZINT(L+1)
529               AREA1=(TWET(L)-273.15)*DZKL
530               IF(IFRZL.EQ.0.AND.TWET(L).GE.273.15)SURFW=SURFW+AREA1
531               IF(IWRML.EQ.0.AND.TWET(L).LE.273.15)SURFC=SURFC+AREA1
532             ENDIF
533      2050   CONTINUE
534             IF(SURFC.LT.-3000.0.OR.   &
535               (AREAS8.LT.-3000.0.AND.SURFW.LT.50.0)) THEN
536     !             TURN ON THE FLAG FOR
537     !             ICE PELLETS = 2
538     !             IF ITS NOT ON ALREADY
539     !             IIP=MOD(IWX(I,J),4)/2
540     !             IF (IIP.LT.1) IWX(I,J)=IWX(I,J)+2
541               IWX=IWX+2
542               GOTO 1900
543             ENDIF
544     !
545             IF(TLMHK.LT.273.15) THEN
546     !             TURN ON THE FLAG FOR
547     !             FREEZING RAIN = 4
548     !             IF ITS NOT ON ALREADY
549     !             IZR=MOD(IWX(K),8)/4
550     !             IF (IZR.LT.1) IWX(K)=IWX(K)+4
551               IWX=IWX+4
552             ELSE
553     !             TURN ON THE FLAG FOR
554     !             RAIN = 8
555     !             IF ITS NOT ON ALREADY
556     !             IRAIN=IWX(K)/8
557     !             IF (IRAIN.LT.1) IWX(K)=IWX(K)+8
558               IWX=IWX+8
559             ENDIF
560           ENDIF
561      1900 CONTINUE
562     !---------------------------------------------------------
563     !      DEALLOCATE (TWET)
564     
565           RETURN
566           END
567     !
568     !
569     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
570     !
571     ! DoPhase is a subroutine written and provided by Jim Ramer at NOAA/FSL
572     !
573     !    Ramer, J, 1993: An empirical technique for diagnosing precipitation
574     !           type from model output.  Preprints, 5th Conf. on Aviation
575     !           Weather Systems, Vienna, VA, Amer. Meteor. Soc., 227-230.
576     !
577     !   CODE ADAPTED FOR WRF POST  24 AUGUST 2005    G MANIKIN
578     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
579     !
580           SUBROUTINE CALWXT_RAMER(lm,lp1,      &
581                               T,Q,PMID,RH,TD,PINT,PREC,PTHRESH,PTYP)
582     
583     !      SUBROUTINE dophase(pq,   !  input pressure sounding mb
584     !     +    t,   !  input temperature sounding K
585     !     +    pmid,   !  input pressure
586     !     +    pint,   !  input interface pressure
587     !     +    q,   !  input spec humidityfraction
588     !     +    lmh,   !  input number of levels in sounding
589     !     +    prec,      ! input amount of precipitation
590     !     +    ptyp) !  output(2) phase 2=Rain, 3=Frzg, 4=Solid,
591     !                                               6=IP     JC  9/16/99
592     !      use params_mod
593     !      use CTLBLK_mod 
594     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
595           implicit none
596     !
597           real,PARAMETER :: LECP=1572.5
598           real,PARAMETER :: twice=266.55,rhprcp=0.80,deltag=1.02,prcpmin=0.3, &
599          &             emelt=0.045,rlim=0.04,slim=0.85
600           real,PARAMETER :: twmelt=273.15,tz=273.15,efac=1.0 ! specify in params now 
601     !
602           INTEGER*4 i, k1, lll, k2, toodry
603     !
604           REAL xxx ,mye, icefrac,flg
605           integer,intent(in) :: lm,lp1
606           real,DIMENSION(LM),intent(in) :: T,Q,PMID,RH,TD
607           real,DIMENSION(LP1),intent(in) :: PINT
608           real,intent(in) :: PREC,PTHRESH
609           real,intent(out) :: PTYP
610     !
611           real,DIMENSION(LM) :: P,TQ,QQ,PQ,RHQ
612           real,DIMENSION(LM) :: tqtmp,pqtmp,rhqtmp
613           real,DIMENSION(LM) :: TWQ
614     !
615           integer J,L,LEV,LNQ,ii
616           real RHMAX,TWMAX,PTOP,dpdrh,twtop,rhtop,wgt1,wgt2,    &
617                rhavg,dtavg,dpk,ptw,rate,pbot,qc, b,qtmp
618           real,external :: xmytw
619     !
620     !  Initialize.
621           icefrac = -9999.
622     !
623     
624           PTYP = 0
625           DO 88 L = 1,LM
626             LEV=LM-L+1
627     !        P(L)=PMID(L)
628     !        QC=PQ0/P(L) * EXP(A2*(T(L)-A3)/(T(L)-A4))
629     !GSM forcing Q (QTMP) to be positive to deal with negative Q values
630     !       causing problems later in this subroutine
631     !        QTMP=AMAX1(H1M12,Q(L))	
632     !        RHQTMP(LEV)=QTMP/QC
633     	RHQ(LEV)=RH(L)
634             PQ(LEV)=PMID(L)/100.
635             TQ(LEV)=T(L)
636        88 CONTINUE
637     
638     !      do 92 L=1,LM 
639     !         TQ(L)=TQTMP(L)
640     !         PQ(L)=PQTMP(L)
641     !         RHQ(L)=RHQTMP(L)
642     !   92 continue
643     
644     !  BIG LOOP
645     !      DO 800 J=JSTA,JEND
646     !      DO 800 I=1,IM
647     !
648     !   SKIP THIS POINT IF NO PRECIP THIS TIME STEP
649     !
650           IF (PREC.LE.PTHRESH) GOTO 800
651     
652     !
653     !
654     !CC   RATE RESTRICTION REMOVED BY JOHN CORTINAS 3/16/99
655     !
656     !     Construct wet-bulb sounding, locate generating level.
657           twmax = -999.0
658           rhmax = 0.0
659           k1 = 0    !  top of precip generating layer
660           k2 = 0    !  layer of maximum rh
661     !
662           IF (rhq(1).lt.rhprcp) THEN
663               toodry = 1
664           ELSE
665               toodry = 0
666           END IF
667     !
668           pbot = pq(1)
669     !      NQ=LM
670           DO 10 L = 1, lm
671     !          xxx = tdofesat(esat(tq(L))*rhq(L))
672     	  xxx = td(l) !HC: use TD consistent with GFS ice physics
673     	  if (xxx .lt. -500.) goto 800
674               twq(L) = xmytw(tq(L),xxx,pq(L))
675               twmax = amax1(twq(L),twmax)
676               IF (pq(L).ge.400.0) THEN
677                   IF (rhq(L).gt.rhmax) THEN
678                       rhmax = rhq(L)
679                       k2 = i
680                   END IF
681     !
682                   IF (L.ne.1) THEN
683                      IF (rhq(L).ge.rhprcp.or.toodry.eq.0) THEN
684                       IF (toodry.ne.0) THEN
685                         dpdrh = alog(pq(L)/pq(L-1)) /              &
686                                (rhq(L)-RHQ(L-1))
687                         pbot = exp(alog(pq(L))+(rhprcp-rhq(L))*dpdrh)
688     !
689                         ptw = pq(L)
690                         toodry = 0
691                         ELSE IF (rhq(L).ge.rhprcp) THEN
692                           ptw = pq(L)
693                         ELSE
694                           toodry = 1
695                           dpdrh = alog(pq(L)/pq(L-1)) /                 &
696                               (rhq(L)-rhq(L-1))
697                           ptw = exp(alog(pq(L))+(rhprcp-rhq(L))*dpdrh)
698     !lin                dpdrh=(Pq(i)-Pq(i-1))/(Rhq(i)-Rhq(i-1))
699     !lin                ptw=Pq(i)+(rhprcp-Rhq(i))*dpdrh
700     !
701                           END IF
702     !
703                           IF (pbot/ptw.ge.deltag) THEN
704     !lin                      If (pbot-ptw.lt.deltag) Goto 2003
705                               k1 = L
706                               ptop = ptw
707                           END IF
708                       END IF
709                   END IF
710               END IF
711     !
712        10 CONTINUE
713     
714     !
715     !     Gross checks for liquid and solid precip which dont require generating level.
716     !
717           IF (twq(1).ge.273.15+2.0) THEN
718               ptyp = 8   ! liquid
719               icefrac = 0.0
720               goto 800 
721           END IF
722     !
723           IF (twmax.le.twice) THEN
724               icefrac = 1.0
725               ptyp = 1   !  solid
726               goto 800 
727           END IF
728     !
729     !     Check to see if we had no success with locating a generating level.
730     !
731           IF (k1.eq.0) goto 800
732     !
733           IF (ptop.eq.pq(k1)) THEN
734               twtop = twq(k1)
735               rhtop = rhq(k1)
736               k2 = k1
737               k1 = k1 - 1
738           ELSE
739               k2 = k1
740               k1 = k1 - 1
741               wgt1 = alog(ptop/pq(k2)) / alog(pq(k1)/pq(k2))
742               wgt2 = 1.0 - wgt1
743               twtop = twq(k1) * wgt1 + twq(k2) * wgt2
744               rhtop = rhq(k1) * wgt1 + rhq(k2) * wgt2
745           END IF
746     !
747     
748     !     Calculate temp and wet-bulb ranges below precip generating level.
749           DO 20 L = 1, k1
750               twmax = amax1(twq(l),twmax)
751        20 CONTINUE
752     !
753     !     Gross check for solid precip, initialize ice fraction.
754           IF (i.eq.1.and.j.eq.1) WRITE (*,*) 'twmax=',twmax,twice,'twtop=',twtop
755           IF (twtop.le.twice) THEN
756               icefrac = 1.0
757               IF (twmax.le.twmelt) THEN     ! gross check for solid precip.
758                   ptyp = 1       !   solid precip
759                   goto 800 
760               END IF
761               lll = 0
762           ELSE
763               icefrac = 0.0
764               lll = 1
765           END IF
766     !
767     !     Loop downward through sounding from highest precip generating level.
768        30 CONTINUE
769     !
770           IF (icefrac.ge.1.0) THEN  !  starting as all ice
771               IF (twq(k1).lt.twmelt) GO TO 40       ! cannot commence melting
772               IF (twq(k1).eq.twtop) GO TO 40        ! both equal twmelt, nothing h
773               wgt1 = (twmelt-twq(k1)) / (twtop-twq(k1))
774               rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
775               dtavg = (twmelt-twq(k1)) / 2
776               dpk = wgt1 * alog(pq(k1)/ptop)        !lin   dpk=wgt1*(Pq(k1)-Ptop)
777     !         mye=emelt*(1.0-(1.0-Rhavg)*efac)
778               mye = emelt * rhavg ** efac
779               icefrac = icefrac + dpk * dtavg / mye
780           ELSE IF (icefrac.le.0.0) THEN     !  starting as all liquid
781               lll = 1
782     !         Goto 1020
783               IF (twq(k1).gt.twice) GO TO 40        ! cannot commence freezing
784               IF (twq(k1).eq.twtop) THEN
785                   wgt1 = 0.5
786               ELSE
787                   wgt1 = (twice-twq(k1)) / (twtop-twq(k1))
788               END IF
789               rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
790               dtavg = twmelt - (twq(k1)+twice) / 2
791               dpk = wgt1 * alog(pq(k1)/ptop)      !lin  dpk=wgt1*(Pq(k1)-Ptop)
792     !         mye=emelt*(1.0-(1.0-Rhavg)*efac)
793               mye = emelt * rhavg ** efac
794               icefrac = icefrac + dpk * dtavg / mye
795           ELSE IF ((twq(k1).le.twmelt).and.(twq(k1).lt.twmelt)) THEN ! mix
796               rhavg = (rhq(k1)+rhtop) / 2
797               dtavg = twmelt - (twq(k1)+twtop) / 2
798               dpk = alog(pq(k1)/ptop)       !lin   dpk=Pq(k1)-Ptop
799     !         mye=emelt*(1.0-(1.0-Rhavg)*efac)
800               mye = emelt * rhavg ** efac
801               icefrac = icefrac + dpk * dtavg / mye           
802           ELSE      ! mix where Tw curve crosses twmelt in layer
803               IF (twq(k1).eq.twtop) GO TO 40   ! both equal twmelt, nothing h
804               wgt1 = (twmelt-twq(k1)) / (twtop-twq(k1))
805               wgt2 = 1.0 - wgt1
806               rhavg = rhtop + wgt2 * (rhq(k1)-rhtop) / 2
807               dtavg = (twmelt-twtop) / 2
808               dpk = wgt2 * alog(pq(k1)/ptop)     !lin   dpk=wgt2*(Pq(k1)-Ptop)
809     !         mye=emelt*(1.0-(1.0-Rhavg)*efac)
810               mye = emelt * rhavg ** efac
811               icefrac = icefrac + dpk * dtavg / mye
812               icefrac = amin1(1.0,amax1(icefrac,0.0))   
813               IF (icefrac.le.0.0) THEN
814     !             Goto 1020
815                   IF (twq(k1).gt.twice) GO TO 40    ! cannot commence freezin
816                   wgt1 = (twice-twq(k1)) / (twtop-twq(k1))
817                   dtavg = twmelt - (twq(k1)+twice) / 2
818               ELSE
819                   dtavg = (twmelt-twq(k1)) / 2
820               END IF
821               rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
822               dpk = wgt1 * alog(pq(k1)/ptop)     !lin  dpk=wgt1*(Pq(k1)-Ptop)
823     !         mye=emelt*(1.0-(1.0-Rhavg)*efac)
824               mye = emelt * rhavg ** efac
825               icefrac = icefrac + dpk * dtavg / mye
826           END IF
827     !
828           icefrac = amin1(1.0,amax1(icefrac,0.0))
829           IF (i.eq.1.and.j.eq.1) WRITE (*,*) 'NEW ICEFRAC:', icefrac, icefrac
830     !
831     !     Get next level down if there is one, loop back.
832        40 IF (k1.gt.1) THEN
833               twtop = twq(k1)
834               ptop = pq(k1)
835               rhtop = rhq(k1)
836               k1 = k1 - 1
837               GO TO 30
838           END IF
839     !
840     !
841     !     Determine precip type based on snow fraction and surface wet-bulb.
842     !
843     !
844           IF (icefrac.ge.slim) THEN
845               IF (lll.ne.0) THEN
846                   ptyp = 2       ! Ice Pellets   JC 9/16/99
847               ELSE
848                   ptyp = 1       !  Snow
849               END IF
850           ELSE IF (icefrac.le.rlim) THEN
851               IF (twq(1).lt.tz) THEN
852                   ptyp = 4       !  Freezing Precip
853               ELSE
854                   ptyp = 8       !  Rain
855               END IF
856           ELSE
857               IF (twq(1).lt.tz) THEN
858     !GSM not sure what to do when 'mix' is predicted;   In previous
859     !GSM   versions of this code for which I had to have an answer,
860     !GSM   I chose sleet.  Here, though, since we have 4 other
861     !GSM   algorithms to provide an answer, I will not declare a
862     !GSM   type from the Ramer in this situation and allow the
863     !GSM   other algorithms to make the call.
864           
865                   ptyp = 0       !  don't know 
866     !              ptyp = 5       !  Mix
867               ELSE
868     !              ptyp = 5       !  Mix
869                   ptyp = 0       !  don't know 
870               END IF
871           END IF
872      800  CONTINUE 
873     
874           RETURN
875     !
876           END
877     !
878     !
879     !--------------------------------------------------------------------------
880     !      REAL*4 FUNCTION mytw(t,td,p)
881           FUNCTION xmytw(t,td,p)
882     !
883           IMPLICIT NONE
884     !
885           INTEGER*4 cflag, l
886     !     REAL*4 f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s,        &
887           REAL   f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s,        &
888          &          de, xmytw
889           DATA f, c0, c1, c2 /0.0006355, 26.66082, 0.0091379024, 6106.3960/
890     !
891     !
892           xmytw = (t+td) / 2
893           IF (td.ge.t) RETURN
894     !
895           IF (t.lt.100.0) THEN
896               k = t + 273.15
897               kd = td + 273.15
898               IF (kd.ge.k) RETURN
899               cflag = 1
900           ELSE
901               k = t
902               kd = td
903               cflag = 0
904           END IF
905     !
906           ed = c0 - c1 * kd - c2 / kd
907           IF (ed.lt.-14.0.or.ed.gt.7.0) RETURN
908           ed = exp(ed)
909           ew = c0 - c1 * k - c2 / k
910           IF (ew.lt.-14.0.or.ew.gt.7.0) RETURN
911           ew = exp(ew)
912           fp = p * f
913           s = (ew-ed) / (k-kd)
914           kw = (k*fp+kd*s) / (fp+s)
915     !
916           DO 10 l = 1, 5
917               ew = c0 - c1 * kw - c2 / kw
918               IF (ew.lt.-14.0.or.ew.gt.7.0) RETURN
919               ew = exp(ew)
920               de = fp * (k-kw) + ed - ew
921               IF (abs(de/ew).lt.1E-5) GO TO 20
922               s = ew * (c1-c2/(kw*kw)) - fp
923               kw = kw - de / s
924        10 CONTINUE
925        20 CONTINUE
926     !
927     !      print *, 'kw ', kw
928           IF (cflag.ne.0) THEN
929               xmytw = kw - 273.15
930           ELSE
931               xmytw = kw
932           END IF
933     !
934           RETURN
935           END
936     !
937     !
938     !$$$  Subprogram documentation block
939     !
940     ! Subprogram: calwxt_bourg    Calculate precipitation type (Bourgouin)
941     !   Prgmmr: Baldwin      Org: np22        Date: 1999-07-06
942     !
943     ! Abstract: This routine computes precipitation type
944     !    using a decision tree approach that uses the so-called
945     !    "energy method" of Bourgouin of AES (Canada) 1992
946     !
947     ! Program history log:
948     !   1999-07-06  M Baldwin
949     !   1999-09-20  M Baldwin  make more consistent with bourgouin (1992)
950     !   2005-08-24  G Manikin  added to wrf post
951     !   2007-06-19  M Iredell  mersenne twister, best practices
952     !   2008-03-03  G Manikin  added checks to prevent stratospheric warming
953     !                           episodes from being seen as "warm" layers
954     !                           impacting precip type
955     !
956     ! Usage:    call calwxt_bourg(im,jm,jsta_2l,jend_2u,jsta,jend,lm,lp1,   &
957     !    &                        iseed,g,pthresh,                          &
958     !    &                        t,q,pmid,pint,lmh,prec,zint,ptype)
959     !   Input argument list:
960     !     im       integer i dimension
961     !     jm       integer j dimension
962     !     jsta_2l  integer j dimension start point (including haloes)
963     !     jend_2u  integer j dimension end point (including haloes)
964     !     jsta     integer j dimension start point (excluding haloes)
965     !     jend     integer j dimension end point (excluding haloes)
966     !     lm       integer k dimension
967     !     lp1      integer k dimension plus 1
968     !     iseed    integer random number seed
969     !     g        real gravity (m/s**2)
970     !     pthresh  real precipitation threshold (m)
971     !     t        real(im,jsta_2l:jend_2u,lm) mid layer temp (K)
972     !     q        real(im,jsta_2l:jend_2u,lm) specific humidity (kg/kg)
973     !     pmid     real(im,jsta_2l:jend_2u,lm) mid layer pressure (Pa)
974     !     pint     real(im,jsta_2l:jend_2u,lp1) interface pressure (Pa)
975     !     lmh      real(im,jsta_2l:jend_2u) max number of layers
976     !     prec     real(im,jsta_2l:jend_2u) precipitation (m)
977     !     zint     real(im,jsta_2l:jend_2u,lp1) interface height (m)
978     !   Output argument list:
979     !     ptype    real(im,jm) instantaneous weather type ()
980     !              acts like a 4 bit binary
981     !                1111 = rain/freezing rain/ice pellets/snow
982     !                where the one's digit is for snow
983     !                      the two's digit is for ice pellets
984     !                      the four's digit is for freezing rain
985     !                  and the eight's digit is for rain
986     !              in other words...
987     !                ptype=1 snow
988     !                ptype=2 ice pellets/mix with ice pellets
989     !                ptype=4 freezing rain/mix with freezing rain
990     !                ptype=8 rain
991     !
992     ! Modules used:
993     !   mersenne_twister pseudo-random number generator
994     !
995     ! Subprograms called:
996     !   random_number    pseudo-random number generator
997     !
998     ! Attributes:
999     !   Language: Fortran 90
1000     !
1001     ! Remarks: vertical order of arrays must be layer   1 = top
1002     !                                       and layer lmh = bottom
1003     !
1004     !$$$
1005           subroutine calwxt_bourg(lm,lp1,rn,g,pthresh,      &
1006          &                        t,q,pmid,pint,prec,zint,ptype)
1007     !      use mersenne_twister
1008           implicit none
1009     !
1010     !    input:
1011           integer,intent(in):: lm,lp1
1012     !      integer,intent(in):: iseed
1013           real,intent(in):: g,pthresh,rn(2)
1014           real,intent(in):: t(lm)
1015           real,intent(in):: q(lm)
1016           real,intent(in):: pmid(lm)
1017           real,intent(in):: pint(lp1)
1018           real,intent(in):: prec
1019           real,intent(in):: zint(lp1)
1020     !
1021     !    output:
1022           real,intent(out):: ptype
1023     !
1024           integer ifrzl,iwrml,l,lhiwrm,lmhk
1025           real pintk1,areane,tlmhk,areape,pintk2,surfw,area1,dzkl,psfck
1026           real r1,r2
1027     !
1028     !     initialize weather type array to zero (ie, off).
1029     !     we do this since we want ptype to represent the
1030     !     instantaneous weather type on return.
1031     !     
1032     !!$omp  parallel do
1033     
1034           ptype = 0
1035     
1036     !
1037     !      call random_number(rn,iseed)
1038     !
1039     !!$omp  parallel do
1040     !!$omp& private(a,lmhk,tlmhk,iwrml,psfck,lhiwrm,pintk1,pintk2,area1,
1041     !!$omp&         areape,dzkl,surfw,r1,r2)
1042     
1043           psfck=pint(lm+1)
1044     !
1045     !   skip this point if no precip this time step 
1046     !
1047           if (prec.le.pthresh) return
1048     !     find the depth of the warm layer based at the surface
1049     !     this will be the cut off point between computing
1050     !     the surface based warm air and the warm air aloft
1051     !
1052     !
1053     !     lowest layer t
1054     !
1055           tlmhk = t(lm)
1056           iwrml = lm + 1
1057           if (tlmhk.ge.273.15) then
1058             do l = lm, 2, -1
1059              if (t(l).ge.273.15.and.t(l-1).lt.273.15.and.           &
1060          &            iwrml.eq.lm+1) iwrml = l
1061               end do
1062           end if
1063     !
1064     !     now find the highest above freezing level
1065     !
1066           lhiwrm = lm + 1
1067           do l = lm, 1, -1
1068     ! gsm  added 250 mb check to prevent stratospheric warming situations
1069     !       from counting as warm layers aloft      
1070               if (t(l).ge.273.15 .and. pmid(l).gt.25000.) lhiwrm = l
1071           end do
1072     
1073     !     energy variables
1074     !     surfw is the positive energy between the ground
1075     !     and the first sub-freezing layer above ground
1076     !     areane is the negative energy between the ground
1077     !     and the highest layer above ground
1078     !     that is above freezing
1079     !     areape is the positive energy "aloft"
1080     !     which is the warm energy not based at the ground
1081     !     (the total warm energy = surfw + areape)
1082     !
1083     !     pintk1 is the pressure at the bottom of the layer
1084     !     pintk2 is the pressure at the top of the layer
1085     !     dzkl is the thickness of the layer
1086     !     ifrzl is a flag that tells us if we have hit
1087     !     a below freezing layer
1088     !
1089           pintk1 = psfck
1090           ifrzl = 0
1091           areane = 0.0
1092           areape = 0.0
1093           surfw = 0.0                                         
1094     
1095           do l = lm, 1, -1
1096               if (ifrzl.eq.0.and.t(l).le.273.15) ifrzl = 1
1097               pintk2=pint(l)
1098               dzkl=zint(l)-zint(l+1)
1099               area1 = alog(t(l)/273.15) * g * dzkl
1100               if (t(l).ge.273.15.and. pmid(l).gt.25000.) then
1101                   if (l.lt.iwrml) areape = areape + area1
1102                   if (l.ge.iwrml) surfw = surfw + area1
1103               else
1104                   if (l.gt.lhiwrm) areane = areane + abs(area1)
1105               end if
1106               pintk1 = pintk2
1107           end do
1108           
1109     !
1110     !     decision tree time
1111     !
1112           if (areape.lt.2.0) then
1113     !         very little or no positive energy aloft, check for
1114     !         positive energy just above the surface to determine rain vs. snow
1115               if (surfw.lt.5.6) then
1116     !             not enough positive energy just above the surface
1117     !             snow = 1
1118                   ptype = 1
1119               else if (surfw.gt.13.2) then
1120     !             enough positive energy just above the surface
1121     !             rain = 8
1122                   ptype = 8
1123               else
1124     !             transition zone, assume equally likely rain/snow
1125     !             picking a random number, if <=0.5 snow
1126                   r1 = rn(1)
1127                   if (r1.le.0.5) then
1128     !                 snow = 1
1129                       ptype = 1
1130                   else
1131     !                 rain = 8
1132                       ptype = 8
1133                   end if
1134               end if
1135     !
1136           else
1137     !         some positive energy aloft, check for enough negative energy
1138     !         to freeze and make ice pellets to determine ip vs. zr
1139               if (areane.gt.66.0+0.66*areape) then
1140     !             enough negative area to make ip,
1141     !             now need to check if there is enough positive energy
1142     !             just above the surface to melt ip to make rain
1143                   if (surfw.lt.5.6) then
1144     !                 not enough energy at the surface to melt ip
1145     !                 ice pellets = 2
1146                       ptype = 2
1147                   else if (surfw.gt.13.2) then
1148     !                 enough energy at the surface to melt ip
1149     !                 rain = 8
1150                       ptype = 8
1151                   else
1152     !                 transition zone, assume equally likely ip/rain
1153     !                 picking a random number, if <=0.5 ip
1154                       r1 = rn(1)
1155                       if (r1.le.0.5) then
1156     !                     ice pellets = 2
1157                           ptype = 2
1158                       else
1159     !                     rain = 8
1160                           ptype = 8
1161                       end if
1162                   end if
1163               else if (areane.lt.46.0+0.66*areape) then
1164     !             not enough negative energy to refreeze, check surface temp
1165     !             to determine rain vs. zr
1166                   if (tlmhk.lt.273.15) then
1167     !                 freezing rain = 4
1168                       ptype = 4
1169                   else
1170     !                 rain = 8
1171                       ptype = 8
1172                   end if
1173               else
1174     !             transition zone, assume equally likely ip/zr
1175     !             picking a random number, if <=0.5 ip
1176                   r1 = rn(1)
1177                   if (r1.le.0.5) then
1178     !                 still need to check positive energy
1179     !                 just above the surface to melt ip vs. rain
1180                       if (surfw.lt.5.6) then
1181     !                     ice pellets = 2
1182                           ptype = 2
1183                       else if (surfw.gt.13.2) then
1184     !                     rain = 8
1185                           ptype = 8
1186                       else
1187     !                     transition zone, assume equally likely ip/rain
1188     !                     picking a random number, if <=0.5 ip
1189                           r2 = rn(2)
1190                           if (r2.le.0.5) then
1191     !                         ice pellets = 2
1192                               ptype = 2
1193                           else
1194     !                         rain = 8
1195                               ptype = 8
1196                           end if
1197                       end if
1198                   else
1199     !                 not enough negative energy to refreeze, check surface temp
1200     !                 to determine rain vs. zr
1201                       if (tlmhk.lt.273.15) then
1202     !                     freezing rain = 4
1203                           ptype = 4
1204                       else
1205     !                     rain = 8
1206                           ptype = 8
1207                       end if
1208                   end if
1209               end if
1210           end if
1211     !      end do
1212     !      end do
1213           return
1214           end
1215     !
1216     !
1217            SUBROUTINE CALWXT_REVISED(LM,LP1,T,Q,PMID,PINT,PREC,  &
1218                      PTHRESH,D608,ROG,EPSQ,    &
1219          &             ZINT,TWET,IWX)
1220     ! 
1221     !     FILE: CALWXT.f
1222     !     WRITTEN: 11 NOVEMBER 1993, MICHAEL BALDWIN
1223     !     REVISIONS:
1224     !               30 SEPT 1994-SETUP NEW DECISION TREE (M BALDWIN)
1225     !               12 JUNE 1998-CONVERSION TO 2-D (T BLACK)
1226     !     01-10-25  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
1227     !     02-01-15  MIKE BALDWIN - WRF VERSION
1228     !     05-07-07  BINBIN ZHOU  - ADD PREC FOR RSM
1229     !     05-08-24  GEOFF MANIKIN - MODIFIED THE AREA REQUIREMENTS
1230     !                TO MAKE AN ALTERNATE ALGORITHM 
1231     !                              
1232     !
1233     !     ROUTINE TO COMPUTE PRECIPITATION TYPE USING A DECISION TREE
1234     !     APPROACH THAT USES VARIABLES SUCH AS INTEGRATED WET BULB TEMP
1235     !     BELOW FREEZING AND LOWEST LAYER TEMPERATURE
1236     !
1237     !     SEE BALDWIN AND CONTORNO PREPRINT FROM 13TH WEATHER ANALYSIS
1238     !     AND FORECASTING CONFERENCE FOR MORE DETAILS
1239     !     (OR BALDWIN ET AL, 10TH NWP CONFERENCE PREPRINT)
1240     !
1241     !     SINCE THE ORIGINAL VERSION OF THE ALGORITHM HAS A HIGH BIAS
1242     !      FOR FREEZING RAIN AND SLEET, THE GOAL IS TO BALANCE THAT BIAS
1243     !      WITH A VERSION MORE LIKELY TO PREDICT SNOW
1244     !
1245     !     use params_mod
1246     !     use ctlblk_mod
1247     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1248           implicit none
1249     !
1250     !  LIST OF VARIABLES NEEDED
1251     !    PARAMETERS:
1252     !      D608,ROG,H1,D00
1253     !HC       PARAMETER(D608=0.608,ROG=287.04/9.8,H1=1.0,D00=0.0)
1254     !
1255     !    INPUT:
1256     !      T,Q,PMID,HTM,LMH,PREC,ZINT
1257           integer,intent(in):: lm,lp1
1258           REAL,dimension(LM),intent(in) ::  T,Q,PMID,TWET
1259           REAL,dimension(LP1),intent(in) ::  PINT,ZINT 
1260           REAL,intent(in) ::  PREC,PTHRESH,D608,ROG,EPSQ
1261     !    OUTPUT:
1262     !      IWX - INSTANTANEOUS WEATHER TYPE.
1263     !        ACTS LIKE A 4 BIT BINARY
1264     !          1111 = RAIN/FREEZING RAIN/ICE PELLETS/SNOW
1265     !          WHERE THE ONE'S DIGIT IS FOR SNOW
1266     !                THE TWO'S DIGIT IS FOR ICE PELLETS
1267     !                THE FOUR'S DIGIT IS FOR FREEZING RAIN
1268     !            AND THE EIGHT'S DIGIT IS FOR RAIN
1269           integer, intent(out) ::  IWX
1270     !    INTERNAL:
1271     !
1272           real, parameter :: D00=0.0  
1273           integer KARR,LICEE
1274           real TCOLD,TWARM
1275     !
1276           integer I,J,L,LMHK,LICE,IFREL,IWRML,IFRZL
1277           real PSFCK,TDCHK,A,TDKL,TDPRE,TLMHK,TWRMK,AREAS8,AREAP4,AREA1,  &
1278                SURFW,SURFC,DZKL,PINTK1,PINTK2,PM150,QKL,TKL,PKL,AREA0,    &
1279                AREAP0
1280     
1281     !    SUBROUTINES CALLED:
1282     !     WETBULB
1283     !     
1284     !
1285     !     INITIALIZE WEATHER TYPE ARRAY TO ZERO (IE, OFF).
1286     !     WE DO THIS SINCE WE WANT IWX TO REPRESENT THE
1287     !     INSTANTANEOUS WEATHER TYPE ON RETURN.
1288     !     
1289     !
1290     !     ALLOCATE LOCAL STORAGE
1291     !
1292     !
1293     !!$omp  parallel do
1294           IWX = 0
1295     
1296     !!$omp  parallel do
1297     !!$omp& private(a,lmhk,pkl,psfck,qkl,tdchk,tdkl,tdpre,tkl)
1298     
1299           LMHK=LM
1300     !
1301     !   SKIP THIS POINT IF NO PRECIP THIS TIME STEP 
1302     !
1303           IF (PREC.LE.PTHRESH) GOTO 800
1304     !
1305     !   FIND COLDEST AND WARMEST TEMPS IN SATURATED LAYER BETWEEN
1306     !   70 MB ABOVE GROUND AND 500 MB
1307     !   ALSO FIND HIGHEST SATURATED LAYER IN THAT RANGE
1308     !
1309     !meb
1310           PSFCK=PINT(LP1)
1311     !meb
1312           TDCHK=2.0
1313       760 TCOLD=T(LMHK)
1314           TWARM=T(LMHK)
1315           LICEE=LMHK
1316     !
1317           DO 775 L=1,LMHK
1318           QKL=Q(L)
1319           QKL=AMAX1(EPSQ,QKL)
1320           TKL=T(L)
1321           PKL=PMID(L)
1322     !
1323     !   SKIP PAST THIS IF THE LAYER IS NOT BETWEEN 70 MB ABOVE GROUND
1324     !       AND 500 MB
1325     !
1326           IF (PKL.LT.50000.0.OR.PKL.GT.PSFCK-7000.0) GOTO 775
1327           A=ALOG(QKL*PKL/(6.1078*(0.378*QKL+0.622)))
1328           TDKL=(237.3*A)/(17.269-A)+273.15
1329           TDPRE=TKL-TDKL
1330           IF (TDPRE.LT.TDCHK.AND.TKL.LT.TCOLD) TCOLD=TKL
1331           IF (TDPRE.LT.TDCHK.AND.TKL.GT.TWARM) TWARM=TKL
1332           IF (TDPRE.LT.TDCHK.AND.L.LT.LICEE) LICEE=L
1333       775 CONTINUE
1334     !
1335     !    IF NO SAT LAYER AT DEW POINT DEP=TDCHK, INCREASE TDCHK
1336     !     AND START AGAIN (BUT DON'T MAKE TDCHK > 6)
1337     !
1338           IF (TCOLD.EQ.T(LMHK).AND.TDCHK.LT.6.0) THEN
1339             TDCHK=TDCHK+2.0
1340             GOTO 760
1341           ENDIF
1342       800 CONTINUE
1343     !
1344     !    LOWEST LAYER T
1345     !
1346           KARR=0
1347           IF (PREC.LE.PTHRESH) GOTO 850
1348           LMHK=LM
1349           TLMHK=T(LMHK)
1350     !
1351     !    DECISION TREE TIME
1352     !
1353           IF (TCOLD.GT.269.15) THEN
1354               IF (TLMHK.LE.273.15) THEN
1355     !             TURN ON THE FLAG FOR
1356     !             FREEZING RAIN = 4
1357     !             IF ITS NOT ON ALREADY
1358     !             IZR=MOD(IWX,8)/4
1359     !             IF (IZR.LT.1) IWX=IWX+4
1360                   IWX=IWX+4
1361                 GOTO 850
1362               ELSE
1363     !             TURN ON THE FLAG FOR
1364     !             RAIN = 8
1365     !             IF ITS NOT ON ALREADY
1366     !             IRAIN=IWX/8
1367     !             IF (IRAIN.LT.1) IWX=IWX+8
1368                   IWX=IWX+8
1369                 GOTO 850
1370               ENDIF
1371           ENDIF
1372           KARR=1
1373       850 CONTINUE
1374     !
1375     !!$omp  parallel do
1376     !!$omp& private(area1,areap4,areap0,areas8,dzkl,ifrzl,iwrml,lice,
1377     !!$omp&         lmhk,pintk1,pintk2,pm150,psfck,surfc,surfw,
1378     !!$omp&         tlmhk,twrmk)
1379     
1380           IF(KARR.GT.0)THEN
1381             LMHK=LM
1382             LICE=LICEE
1383     !meb
1384             PSFCK=PINT(LP1)
1385     !meb
1386             TLMHK=T(LMHK)
1387             TWRMK=TWARM
1388     !
1389     !    TWET AREA VARIABLES
1390     !     CALCULATE ONLY WHAT IS NEEDED
1391     !      FROM GROUND TO 150 MB ABOVE SURFACE
1392     !      FROM GROUND TO TCOLD LAYER
1393     !      AND FROM GROUND TO 1ST LAYER WHERE WET BULB T < 0.0
1394     !
1395     !     PINTK1 IS THE PRESSURE AT THE BOTTOM OF THE LAYER
1396     !     PINTK2 IS THE PRESSURE AT THE TOP OF THE LAYER
1397     !
1398     !     AREAP4 IS THE AREA OF TWET ABOVE -4 C BELOW HIGHEST SAT LYR 
1399     !     AREAP0 IS THE AREA OF TWET ABOVE 0 C BELOW HIGHEST SAT LYR
1400     !
1401             AREAS8=D00
1402             AREAP4=D00
1403     	AREAP0=D00
1404             SURFW =D00
1405             SURFC =D00
1406             
1407     !
1408             DO 1945 L=LMHK,LICE,-1
1409             DZKL=ZINT(L)-ZINT(L+1)
1410             AREA1=(TWET(L)-269.15)*DZKL
1411             AREA0=(TWET(L)-273.15)*DZKL
1412             IF (TWET(L).GE.269.15) AREAP4=AREAP4+AREA1
1413             IF (TWET(L).GE.273.15) AREAP0=AREAP0+AREA0
1414      1945   CONTINUE
1415     !
1416     !        IF (AREAP4.LT.3000.0) THEN
1417     !             TURN ON THE FLAG FOR
1418     !             SNOW = 1
1419     !             IF ITS NOT ON ALREADY
1420     !             ISNO=MOD(IWX,2)
1421     !             IF (ISNO.LT.1) IWX=IWX+1
1422     !          IWX=IWX+1
1423     !          GO TO 1900
1424     !        ENDIF
1425             IF (AREAP0.LT.350.0) THEN
1426     !             TURN ON THE FLAG FOR
1427     !             SNOW = 1
1428                   IWX=IWX+1
1429                 GOTO 1900
1430            ENDIF
1431     !
1432     !     AREAS8 IS THE NET AREA OF TWET W.R.T. FREEZING IN LOWEST 150MB
1433     !
1434             PINTK1=PSFCK
1435             PM150=PSFCK-15000.
1436     !
1437             DO 1955 L=LMHK,1,-1
1438             PINTK2=PINT(L)
1439             IF(PINTK1.LT.PM150)GO TO 1950
1440             DZKL=ZINT(L)-ZINT(L+1)
1441     !
1442     !    SUM PARTIAL LAYER IF IN 150 MB AGL LAYER
1443     !
1444             IF(PINTK2.LT.PM150)                                   &
1445               DZKL=T(L)*(Q(L)*D608+1.0)*ROG*               &
1446                    ALOG(PINTK1/PM150)
1447             AREA1=(TWET(L)-273.15)*DZKL
1448             AREAS8=AREAS8+AREA1
1449      1950   PINTK1=PINTK2
1450      1955   CONTINUE
1451     !
1452     !     SURFW IS THE AREA OF TWET ABOVE FREEZING BETWEEN THE GROUND
1453     !       AND THE FIRST LAYER ABOVE GROUND BELOW FREEZING
1454     !     SURFC IS THE AREA OF TWET BELOW FREEZING BETWEEN THE GROUND
1455     !       AND THE WARMEST SAT LAYER
1456     !
1457             IFRZL=0
1458             IWRML=0
1459     !
1460             DO 2050 L=LMHK,1,-1
1461             IF (IFRZL.EQ.0.AND.T(L).LT.273.15) IFRZL=1
1462             IF (IWRML.EQ.0.AND.T(L).GE.TWRMK) IWRML=1
1463     !
1464             IF (IWRML.EQ.0.OR.IFRZL.EQ.0) THEN
1465     !	  if(pmid(l) .lt. 50000.)print*,'twet needed above 500mb'
1466               DZKL=ZINT(L)-ZINT(L+1)
1467               AREA1=(TWET(L)-273.15)*DZKL
1468               IF(IFRZL.EQ.0.AND.TWET(L).GE.273.15)SURFW=SURFW+AREA1
1469               IF(IWRML.EQ.0.AND.TWET(L).LE.273.15)SURFC=SURFC+AREA1
1470             ENDIF
1471      2050   CONTINUE
1472             IF(SURFC.LT.-3000.0.OR.                                    &
1473          &    (AREAS8.LT.-3000.0.AND.SURFW.LT.50.0)) THEN
1474     !             TURN ON THE FLAG FOR
1475     !             ICE PELLETS = 2
1476     !             IF ITS NOT ON ALREADY
1477     !             IIP=MOD(IWX,4)/2
1478     !             IF (IIP.LT.1) IWX=IWX+2
1479               IWX=IWX+2
1480               GOTO 1900
1481             ENDIF
1482     !
1483             IF(TLMHK.LT.273.15) THEN
1484     !             TURN ON THE FLAG FOR
1485     !             FREEZING RAIN = 4
1486     !             IF ITS NOT ON ALREADY
1487     !             IZR=MOD(IWX(K),8)/4
1488     !             IF (IZR.LT.1) IWX(K)=IWX(K)+4
1489               IWX=IWX+4
1490             ELSE
1491     !             TURN ON THE FLAG FOR
1492     !             RAIN = 8
1493     !             IF ITS NOT ON ALREADY
1494     !             IRAIN=IWX(K)/8
1495     !             IF (IRAIN.LT.1) IWX(K)=IWX(K)+8
1496               IWX=IWX+8
1497             ENDIF
1498           ENDIF
1499      1900 CONTINUE
1500     !      print *, 'revised check ', IWX(500,800)
1501     
1502           RETURN
1503           END
1504     !
1505     !
1506           SUBROUTINE CALWXT_EXPLICIT(LM,PTHRESH,TSKIN,PREC,SR,F_RIMEF,IWX)
1507     ! 
1508     !     FILE: CALWXT.f
1509     !     WRITTEN: 24 AUGUST 2005, G MANIKIN and B FERRIER 
1510     !
1511     !     ROUTINE TO COMPUTE PRECIPITATION TYPE USING EXPLICIT FIELDS
1512     !       FROM THE MODEL MICROPHYSICS
1513     
1514     !      use params_mod
1515     !      use ctlblk_mod
1516     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1517           implicit none
1518     !
1519     !  LIST OF VARIABLES NEEDED
1520     !    PARAMETERS:
1521     !
1522     !    INPUT:
1523           integer, intent(in):: lm
1524           real,intent(in)::  TSKIN, PREC, SR,PTHRESH
1525           REAL,intent(in):: F_RimeF(LM)
1526           integer,intent(out) :: IWX
1527           integer I,J,LMHK
1528           real PSFC,SNOW
1529     !
1530     !     ALLOCATE LOCAL STORAGE
1531     !
1532     !!$omp  parallel do
1533           IWX = 0
1534     
1535     !GSM  THE RSM IS CURRENTLY INCOMPATIBLE WITH THIS ROUTINE
1536     !GSM   ACCORDING TO B FERRIER, THERE MAY BE A WAY TO WRITE
1537     !GSM   A VERSION OF THIS ALGORITHM TO WORK WITH THE RSM
1538     !GSM   MICROPHYSICS, BUT IT DOESN'T EXIST AT THIS TIME
1539     !!$omp  parallel do
1540     !!$omp& private(lmhk,psfc,tskin)
1541     
1542     !   SKIP THIS POINT IF NO PRECIP THIS TIME STEP 
1543     !
1544           IF (PREC.LE.PTHRESH) GOTO 800
1545     !
1546     !  A SNOW RATIO LESS THAN 0.5 ELIMINATES SNOW AND SLEET
1547     !   USE THE SKIN TEMPERATURE TO DISTINGUISH RAIN FROM FREEZING RAIN
1548     !   NOTE THAT 2-M TEMPERATURE MAY BE A BETTER CHOICE IF THE MODEL
1549     !   HAS A COLD BIAS FOR SKIN TEMPERATURE
1550     ! 
1551           IF (SR.LT.0.5) THEN
1552     !        SURFACE (SKIN) POTENTIAL TEMPERATURE AND TEMPERATURE.
1553     !         PSFC=PMID(LM)
1554     !         TSKIN=THS*(PSFC/P1000)**CAPA 
1555     
1556              IF (TSKIN.LT.273.15) THEN
1557     !          FREEZING RAIN = 4
1558                IWX=IWX+4
1559              ELSE
1560     !          RAIN = 8
1561                IWX=IWX+8
1562              ENDIF
1563           ELSE
1564     !  
1565     !  DISTINGUISH SNOW FROM SLEET WITH THE RIME FACTOR
1566     ! 
1567             IF(F_RimeF(LM).GE.10) THEN
1568     !          SLEET = 2
1569                IWX=IWX+2
1570             ELSE
1571                SNOW = 1
1572                IWX=IWX+1 
1573             ENDIF
1574           ENDIF
1575      800  CONTINUE
1576      810  RETURN 
1577           END
1578     !
1579     !
1580            SUBROUTINE CALWXT_DOMINANT(NALG,PREC,PTHRESH,RAIN,FREEZR,SLEET,SNOW,     &
1581          &         DOMR,DOMZR,DOMIP,DOMS)
1582     !
1583     !     WRITTEN: 24 AUGUST 2005, G MANIKIN 
1584     !      
1585     !     THIS ROUTINE TAKES THE PRECIP TYPE SOLUTIONS FROM DIFFERENT
1586     !       ALGORITHMS AND SUMS THEM UP TO GIVE A DOMINANT TYPE
1587     !
1588     !      use params_mod
1589     !      use ctlblk_mod
1590     !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1591           implicit none
1592     !
1593     !    INPUT:
1594           integer,intent(in) :: NALG
1595           REAL, intent(in) :: PREC,PTHRESH
1596           real,intent(out) ::  DOMS,DOMR,DOMZR,DOMIP
1597           real,DIMENSION(NALG),intent(in) ::  RAIN,SNOW,SLEET,FREEZR
1598           integer I,J,L
1599           real TOTSN,TOTIP,TOTR,TOTZR
1600     !--------------------------------------------------------------------------
1601     !      print* , 'into dominant'
1602     !!$omp  parallel do
1603           DOMR = 0.
1604           DOMS = 0.
1605           DOMZR = 0.
1606           DOMIP = 0.
1607     !
1608     !!$omp  parallel do
1609     !!$omp& private(totsn,totip,totr,totzr)
1610     !   SKIP THIS POINT IF NO PRECIP THIS TIME STEP
1611           IF (PREC.LE.PTHRESH) GOTO 800
1612           TOTSN = 0
1613           TOTIP = 0
1614           TOTR  = 0
1615           TOTZR = 0 
1616     !   LOOP OVER THE NUMBER OF DIFFERENT ALGORITHMS THAT ARE USED
1617           DO 820 L = 1, NALG
1618             IF (RAIN(L).GT. 0) THEN
1619                TOTR = TOTR + 1
1620                GOTO 830
1621             ENDIF
1622     
1623             IF (SNOW(L).GT. 0) THEN
1624                TOTSN = TOTSN + 1
1625                GOTO 830
1626             ENDIF
1627     
1628             IF (SLEET(L).GT. 0) THEN
1629                TOTIP = TOTIP + 1
1630                GOTO 830
1631             ENDIF
1632     
1633             IF (FREEZR(L).GT. 0) THEN
1634                TOTZR = TOTZR + 1
1635                GOTO 830
1636             ENDIF
1637      830    CONTINUE
1638      820  CONTINUE
1639     
1640     !   TIES ARE BROKEN TO FAVOR THE MOST DANGEROUS FORM OF PRECIP
1641     !     FREEZING RAIN > SNOW > SLEET > RAIN 
1642           IF (TOTSN .GT. TOTIP) THEN
1643             IF (TOTSN .GT. TOTZR) THEN
1644               IF (TOTSN .GE. TOTR) THEN
1645                DOMS = 1
1646                GOTO 800 
1647               ELSE
1648                DOMR = 1 
1649                GOTO 800 
1650               ENDIF
1651             ELSE IF (TOTZR .GE. TOTR) THEN
1652               DOMZR = 1
1653               GOTO 800 
1654             ELSE
1655               DOMR = 1
1656               GOTO 800 
1657             ENDIF 
1658           ELSE IF (TOTIP .GT. TOTZR) THEN
1659             IF (TOTIP .GE. TOTR) THEN
1660               DOMIP = 1
1661               GOTO 800 
1662             ELSE
1663               DOMR = 1
1664               GOTO 800 
1665             ENDIF
1666           ELSE IF (TOTZR .GE. TOTR) THEN
1667              DOMZR = 1
1668              GOTO 800 
1669           ELSE
1670               DOMR = 1
1671               GOTO 800 
1672           ENDIF
1673      800  CONTINUE 
1674           RETURN
1675           END
1676     
1677           
1678                 
1679           
1680           
1681