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

1            SUBROUTINE PRECPD (IM,IX,KM,DT,DEL,PRSL,PS,Q,CWM,T,RN
2          &,                   rainp,u00k,psautco,prautco,evpco,lprnt,jpr)
3     !
4     !
5     !     ******************************************************************
6     !     *                                                                *
7     !     *           SUBROUTINE FOR PRECIPITATION PROCESSES               *
8     !     *           FROM SUSPENDED CLOUD WATER/ICE                       *
9     !     *                                                                *
10     !     ******************************************************************
11     !     *                                                                *
12     !     *  Originally CREATED BY  Q. ZHAO                JAN. 1995       *
13     !     *                         -------                                *    
14     !     *  Modified and rewritten by Shrinivas Moorthi   Oct. 1998       *
15     !     *                            -----------------                   *
16     !     *  and                       Hua-Lu Pan                          *
17     !     *                            ----------                          *
18     !     *                                                                *
19     !     *  References:                                                   *
20     !     *                                                                *
21     !     *  Zhao and Carr (1997), Monthly Weather Review (August)         *
22     !     *  Sundqvist et al., (1989) Monthly Weather review. (August)     *
23     !     *                                                                *
24     !     ******************************************************************
25     !
26     !     In this code vertical indexing runs from surface to top of the
27     !     model
28     !
29     !     Argument List:
30     !     --------------
31     !       IM         : Inner dimension over which calculation is made
32     !       IX         : Maximum inner dimension
33     !       KM         : Number of vertical levels
34     !       DT         : Time step in seconds
35     !       DEL(KM)    : Pressure layer thickness (Bottom to top)
36     !       PRSL(KM)   : Pressure values for model layers (bottom to top)
37     !       PS(IM)     : Surface pressure (centibars)
38     !       Q(IX,KM)   : Specific Humidity (Updated in the code)
39     !       CWM(IX,KM) : Condensate mixing ratio (Updated in the code)
40     !       T(IX,KM)   : Temperature       (Updated in the code)
41     !       RN(IM)     : Precipitation over one time-step DT (m/DT)
42     !       SR(IM)     : Index (=-1 Snow, =0 Rain/Snow, =1 Rain)
43     !       TCW(IM)    : Vertically integrated liquid water (Kg/m**2)
44     !       CLL(IX,KM) : Cloud cover
45     !hchuang RN(IM) unit in m per time step
46     !        precipitation rate conversion 1 mm/s = 1 kg/m2/s
47     !
48           USE MACHINE , ONLY : kind_phys
49           USE FUNCPHYS , ONLY : fpvs
50           USE PHYSCONS, grav => con_g, HVAP => con_HVAP, HFUS => con_HFUS
51          &,             TTP => con_TTP, CP => con_CP
52          &,             EPS => con_eps, EPSM1 => con_epsm1
53           implicit none
54     !     include 'constant.h'
55     !
56           real (kind=kind_phys) G,      H1,    H2,   H1000
57          &,                     H1000G, D00,   D125, D5
58          &,                     ELWV,   ELIV,  ROW
59          &,                     EPSQ,   DLDT,  TM10, ELIW
60          &,                     RCP,    RROW
61            PARAMETER (G=grav,         H1=1.E0,     H2=2.E0,     H1000=1000.0
62          &,           H1000G=H1000/G, D00=0.E0,    D125=.125E0, D5=0.5E0
63          &,           ELWV=HVAP,      ELIV=HVAP+HFUS,   ROW=1.E3
64          &,           EPSQ=2.E-12,    DLDT=2274.E0,TM10=TTP-10.0
65          &,           ELIW=ELIV-ELWV, RCP=H1/CP,   RROW=H1/ROW)
66     !
67           real(kind=kind_phys), parameter :: cons_0=0.0,     cons_p01=0.01
68          &,                                  cons_20=20.0
69          &,                                  cons_m30=-30.0, cons_50=50.0
70     !
71           integer IM, IX, KM, LAT, jpr
72           real (kind=kind_phys) Q(IX,KM),   T(IX,KM),    CWM(IX,KM)
73          &,                                 DEL(IX,KM),  PRSL(IX,KM)
74     !    &,                     CLL(IM,KM), DEL(IX,KM),  PRSL(IX,KM)
75          &,                     PS(IM),     RN(IM),      SR(IM)
76          &,                     TCW(IM),    DT
77     !hchuang code change [+1L] : add record to record information in vertical in
78     !                       addition to total column PRECRL
79          &,                     RAINP(IM,KM), RNP(IM),
80          &                      psautco, prautco, evpco
81     !
82     !
83           real (kind=kind_phys) ERR(IM),      ERS(IM),     PRECRL(IM)
84          &,                     PRECSL(IM),   PRECRL1(IM), PRECSL1(IM)
85          &,                     RQ(IM),       CONDT(IM)
86          &,                     CONDE(IM),    RCONDE(IM),  TMT0(IM)
87          &,                     WMIN(IM,KM),  WMINK(IM),   PRES(IM)
88          &,                     WMINI(IM,KM), CCR(IM),     CCLIM(KM)
89          &,                     TT(IM),       QQ(IM),      WW(IM)
90          &,                     WFIX(KM),     U00K(IM,KM), ES(IM)
91          &,                     Zaodt
92     !
93           integer IW(IM,KM), IPR(IM), IWL(IM),     IWL1(IM)
94     !
95            LOGICAL COMPUT(IM)
96            logical lprnt
97     !
98           real (kind=kind_phys) ke,   rdt,  us, cclimit, climit, cws, csm1
99          &,                     crs1, crs2, cr, aa2,     dtcp,   c00, cmr
100          &,                     tem,  c1,   c2, wwn
101     !    &,                     tem,  c1,   c2, u00b,    u00t,   wwn
102          &,                     precrk, precsk, pres1,   qk,     qw,  qi
103          &,                     ai,     bi, qint, fiw, wws, cwmk, expf
104          &,                     psaut, psaci, amaxcm, tem1, tem2
105          &,                     tmt0k, tmt15, psm1, psm2, ppr
106          &,                     rprs,  erk,   pps, sid, rid, amaxps
107          &,                     praut, pracw, fi, qc, amaxrq, rqkll
108           integer i, k, ihpr, n
109     !
110     !-----------------------Preliminaries ---------------------------------
111     !
112     !     DO K=1,KM
113     !       DO I=1,IM
114     !         CLL(I,K) = 0.0
115     !       ENDDO
116     !     ENDDO
117     !
118           RDT     = H1 / DT
119     !     KE      = 2.0E-5  ! commented on 09/10/99  -- OPR value
120     !     KE      = 2.0E-6
121     !     KE      = 1.0E-5
122     !!!   KE      = 5.0E-5
123     !!    KE      = 7.0E-5
124           KE      = evpco
125     !     KE      = 7.0E-5
126           US      = H1
127           CCLIMIT = 1.0E-3
128           CLIMIT  = 1.0E-20
129           CWS     = 0.025
130     !
131           zaodt   = 800.0 * RDT
132     !
133           CSM1    = 5.0000E-8   * zaodt
134           CRS1    = 5.00000E-6  * zaodt
135           CRS2    = 6.66600E-10 * zaodt
136           CR      = 5.0E-4      * zaodt
137           AA2     = 1.25E-3     * zaodt
138     !
139           ke      = ke * sqrt(rdt)
140     !     ke      = ke * sqrt(zaodt)
141     !
142           DTCP    = DT * RCP
143     !
144     !     C00 = 1.5E-1 * DT
145     !     C00 = 10.0E-1 * DT
146     !     C00 = 3.0E-1 * DT          !05/09/2000
147     !     C00 = 1.0E-4 * DT          !05/09/2000
148           C00 = prautco * DT         !05/09/2000
149           CMR = 1.0 / 3.0E-4
150     !     CMR = 1.0 / 5.0E-4
151     !     C1  = 100.0
152           C1  = 300.0
153           C2  = 0.5
154     !
155     !
156     !--------CALCULATE C0 AND CMR USING LC AT PREVIOUS STEP-----------------
157     !
158           DO K=1,KM
159             DO I=1,IM
160               tem   = (prsl(i,k)*0.00001)
161     !         tem   = sqrt(tem)
162               IW(I,K)    = 0.0
163               wmin(i,k)  = 1.0e-5 * tem
164               wmini(i,k) = 1.0e-5 * tem       ! Testing for RAS
165     !         wmin(i,k)  = 5.0e-6 * tem       ! Testing 
166     !         wmini(i,k) = 5.0e-6 * tem       ! Testing
167     !         wmin(i,k)  = 3.0e-6 * tem       ! Testing 
168     !         wmini(i,k) = 3.0e-6 * tem       ! Testing
169     !         wmini(i,k) = 1.0e-6 * tem       ! for SAS
170     
171               rainp(i,k) = 0.0
172     
173             ENDDO
174           ENDDO
175           DO I=1,IM
176     !       C0(I)  = 1.5E-1
177     !       CMR(I) = 3.0E-4
178     !
179             IWL1(I)    = 0
180             PRECRL1(I) = D00
181             PRECSL1(I) = D00
182             COMPUT(I)  = .FALSE.
183             RN(I)      = D00
184             SR(I)      = D00
185             ccr(i)     = D00
186     !
187             RNP(I)     = D00
188           ENDDO
189     !------------SELECT COLUMNS WHERE RAIN CAN BE PRODUCED--------------
190           DO K=1, KM-1
191             DO I=1,IM
192               tem = min(wmin(i,k), wmini(i,k))
193               IF (CWM(I,K) .GT. tem) COMPUT(I) = .TRUE.
194             ENDDO
195           ENDDO
196           IHPR = 0
197           DO I=1,IM
198             IF (COMPUT(I)) THEN
199                IHPR      = IHPR + 1
200                IPR(IHPR) = I
201             ENDIF
202           ENDDO
203     !***********************************************************************
204     !-----------------BEGINING OF PRECIPITATION CALCULATION-----------------
205     !***********************************************************************
206     !     DO K=KM-1,2,-1
207           DO K=KM,1,-1
208             DO N=1,IHPR
209               PRECRL(N) = PRECRL1(N)
210               PRECSL(N) = PRECSL1(N)
211               ERR  (N)  = D00
212               ERS  (N)  = D00
213               IWL  (N)  = 0
214     !
215               I         = IPR(N)
216               TT(N)     = T(I,K)
217               QQ(N)     = Q(I,K)
218               WW(N)     = CWM(I,K)
219               WMINK(N)  = WMIN(I,K)
220               PRES(N)   = prSL(I,K)
221     !
222               PRECRK = MAX(cons_0,    PRECRL1(N))
223               PRECSK = MAX(cons_0,    PRECSL1(N))
224               WWN    = MAX(WW(N), CLIMIT)
225     !         IF (WWN .GT. WMINK(N) .OR. (PRECRK+PRECSK) .GT. D00) THEN
226               IF (WWN .GT. CLIMIT .OR. (PRECRK+PRECSK) .GT. D00) THEN
227                 COMPUT(N) = .TRUE.
228               ELSE
229                 COMPUT(N) = .FALSE.
230               ENDIF
231             ENDDO
232     !
233     !       es(1:IHPR) = fpvs(TT(1:IHPR))
234             DO N=1,IHPR
235               IF (COMPUT(N)) THEN
236                 I = IPR(N)
237                 CONDE(N)  = (DT/G) * DEL(I,K)
238                 CONDT(N)  = CONDE(N) * RDT
239                 RCONDE(N) = H1 / CONDE(N)
240                 QK        = MAX(EPSQ,  QQ(N))
241                 TMT0(N)   = TT(N) - 273.16
242                 WWN       = MAX(WW(N), CLIMIT)
243     !
244     !           PL = PRES(N) * 0.01
245     !           CALL QSATD(TT(N), PL, QC)
246     !           RQ(N) = MAX(QQ(N), EPSQ) / MAX(QC, 1.0E-10)
247     !           RQ(N) = MAX(1.0E-10, RQ(N))           ! -- RELATIVE HUMIDITY---
248     !
249     !  the global qsat computation is done in Pa
250                 pres1   = pres(n) 
251     !           QW      = es(N)
252                 QW      = min(pres1, fpvs(TT(N)))
253                 QW      = EPS * QW / (PRES1 + EPSM1 * QW)
254                 QW      = MAX(QW,EPSQ)
255     !
256     !           TMT15 = MIN(TMT0(N), cons_m15)
257     !           AI    = 0.008855
258     !           BI    = 1.0
259     !           IF (TMT0(N) .LT. -20.0) THEN
260     !             AI = 0.007225
261     !             BI = 0.9674
262     !           ENDIF
263     !           QI   = QW * (BI + AI*MIN(TMT0(N),cons_0))
264     !           QINT = QW * (1.-0.00032*TMT15*(TMT15+15.))
265     !
266                 qi   = qw
267                 qint = qw
268     !           IF (TMT0(N).LE.-40.) QINT = QI
269     !
270     !-------------------ICE-WATER ID NUMBER IW------------------------------
271                 IF(TMT0(N).LT.-15.) THEN
272                    FI = QK - U00K(I,K)*QI
273                    IF(FI.GT.D00.OR.WWN.GT.CLIMIT) THEN
274                       IWL(N) = 1
275                    ELSE
276                       IWL(N) = 0
277                    ENDIF
278     !           ENDIF
279                 ELSEIF (TMT0(N).GE.0.) THEN
280                    IWL(N) = 0
281     !
282     !           IF(TMT0(N).LT.0.0.AND.TMT0(N).GE.-15.0) THEN
283                 ELSE
284                   IWL(N) = 0
285                   IF(IWL1(N).EQ.1.AND.WWN.GT.CLIMIT) IWL(N)=1
286                 ENDIF
287     !
288     !           IF(TMT0(N).GE.0.) THEN
289     !              IWL(N) = 0
290     !           ENDIF
291     !----------------THE SATUATION SPECIFIC HUMIDITY------------------------
292                 FIW   = FLOAT(IWL(N))
293                 QC    = (H1-FIW)*QINT + FIW*QI
294     !----------------THE RELATIVE HUMIDITY----------------------------------
295                 IF(QC .LE. 1.0E-10) THEN
296                    RQ(N) = D00
297                 ELSE
298                    RQ(N) = QK / QC
299                 ENDIF
300     !----------------CLOUD COVER RATIO CCR----------------------------------
301                 IF(RQ(N).LT.U00K(I,K)) THEN
302                        CCR(N)=D00
303                 ELSEIF(RQ(N).GE.US) THEN
304                        CCR(N)=US
305                 ELSE
306                      RQKLL=MIN(US,RQ(N))
307                      CCR(N)= H1-SQRT((US-RQKLL)/(US-U00K(I,K)))
308                 ENDIF
309     !
310               ENDIF
311             ENDDO
312     !-------------------ICE-WATER ID NUMBER IWL------------------------------
313     !       DO N=1,IHPR
314     !         IF (COMPUT(N) .AND.  (WW(N) .GT. CLIMIT)) THEN
315     !           IF (TMT0(N) .LT. -15.0
316     !    *         .OR. (TMT0(N) .LT. 0.0 .AND. IWL1(N) .EQ. 1))
317     !    *                                      IWL(N) = 1
318     !             CLL(IPR(N),K) = 1.0                           ! Cloud Cover!
319     !             CLL(IPR(N),K) = MIN(1.0, WW(N)*CCLIM(K))      ! Cloud Cover!
320     !         ENDIF
321     !       ENDDO
322     !
323     !---   PRECIPITATION PRODUCTION --  Auto Conversion and Accretion
324     !
325             DO N=1,IHPR
326               IF (COMPUT(N) .AND. CCR(N) .GT. 0.0) THEN
327                 WWS    = WW(N)
328                 CWMK   = MAX(cons_0, WWS)
329     !           AMAXCM = MAX(cons_0, CWMK - WMINK(N))
330                 IF (IWL(N) .EQ. 1) THEN                 !  Ice Phase
331                    AMAXCM = MAX(cons_0, CWMK - WMINI(IPR(N),K))
332                    EXPF      = DT * EXP(0.025*TMT0(N))
333                    PSAUT     = MIN(CWMK, psautco*EXPF*AMAXCM)
334     
335     !              PSAUT     = MIN(CWMK, 2.0E-3*EXPF*AMAXCM)
336     !              PSAUT     = MIN(CWMK, 1.0E-3*EXPF*AMAXCM)
337     !              PSAUT     = MIN(CWMK, 7.5E-4*EXPF*AMAXCM)
338     !!!!!!!        PSAUT     = MIN(CWMK, 7.0E-4*EXPF*AMAXCM)
339     !b             PSAUT     = MIN(CWMK, 6.5E-4*EXPF*AMAXCM)
340     !!!!           PSAUT     = MIN(CWMK, 6.0E-4*EXPF*AMAXCM)
341     !              PSAUT     = MIN(CWMK, 5.0E-4*EXPF*AMAXCM)
342     !              PSAUT     = MIN(CWMK, 4.0E-4*EXPF*AMAXCM)
343     
344                    WW(N)     = WW(N) - PSAUT
345                    CWMK      = MAX(cons_0, WW(N))
346     !              CWMK      = MAX(cons_0, WW(N)-wmini(ipr(n),k))
347                    PSACI     = MIN(CWMK, AA2*EXPF*PRECSL1(N)*CWMK)
348     
349                    WW(N)     = WW(N) - PSACI
350      
351                    PRECSL(N) = PRECSL(N) + (WWS - WW(N)) * CONDT(N)
352                 ELSE                                    !  Liquid Water
353     !
354     !          For using Sundqvist precip formulation of rain
355     !
356     !              AMAXCM    = MAX(cons_0, CWMK - WMINK(N))
357                    AMAXCM    = CWMK
358                    TEM1      = PRECSL1(N) + PRECRL1(N)
359                    TEM2      = MIN(MAX(cons_0, 268.0-TT(N)), cons_20)
360                    TEM       = (1.0+C1*SQRT(TEM1*RDT)) * (1+C2*SQRT(TEM2))
361     !
362                    TEM2      = AMAXCM * CMR * TEM / max(CCR(N),cons_p01)
363                    TEM2      = MIN(cons_50, TEM2*TEM2)
364                    PRAUT     = C00  * TEM * AMAXCM * (1.0-EXP(-TEM2))
365                    PRAUT     = MIN(PRAUT, CWMK)
366                    WW(N)     = WW(N) - PRAUT
367     !
368     !          Below is for Zhao's precip formulation (water)
369     !
370     !              AMAXCM    = MAX(cons_0, CWMK - WMINK(N))
371     !              PRAUT     = MIN(CWMK, C00*AMAXCM*AMAXCM)
372     !              WW(N)     = WW(N) - PRAUT
373     !
374     !              CWMK      = MAX(cons_0, WW(N))
375     !              TEM1      = PRECSL1(N) + PRECRL1(N)
376     !              PRACW     = MIN(CWMK, CR*DT*TEM1*CWMK)
377     !              WW(N)     = WW(N) - PRACW
378     !
379                    PRECRL(N) = PRECRL(N) + (WWS - WW(N)) * CONDT(N)
380     !
381     !hchuang code change [+1L] : add record to record information in vertical
382     ! TURN RNP in unit of WW (CWM and Q, kg/kg ???)
383                    RNP(N) = RNP(N) + (WWS - WW(N))
384                 ENDIF
385               ENDIF
386             ENDDO
387     !
388     !-----EVAPORATION OF PRECIPITATION-------------------------
389     !**** ERR & ERS POSITIVE--->EVAPORATION-- NEGTIVE--->CONDENSATION
390     !
391             DO N=1,IHPR
392               IF (COMPUT(N)) THEN
393                 I      = IPR(N)
394                 QK     = MAX(EPSQ,  QQ(N))
395                 TMT0K  = MAX(cons_m30, TMT0(N))
396                 PRECRK = MAX(cons_0,    PRECRL(N))
397                 PRECSK = MAX(cons_0,    PRECSL(N))
398                 AMAXRQ = MAX(cons_0,    U00K(I,K)-RQ(N)) * CONDE(N)
399     !----------------------------------------------------------------------
400     ! INCREASE THE EVAPORATION FOR STRONG/LIGHT PREC
401     !----------------------------------------------------------------------
402                 PPR    = KE * AMAXRQ * SQRT(PRECRK)
403     !           PPR    = KE * AMAXRQ * SQRT(PRECRK*RDT)
404                 IF (TMT0(N) .GE. 0.) THEN
405                   PPS = 0.
406                 ELSE
407                   PPS = (CRS1+CRS2*TMT0K) * AMAXRQ * PRECSK / U00K(I,K)
408                 END IF
409     !---------------CORRECT IF OVER-EVAPO./COND. OCCURS--------------------
410                 ERK=PRECRK+PRECSK
411                 IF(RQ(N).GE.1.0E-10)  ERK = AMAXRQ * QK * RDT / RQ(N)
412                 IF (PPR+PPS .GT. ABS(ERK)) THEN
413                    RPRS   = ERK / (PRECRK+PRECSK)
414                    PPR    = PRECRK * RPRS
415                    PPS    = PRECSK * RPRS
416                 ENDIF
417                 PPR       = MIN(PPR, PRECRK)
418                 PPS       = MIN(PPS, PRECSK)
419                 ERR(N)    = PPR * RCONDE(N)
420                 ERS(N)    = PPS * RCONDE(N)
421                 PRECRL(N) = PRECRL(N) - PPR
422     !hchuang code change [+1L] : add record to record information in vertical
423     ! Use ERR for kg/kg/DT not the PPR (mm/DT=kg/m2/DT)
424     !
425                 RNP(N) = RNP(N) - ERR(N)
426     !
427                 PRECSL(N) = PRECSL(N) - PPS
428               ENDIF
429             ENDDO
430     !--------------------MELTING OF THE SNOW--------------------------------
431             DO N=1,IHPR
432               IF (COMPUT(N)) THEN
433                 IF (TMT0(N) .GT. 0.) THEN
434                    AMAXPS = MAX(cons_0,    PRECSL(N))
435                    PSM1   = CSM1 * TMT0(N) * TMT0(N) * AMAXPS
436                    PSM2   = CWS * CR * MAX(cons_0, WW(N)) * AMAXPS
437                    PPR    = (PSM1 + PSM2) * CONDE(N)
438                    IF (PPR .GT. AMAXPS) THEN
439                      PPR  = AMAXPS
440                      PSM1 = AMAXPS * RCONDE(N)
441                    ENDIF
442                    PRECRL(N) = PRECRL(N) + PPR
443     !
444     !hchuang code change [+1L] : add record to record information in vertical
445     ! TURN PPR (mm/DT=kg/m2/DT) to kg/kg/DT -> PPR/air density (kg/m3)
446                    RNP(N) = RNP(N) + PPR * RCONDE(N)
447     !
448                    PRECSL(N) = PRECSL(N) - PPR
449                 ELSE
450                    PSM1 = D00
451                 ENDIF
452     !
453     !---------------UPDATE T AND Q------------------------------------------
454                 TT(N) = TT(N) - DTCP * (ELWV*ERR(N)+ELIV*ERS(N)+ELIW*PSM1)
455                 QQ(N) = QQ(N) + DT * (ERR(N)+ERS(N))
456               ENDIF
457             ENDDO
458     !
459             DO N=1,IHPR
460               IWL1(N)    = IWL(N)
461               PRECRL1(N) = MAX(cons_0, PRECRL(N))
462               PRECSL1(N) = MAX(cons_0, PRECSL(N))
463               I          = IPR(N)
464               T(I,K)     = TT(N)
465               Q(I,K)     = QQ(N)
466               CWM(I,K)   = WW(N)
467               IW(I,K)    = IWL(N)
468     !hchuang code change [+1L] : add record to record information in vertical
469     ! RNP = PRECRL1*RCONDE(N) unit in kg/kg/DT
470     !
471               RAINP(I,K) = RNP(N)
472             ENDDO
473     !
474     !  move water from vapor to liquid should the liquid amount be negative
475     !
476             do i = 1, im
477               if (cwm(i,k) < 0.) then
478                 tem      = q(i,k) + cwm(i,k)
479                 if (tem >= 0.0) then
480                   q(i,k)   = tem
481                   t(i,k)   = t(i,k) - elwv * rcp * cwm(i,k)
482                   cwm(i,k) = 0.
483                 elseif (q(i,k) > 0.0) then
484                   cwm(i,k) = tem
485                   t(i,k)   = t(i,k) + elwv * rcp * q(i,k)
486                   q(i,k)   = 0.0
487                 endif
488               endif
489             enddo
490     !
491           ENDDO                               ! K loop ends here!
492     !**********************************************************************
493     !-----------------------END OF PRECIPITATION PROCESSES-----------------
494     !**********************************************************************
495     !
496           DO N=1,IHPR
497             I = IPR(N)
498             RN(I) = (PRECRL1(N)  + PRECSL1(N)) * RROW  ! Precip at surface
499     !
500     !----SR=1 IF SFC PREC IS RAIN ; ----SR=-1 IF SFC PREC IS SNOW
501     !----SR=0 FOR BOTH OF THEM OR NO SFC PREC
502     !
503             RID = 0.
504             SID = 0.
505             IF (PRECRL1(N) .GE. 1.E-13) RID = 1.
506             IF (PRECSL1(N) .GE. 1.E-13) SID = -1.
507             SR(I) = RID + SID  ! SR=1 --> Rain, SR=-1 -->Snow, SR=0 -->Both
508           ENDDO
509     !
510           RETURN
511           END
512