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

1           SUBROUTINE LRGSCL(IX,IM,KM,DT,T1,Q1,PRSL,DEL,PRSLK,RAIN,CLW)
2     !
3           USE MACHINE     , ONLY : kind_phys
4           USE FUNCPHYS , ONLY : fpvs, ftdp, fthe, stma, ftlcl
5           USE PHYSCONS, HVAP => con_HVAP, CP => con_CP, RV => con_RV
6          &,             EPS => con_eps, EPSM1 => con_epsm1, ROCP => con_ROCP
7          &,             grav => con_g
8           implicit none
9     !
10     !     include 'constant.h'
11     !
12           integer IX , IM, KM
13           real(kind=kind_phys) T1(IX,KM),  Q1(IX,KM),    PRSL(IX,KM),
14          &                     DEL(IX,KM), PRSLK(IX,KM), RAIN(IM),
15          &                     CLW(IM,KM), DT
16     !
17     !  LOCAL VARIABLES
18     !
19           integer              k, kmax, I
20           real(kind=kind_phys) dpovg,  EI,    el2orc,
21          &                     elocp,
22          &                     pk,     qcond, qevap,
23          &                     rnevap, SLKLCL,TDPD,
24          &                     THELCL, TLCL,  val0,  val1
25     !
26     !
27     !  PHYSICAL PARAMETERS
28           PARAMETER(ELOCP=HVAP/CP,   EL2ORC=HVAP*HVAP/(RV*CP))
29     !
30     !
31           real(kind=kind_phys) TO(IM,KM),   QO(IM,KM),   QS(IM,KM),
32          &                     THE(IM,KM),  DQ(IM,KM),   RAINLVL(IM,KM),
33          &                     ES(IM,KM),   DQINT(IM),   PINT(IM),
34          &                     DELQBAR(IM), DELTBAR(IM), THEBAR(IM),
35          &                     THEINT(IM)
36           integer              KMLEV(IM,KM), KE(IM), KK(IM), KS(IM)
37           LOGICAL FLG(IM), TOPFLG(IM), TOTFLG
38     cc
39     cc--------------------------------------------------------------------
40     cc
41           real(kind=kind_phys) cons_0               !constant
42           real(kind=kind_phys) cons_1pdm8           !constant
43     cc
44           cons_0          =         0.d0            !constant
45           cons_1pdm8      =         1.d-8           !constant
46     cc
47     cc--------------------------------------------------------------------
48     cc
49           KMAX = KM
50           DO K = 1, KM
51             do i=1,im
52               IF (PRSL(I,K) .GT. 6000.0) KMAX = K + 1
53             enddo
54           ENDDO
55     C
56     C   SURFACE PRESSURE UNIT IS CB
57     C
58           DO I=1,IM
59     !       PSK(I)     = FPKAP(PS(I))
60             RAIN(I)    = 0.
61             DELTBAR(I) = 0.
62             DELQBAR(I) = 0.
63             FLG(I)     = .FALSE.
64             TOPFLG(I)  = .FALSE.
65             KE(I)      = kmax + 1
66             KS(I)      = 0
67           ENDDO
68           TOTFLG = .FALSE.
69     C
70     C  COLUMN VARIABLES
71     C  PRSL IS PRESSURE OF THE LAYER (Pa)
72     C  TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
73     C  QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
74     C
75           DO K = 1, KMAX
76            DO I=1,IM
77     !       PFLD(I,k) = PS(I) * SL(K)
78             TO(I,k) = T1(I,k)
79             QO(I,k) = Q1(I,k)
80             ENDDO
81           ENDDO
82     C
83     C  MODEL CONSISTENT SATURATION MIXING RATIO
84     C
85     !     es(:,:) = 0.001 * fpvs(t1(1:IM,:))        ! fpvs in Pa
86           DO K = 1, KMAX
87             DO I=1,IM
88               es(I,k) = min(PRSL(I,k), fpvs(t1(I,k)))   ! fpvs in Pa
89               QS(I,k) = EPS * ES(I,k) / (PRSL(I,k) + EPSM1*ES(I,k))
90               QS(I,k) = MAX(QS(I,k),cons_1pdm8)     !constant
91             ENDDO
92           ENDDO
93           DO K = 1, KMAX
94             DO I=1,IM
95               IF(QO(I,k).GT.QS(I,k)) FLG(I) = .TRUE.
96             ENDDO
97           ENDDO
98     !!
99           DO I=1,IM
100             IF(FLG(I)) TOTFLG = .TRUE.
101           ENDDO
102           IF(.NOT.TOTFLG) RETURN
103     !!
104           DO K = 1, KMAX
105             DO I = 1, IM
106               DQ(I,k) = 0.
107               THE(I,k) = TO(I,k)
108             ENDDO
109           ENDDO
110     C
111     C  COMPUTE THETA-E
112     C
113           DO K = 1, KMAX
114             DO I = 1, IM
115               IF(FLG(I)) THEN
116     !           PK = PSK(I) * SLK(K)
117                 PK = PRSLK(I,K)
118                 THE(I,k) = FTHE(TO(I,k),PK)
119                 IF(THE(I,k).EQ.0.) THEN
120                   THE(I,k) = TO(I,k) / PK
121                 ENDIF
122     C         THE(I,k) = TO(I,k) * ((PRSL(I,k)-ES(I,k))*.01) ** (-ROCP)
123     C    &             * EXP(ELOCP * QS(I,k) / TO(I,k))
124                 DQ(I,k) = QO(I,k)- QS(I,k)
125     C
126     C  MODIFICATION OF THETA-E FOR SUPER-SATURATION
127     C
128                 THE(I,k)= THE(I,k) * (1. + HVAP*MAX(DQ(I,k),cons_0)     !constant
129          &              /(CP*TO(I,k)))
130               ENDIF
131             ENDDO
132           ENDDO
133           DO K = 1, KMAX
134             DO I = 1, IM
135               KMLEV(I,k) = 0
136               RAINLVL(I,k) = 0.
137             ENDDO
138           ENDDO
139     C
140     C  STARTING POINT OF ADJUSTMENT
141     C
142           K = 1
143           DO I = 1, IM
144             KK(I)     = 0
145             DQINT(I)  = 0.
146             THEINT(I) = 0.
147             THEBAR(I) = 0.
148             PINT(I)   = 0.
149     C
150     C  FOR CONDITIONALLY UNSTABLE AND SUPERSATURATED LAYERS,
151     C    OBTAIN INTEGRATED THETA AND Q-QS
152     C
153     C  KMLEV KEEPS TRACK OF THE NUMBER OF LAYERS THAT SATISFIES
154     C    THE CONDITION FOR ADJUSTMENT
155     C
156             IF(DQ(I,k).GT.0..AND.THE(I,k).GE.THE(I,K+1).AND.FLG(I)) THEN
157               DQINT(I)   = DQINT(I)  + DQ(I,k) * DEL(I,K)
158               THEINT(I)  = THEINT(I) + THE(I,k) * DEL(I,K)
159               PINT(I)    = PINT(I)   + DEL(I,K)
160               KK(I)      = KK(I)     + 1
161               KMLEV(I,k) = KK(I)
162             ENDIF
163           ENDDO
164           DO K = 2, KMAX - 1
165             DO I = 1, IM
166               IF(DQ(I,k).GT.0..AND.THE(I,k).GE.THE(I,K+1).AND.FLG(I)) THEN
167                 DQINT(I)   = DQINT(I)  + DQ(I,k) * DEL(I,K)
168                 THEINT(I)  = THEINT(I) + THE(I,k) * DEL(I,K)
169                 PINT(I)    = PINT(I)   + DEL(I,K)
170                 KK(I)      = KK(I)     + 1
171                 KMLEV(I,k) = KK(I)
172               ENDIF
173             ENDDO
174             DO I = 1, IM
175               IF(PINT(I).GT.0.)THEBAR(I) = THEINT(I) / PINT(I)
176     C
177     C  IF THE LAYER BELOW SATISFIES THE CONDITION AND THE PRESENT
178     C    LAYER IS COLDER THAN THE ADJSUTED THETA-E,
179     C    THE LAYER IS INCLUDED IF THE INTEGRATED MOISTURE EXCESS
180     C    CAN BE MAINTAINED
181     C
182               IF(KMLEV(I,k).EQ.0.AND.KMLEV(I,K-1).GT.0.AND.
183          &       THEBAR(I).GE.THE(I,k).AND..NOT.TOPFLG(I)) THEN
184                    DQINT(I) = DQINT(I) + DQ(I,k) * DEL(I,K)
185               ENDIF
186               IF(KMLEV(I,k).EQ.0.AND.KMLEV(I,K-1).GT.0.AND.
187          &       THEBAR(I).GE.THE(I,k).AND.DQINT(I).GT.0.
188          &       .AND..NOT.TOPFLG(I)) THEN
189                 KK(I)      = KK(I) + 1
190                 KMLEV(I,k) = KK(I)
191                 TOPFLG(I) = .TRUE.
192     !           PK        = PSK(I) * SLK(K)
193                 EI        = PRSL(I,k) * QO(I,k)
194          &                                     / (EPS - EPSM1*QO(I,k))
195                 EI        = MIN(MAX(EI,cons_1pdm8),ES(I,k))    !constant
196                 TDPD      = MAX(TO(I,k)-FTDP(EI),cons_0)       !constant
197                 TLCL      = FTLCL(TO(I,k), TDPD)
198                 SLKLCL    = PRSLK(I,K) * TLCL / TO(I,k)
199                 THELCL    = FTHE(TLCL,SLKLCL)
200                 IF(THELCL.NE.0.) THEN
201                   THE(I,k) = THELCL
202     C             THE(I,k) = TO(I,k) * ((PRSL(I,k) - EI)*.01) ** (-ROCP)
203     C    &             * EXP(ELOCP * MAX(QO(I,k),1.E-8) / TO(I,k))
204                 ENDIF
205                 THEINT(I) = THEINT(I) + THE(I,k) * DEL(I,K)
206                 PINT(I)   = PINT(I) + DEL(I,K)
207               ENDIF
208             ENDDO
209     C
210     C  RESET THE INTEGRAL IF THE LAYER IS NOT IN THE CLOUD
211     C
212             DO I = 1, IM
213               IF(KMLEV(I,k).EQ.0.AND.KMLEV(I,K-1).GT.0) THEN
214                 THEBAR(I) = THEINT(I) / PINT(I)
215                 DQINT(I) = 0.
216                 THEINT(I) = 0.
217                 PINT(I) = 0.
218                 KK(I) = 0
219                 KS(I) = k - 1
220                 KE(I) = KS(I) - KMLEV(I,k-1) + 1
221                 FLG(I) = .false.
222               ENDIF
223             ENDDO
224           enddo
225     C
226     C  When within A CLOUD LAYER, COMPUTE THE MOIST-ADIABATIC
227     C    (TO AND QO) USING THE AVERAGED THETA-E AND THE RESULTANT RAIN
228     C
229           do k = 1, kmax
230             DO I = 1, IM
231               if(k.ge.KE(I).and.k.le.KS(I)) then
232     !           PK = PSK(I) * SLK(K)
233                 PK = PRSLK(I,K)
234     !           TO(I,k)  = FTMA(THEBAR(I),PK,QO(I,k))
235                 CALL STMA(THEBAR(i),PK,TO(I,k),QO(I,k))
236                 THE(I,k) = THEBAR(I)
237                 QS(I,k)  = QO(I,k)
238                 DPOVG    = DEL(I,K) * (1.0/grav)
239                 RAINLVL(I,k) = (Q1(I,k) - QO(I,k)) * dpovg
240                 DELTBAR(I)   = DELTBAR(I) + (TO(I,k) - T1(I,k)) * dpovg / PK
241                 DELQBAR(I)   = DELQBAR(I) + (QO(I,k) - Q1(I,k)) * dpovg
242               ENDIF
243     C
244     C  THIS STEP TAKES CARE OF STABLE HEATING
245     C
246               IF(KMLEV(I,k).EQ.0.AND.DQ(I,k).GT.0.) THEN
247                 QCOND   = (QO(I,k)-QS(I,k)) /
248          &                (1.+EL2ORC*QS(I,k)/(TO(I,K)*TO(I,K)))
249                 QO(I,k) = QO(I,k) - QCOND
250                 TO(I,k) = TO(I,k) + QCOND * ELOCP
251     !           PK = PSK(I) * SLK(K)
252                 PK = PRSLK(I,K)
253     C           TO(I,k) = FTMA(THE(I,k),PK,QO(I,k))
254                 DPOVG    = DEL(I,K) * (1.0/grav)
255                 RAINLVL(I,k) = (Q1(I,k) - QO(I,k)) * dpovg
256                 DELTBAR(I)   = DELTBAR(I) + (TO(I,k) - T1(I,k)) * dpovg / PK
257                 DELQBAR(I)   = DELQBAR(I) + (QO(I,k) - Q1(I,k)) * dpovg
258                 QS(I,k) = QO(I,k)
259               ENDIF
260             ENDDO
261           ENDDO
262     C
263     C  EVAPORATION OF FALLING RAIN
264     C
265           DO K = KMAX, 1, -1
266             DO I = 1, IM
267               T1(I,k) = TO(I,k)
268               Q1(I,k) = QO(I,k)
269               DPOVG   = DEL(I,K) * (1.0/grav)
270               RAIN(I) = RAIN(I) + RAINLVL(I,k) + CLW(I,k) * DPOVG
271               DQ(I,k) = (QO(I,k) - QS(I,k)) /
272          &              (1. + EL2ORC*QS(I,k)/(TO(I,K)*TO(I,K)))
273               IF(RAIN(I).GT.0..AND.RAINLVL(I,k).LE.0.) THEN
274                 QEVAP      = -DQ(I,k)*(1.-EXP(-0.32*SQRT(DT*RAIN(I))))
275                 RNEVAP     = MIN(QEVAP*DPOVG,RAIN(I))
276                 Q1(I,k)    = Q1(I,k)+RNEVAP/DPOVG
277                 T1(I,k)    = T1(I,k)-RNEVAP/DPOVG*ELOCP
278                 RAIN(I)    = RAIN(I)-RNEVAP
279                 DELTBAR(I) = DELTBAR(I) - RNEVAP * ELOCP
280                 DELQBAR(I) = DELQBAR(I) + RNEVAP
281               ENDIF
282             ENDDO
283           ENDDO
284           DO I = 1, IM
285             RAIN(I) = MAX(RAIN(I),cons_0)     !constant
286           ENDDO
287     !!
288           RETURN
289           END
290