source: trunk/WRF.COMMON/WRFV2/phys/module_cu_kfeta.F @ 3567

Last change on this file since 3567 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 101.2 KB
Line 
1MODULE module_cu_kfeta
2
3   USE module_wrf_error
4
5!--------------------------------------------------------------------
6! Lookup table variables:
7      INTEGER, PARAMETER :: KFNT=250,KFNP=220
8      REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
9      REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
10      REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
11      REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
12! Note:  KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
13!        TPMIX2DD, ENVIRTHT
14! End of Lookup table variables:
15
16CONTAINS
17
18   SUBROUTINE KF_eta_CPS(                                    &
19              ids,ide, jds,jde, kds,kde                      &
20             ,ims,ime, jms,jme, kms,kme                      &
21             ,its,ite, jts,jte, kts,kte                      &
22             ,DT,KTAU,DX                                     &
23             ,rho,RAINCV,NCA                                 &
24             ,U,V,TH,T,W,dz8w,Pcps,pi                        &
25             ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1           &
26             ,EP2,SVP1,SVP2,SVP3,SVPT0                       &
27             ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT       &
28             ,QV                                             &
29            ! optionals
30             ,F_QV    ,F_QC    ,F_QR    ,F_QI    ,F_QS       &
31             ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN            &
32             ,RQICUTEN,RQSCUTEN                              &
33                                                             )
34!
35!-------------------------------------------------------------
36   IMPLICIT NONE
37!-------------------------------------------------------------
38   INTEGER,      INTENT(IN   ) ::                            &
39                                  ids,ide, jds,jde, kds,kde, &
40                                  ims,ime, jms,jme, kms,kme, &
41                                  its,ite, jts,jte, kts,kte
42
43   INTEGER,      INTENT(IN   ) :: STEPCU
44   LOGICAL,      INTENT(IN   ) :: warm_rain
45
46   REAL,         INTENT(IN   ) :: XLV0,XLV1,XLS0,XLS1
47   REAL,         INTENT(IN   ) :: CP,R,G,EP1,EP2
48   REAL,         INTENT(IN   ) :: SVP1,SVP2,SVP3,SVPT0
49
50   INTEGER,      INTENT(IN   ) :: KTAU           
51
52   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
53          INTENT(IN   ) ::                                   &
54                                                          U, &
55                                                          V, &
56                                                          W, &
57                                                         TH, &
58                                                          T, &
59                                                         QV, &
60                                                       dz8w, &
61                                                       Pcps, &
62                                                        rho, &
63                                                         pi
64!
65   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
66          INTENT(INOUT) ::                                   &
67                                                      W0AVG
68
69   REAL,  INTENT(IN   ) :: DT, DX
70!
71   REAL, DIMENSION( ims:ime , jms:jme ),                     &
72          INTENT(INOUT) ::                           RAINCV
73
74   REAL,    DIMENSION( ims:ime , jms:jme ),                  &
75            INTENT(INOUT) ::                            NCA
76
77   REAL, DIMENSION( ims:ime , jms:jme ),                     &
78          INTENT(OUT) ::                              CUBOT, &
79                                                      CUTOP   
80
81   LOGICAL, DIMENSION( ims:ime , jms:jme ),                  &
82          INTENT(INOUT) :: CU_ACT_FLAG
83
84!
85! Optional arguments
86!
87
88   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),           &
89         OPTIONAL,                                           &
90         INTENT(INOUT) ::                                    &
91                                                   RTHCUTEN, &
92                                                   RQVCUTEN, &
93                                                   RQCCUTEN, &
94                                                   RQRCUTEN, &
95                                                   RQICUTEN, &
96                                                   RQSCUTEN
97
98!
99! Flags relating to the optional tendency arrays declared above
100! Models that carry the optional tendencies will provdide the
101! optional arguments at compile time; these flags all the model
102! to determine at run-time whether a particular tracer is in
103! use or not.
104!
105   LOGICAL, OPTIONAL ::                                      &
106                                                   F_QV      &
107                                                  ,F_QC      &
108                                                  ,F_QR      &
109                                                  ,F_QI      &
110                                                  ,F_QS
111
112
113! LOCAL VARS
114
115   LOGICAL :: flag_qr, flag_qi, flag_qs
116
117   REAL, DIMENSION( kts:kte ) ::                             &
118                                                        U1D, &
119                                                        V1D, &
120                                                        T1D, &
121                                                       DZ1D, &
122                                                       QV1D, &
123                                                        P1D, &
124                                                      RHO1D, &
125                                                    W0AVG1D
126
127   REAL, DIMENSION( kts:kte )::                              &
128                                                       DQDT, &
129                                                      DQIDT, &
130                                                      DQCDT, &
131                                                      DQRDT, &
132                                                      DQSDT, &
133                                                       DTDT
134
135   REAL    ::         TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp,RTHCUMAX
136
137   INTEGER :: i,j,k,NTST,ICLDCK
138!
139   DXSQ=DX*DX
140
141!----------------------
142   NTST=STEPCU
143   TST=float(NTST*2)
144   flag_qr = .FALSE.
145   flag_qi = .FALSE.
146   flag_qs = .FALSE.
147   IF ( PRESENT(F_QR) ) flag_qr = F_QR
148   IF ( PRESENT(F_QI) ) flag_qi = F_QI
149   IF ( PRESENT(F_QS) ) flag_qs = F_QS
150!
151  DO J = jts,jte
152      DO K=kts,kte
153         DO I= its,ite
154!            SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
155!            TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
156!            RHOE=Pcps(I,K,J)/(R*TV)
157!            W0=-101.9368*SCR1/RHOE
158            W0=0.5*(w(I,K,J)+w(I,K+1,J))
159            W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
160         ENDDO
161      ENDDO
162   ENDDO
163!
164!...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
165!
166!----------------------
167   ICLDCK=MOD(KTAU,NTST)
168   IF(ICLDCK.EQ.0 .or. KTAU .eq. 1) then
169!
170     DO J = jts,jte
171     DO I= its,ite
172        CU_ACT_FLAG(i,j) = .true.
173     ENDDO
174     ENDDO
175
176     DO J = jts,jte
177       DO I=its,ite
178
179         IF ( NINT(NCA(I,J)) .gt. 0 ) then
180            CU_ACT_FLAG(i,j) = .false.
181         ELSE
182
183            DO k=kts,kte
184               DQDT(k)=0.
185               DQIDT(k)=0.
186               DQCDT(k)=0.
187               DQRDT(k)=0.
188               DQSDT(k)=0.
189               DTDT(k)=0.
190            ENDDO
191            RAINCV(I,J)=0.
192            CUTOP(I,J)=KTS
193            CUBOT(I,J)=KTE+1
194!
195! assign vars from 3D to 1D
196
197            DO K=kts,kte
198               U1D(K) =U(I,K,J)
199               V1D(K) =V(I,K,J)
200               T1D(K) =T(I,K,J)
201               RHO1D(K) =rho(I,K,J)
202               QV1D(K)=QV(I,K,J)
203               P1D(K) =Pcps(I,K,J)
204               W0AVG1D(K) =W0AVG(I,K,J)
205               DZ1D(k)=dz8w(I,K,J)
206            ENDDO
207            CALL KF_eta_PARA(I, J,                  &
208                 U1D,V1D,T1D,QV1D,P1D,DZ1D,         &
209                 W0AVG1D,DT,DX,DXSQ,RHO1D,          &
210                 XLV0,XLV1,XLS0,XLS1,CP,R,G,        &
211                 EP2,SVP1,SVP2,SVP3,SVPT0,          &
212                 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
213                 RAINCV,NCA,NTST,                   &
214                 flag_QI,flag_QS,warm_rain,         &
215                 CUTOP,CUBOT,                       &
216                 ids,ide, jds,jde, kds,kde,         &
217                 ims,ime, jms,jme, kms,kme,         &
218                 its,ite, jts,jte, kts,kte)
219            IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
220              DO K=kts,kte
221                 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
222!              RTHCUMAX=max(abs(RTHCUTEN(I,K,J)),RTHCUMAX)
223                 RQVCUTEN(I,K,J)=DQDT(K)
224              ENDDO
225            ENDIF
226
227            IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
228              IF( F_QR )THEN
229                DO K=kts,kte
230                   RQRCUTEN(I,K,J)=DQRDT(K)
231                   RQCCUTEN(I,K,J)=DQCDT(K)
232                ENDDO
233              ELSE
234! This is the case for Eta microphysics without 3d rain field
235                DO K=kts,kte
236                   RQRCUTEN(I,K,J)=0.
237                   RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
238                ENDDO
239              ENDIF
240            ENDIF
241
242!......     QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
243
244
245            IF(PRESENT( rqicuten )) THEN
246              IF ( F_QI ) THEN
247                DO K=kts,kte
248                   RQICUTEN(I,K,J)=DQIDT(K)
249                ENDDO
250              ENDIF
251            ENDIF
252
253            IF(PRESENT( rqscuten )) THEN
254              IF ( F_QS ) THEN
255                DO K=kts,kte
256                   RQSCUTEN(I,K,J)=DQSDT(K)
257                ENDDO
258              ENDIF
259            ENDIF
260!
261         ENDIF
262       ENDDO
263     ENDDO
264   ENDIF
265!
266   END SUBROUTINE KF_eta_CPS
267! ****************************************************************************
268!-----------------------------------------------------------
269   SUBROUTINE KF_eta_PARA (I, J,                           &
270                      U0,V0,T0,QV0,P0,DZQ,W0AVG1D,         &
271                      DT,DX,DXSQ,rhoe,                     &
272                      XLV0,XLV1,XLS0,XLS1,CP,R,G,          &
273                      EP2,SVP1,SVP2,SVP3,SVPT0,            &
274                      DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT,   &
275                      RAINCV,NCA,NTST,                     &
276                      F_QI,F_QS,warm_rain,                 &
277                      CUTOP,CUBOT,                         &
278                      ids,ide, jds,jde, kds,kde,           &
279                      ims,ime, jms,jme, kms,kme,           &
280                      its,ite, jts,jte, kts,kte)
281!-----------------------------------------------------------
282!***** The KF scheme that is currently used in experimental runs of EMCs
283!***** Eta model....jsk 8/00
284!
285      IMPLICIT NONE
286!-----------------------------------------------------------
287      INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde, &
288                                ims,ime, jms,jme, kms,kme, &
289                                its,ite, jts,jte, kts,kte, &
290                                I,J,NTST
291          ! ,P_QI,P_QS,P_FIRST_SCALAR
292
293      LOGICAL, INTENT(IN   ) :: F_QI, F_QS
294
295      LOGICAL, INTENT(IN   ) :: warm_rain
296!
297      REAL, DIMENSION( kts:kte ),                          &
298            INTENT(IN   ) ::                           U0, &
299                                                       V0, &
300                                                       T0, &
301                                                      QV0, &
302                                                       P0, &
303                                                     rhoe, &
304                                                      DZQ, &
305                                                  W0AVG1D
306!
307      REAL,  INTENT(IN   ) :: DT,DX,DXSQ
308!
309
310      REAL,  INTENT(IN   ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
311      REAL,  INTENT(IN   ) :: EP2,SVP1,SVP2,SVP3,SVPT0
312
313!
314      REAL, DIMENSION( kts:kte ), INTENT(INOUT) ::         &
315                                                     DQDT, &
316                                                    DQIDT, &
317                                                    DQCDT, &
318                                                    DQRDT, &
319                                                    DQSDT, &
320                                                     DTDT
321
322      REAL,    DIMENSION( ims:ime , jms:jme ),             &
323            INTENT(INOUT) ::                          NCA
324
325      REAL, DIMENSION( ims:ime , jms:jme ),                &
326            INTENT(INOUT) ::                       RAINCV
327
328      REAL, DIMENSION( ims:ime , jms:jme ),                &
329            INTENT(OUT) ::                          CUBOT, &
330                                                    CUTOP
331!
332!...DEFINE LOCAL VARIABLES...
333!
334      REAL, DIMENSION( kts:kte ) ::                        &
335            Q0,Z0,TV0,TU,TVU,QU,TZ,TVD,                    &
336            QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD,      &
337            UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2,             &
338            UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE,          &
339            THTAU,THETEU,THTAD,THETED,QLIQ,QICE,           &
340            QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC,       &
341            DETLQ2,DETIC2,RATIO,RATIO2
342
343
344      REAL, DIMENSION( kts:kte ) ::                        &
345            DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD,              &
346            QDT,FXM,THTAG,THPA,THFXOUT,                    &
347            THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN,           &
348            QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA,              &
349            QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT,            &
350            QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
351
352
353      REAL, DIMENSION( kts:kte+1 ) :: OMG
354      REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
355      REAL, DIMENSION( kts:kte ) ::                        &
356            CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
357
358! LOCAL VARS
359
360      REAL    :: P00,T00,RLF,RHIC,RHBC,PIE,         &
361                 TTFRZ,TBFRZ,C5,RATE
362      REAL    :: GDRY,ROCP,ALIQ,BLIQ,                      &
363                 CLIQ,DLIQ
364      REAL    :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX,   &
365                 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL,     &
366                 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR,   &
367                 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
368                 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD,   &
369                 UPNEW,ABE,WKLCL,TTEMP,FRC1,   &
370                 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
371                 DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2,         &
372                 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1,  &
373                 UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT,           &
374                 THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
375                 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN,   &
376                 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1,     &
377                 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF,   &
378                 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF,     &
379                 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
380                 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1,     &
381                 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF,   &
382                 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR,     &
383                 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2,    &
384                 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
385                 DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
386   REAL    ::    ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
387                 QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
388!
389      INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
390   REAL    :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
391   REAL    :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
392
393      INTEGER :: KX,K,KL
394!
395      INTEGER :: NCHECK
396      INTEGER, DIMENSION (kts:kte) :: KCHECK
397
398      INTEGER :: ISTOP,ML,L5,KMIX,LOW,                     &
399                 LC,MXLAYR,LLFC,NLAYRS,NK,                 &
400                 KPBL,KLCL,LCL,LET,IFLAG,                  &
401                 NK1,LTOP,NJ,LTOP1,                        &
402                 LTOPM1,LVF,KSTART,KMIN,LFS,               &
403                 ND,NIC,LDB,LDT,ND1,NDK,                   &
404                 NM,LMAX,NCOUNT,NOITR,                     &
405                 NSTEP,NTC,NCHM,ISHALL,NSHALL
406      LOGICAL :: IPRNT
407      CHARACTER*1024 message
408!
409      DATA P00,T00/1.E5,273.16/
410      DATA RLF/3.339E5/
411      DATA RHIC,RHBC/1.,0.90/
412      DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
413      DATA RATE/0.03/
414!-----------------------------------------------------------
415      IPRNT=.FALSE.
416      GDRY=-G/CP
417      ROCP=R/CP
418      NSHALL = 0
419      KL=kte
420      KX=kte
421!
422!     ALIQ = 613.3
423!     BLIQ = 17.502
424!     CLIQ = 4780.8
425!     DLIQ = 32.19
426      ALIQ = SVP1*1000.
427      BLIQ = SVP2
428      CLIQ = SVP2*SVPT0
429      DLIQ = SVP3
430!
431!
432!****************************************************************************
433!                                                      ! PPT FB MODS
434!...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER    ! PPT FB MODS
435!...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)     ! PPT FB MODS
436!...FIELD.  "FBFRC" IS THE FRACTION OF AVAILABLE       ! PPT FB MODS
437!...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...        ! PPT FB MODS
438      FBFRC=0.0                                        ! PPT FB MODS
439!...mods to allow shallow convection...
440      NCHM = 0
441      ISHALL = 0
442      DPMIN = 5.E3
443!...
444      P300=P0(1)-30000.
445!
446!...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
447!...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
448!...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
449!
450!...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
451!...FROM BOTTOM-UP IN THE KF SCHEME...
452!
453      ML=0
454!SUE  tmprpsb=1./PSB(I,J)
455!SUE  CELL=PTOP*tmprpsb
456!
457      DO K=1,KX
458!
459!...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
460!
461         ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
462         QES(K)=0.622*ES/(P0(K)-ES)
463         Q0(K)=AMIN1(QES(K),QV0(K))
464         Q0(K)=AMAX1(0.000001,Q0(K))
465         QL0(K)=0.
466         QI0(K)=0.
467         QR0(K)=0.
468         QS0(K)=0.
469         RH(K) = Q0(K)/QES(K)
470         DILFRC(K) = 1.
471         TV0(K)=T0(K)*(1.+0.608*Q0(K))
472!        RHOE(K)=P0(K)/(R*TV0(K))
473!   DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
474         DP(K)=rhoe(k)*g*DZQ(k)
475! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
476! use it for shallow convection...For now, assume it is not available....
477!         TKE(K) = Q2(I,J,NK)
478         TKE(K) = 0.
479         CLDHGT(K) = 0.
480!        IF(P0(K).GE.500E2)L5=K
481         IF(P0(K).GE.0.5*P0(1))L5=K
482         IF(P0(K).GE.P300)LLFC=K
483         IF(T0(K).GT.T00)ML=K
484      ENDDO
485!
486!...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
487        Z0(1)=.5*DZQ(1)
488!cdir novector
489        DO K=2,KL
490          Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
491          DZA(K-1)=Z0(K)-Z0(K-1)
492        ENDDO   
493        DZA(KL)=0.
494!
495!
496!  To save time, specify a pressure interval to move up in sequential
497!  check of different ~50 mb deep groups of adjacent model layers in
498!  the process of identifying updraft source layer (USL).  Note that
499!  this search is terminated as soon as a buoyant parcel is found and
500!  this parcel can produce a cloud greater than specifed minimum depth
501!  (CHMIN)...For now, set interval at 15 mb...
502!
503       NCHECK = 1
504       KCHECK(NCHECK)=1
505       PM15 = P0(1)-15.E2
506       DO K=2,LLFC
507         IF(P0(K).LT.PM15)THEN
508           NCHECK = NCHECK+1
509           KCHECK(NCHECK) = K
510           PM15 = PM15-15.E2
511         ENDIF
512       ENDDO
513!
514       NU=0
515       NUCHM=0
516usl:   DO
517           NU = NU+1
518           IF(NU.GT.NCHECK)THEN
519             IF(ISHALL.EQ.1)THEN
520               CHMAX = 0.
521               NCHM = 0
522               DO NK = 1,NCHECK
523                 NNN=KCHECK(NK)
524                 IF(CLDHGT(NNN).GT.CHMAX)THEN
525                   NCHM = NNN
526                   NUCHM = NK
527                   CHMAX = CLDHGT(NNN)
528                 ENDIF
529               ENDDO
530               NU = NUCHM-1
531               FBFRC=1.
532               CYCLE usl
533             ELSE
534               RETURN
535             ENDIF
536           ENDIF     
537           KMIX = KCHECK(NU)
538           LOW=KMIX
539!...
540           LC = LOW
541!
542!...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
543!...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
544!...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
545!...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
546!   
547           NLAYRS=0
548           DPTHMX=0.
549           NK=LC-1
550           IF ( NK+1 .LT. KTS ) THEN
551             WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
552             CALL wrf_message (TRIM(message))
553           ELSE
554             DO
555               NK=NK+1   
556               IF ( NK .GT. KTE ) THEN
557                 WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
558                 CALL wrf_message (TRIM(message))
559                 EXIT
560               ENDIF
561               DPTHMX=DPTHMX+DP(NK)
562               NLAYRS=NLAYRS+1
563               IF(DPTHMX.GT.DPMIN)THEN
564                 EXIT
565               ENDIF
566             END DO   
567           ENDIF
568           IF(DPTHMX.LT.DPMIN)THEN
569             RETURN
570           ENDIF
571           KPBL=LC+NLAYRS-1   
572!
573!...********************************************************
574!...for computational simplicity without much loss in accuracy,
575!...mix temperature instead of theta for evaluating convective
576!...initiation (triggering) potential...
577!          THMIX=0.
578           TMIX=0.
579           QMIX=0.
580           ZMIX=0.
581           PMIX=0.
582!
583!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
584!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
585!...LAYERS...
586!
587!cdir novector
588           DO NK=LC,KPBL
589             TMIX=TMIX+DP(NK)*T0(NK)
590             QMIX=QMIX+DP(NK)*Q0(NK)
591             ZMIX=ZMIX+DP(NK)*Z0(NK)
592             PMIX=PMIX+DP(NK)*P0(NK)
593           ENDDO   
594!         THMIX=THMIX/DPTHMX
595          TMIX=TMIX/DPTHMX
596          QMIX=QMIX/DPTHMX
597          ZMIX=ZMIX/DPTHMX
598          PMIX=PMIX/DPTHMX
599          EMIX=QMIX*PMIX/(0.622+QMIX)
600!
601!...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
602!
603!        TLOG=ALOG(EMIX/ALIQ)
604! ...calculate dewpoint using lookup table...
605!
606          astrt=1.e-3
607          ainc=0.075
608          a1=emix/aliq
609          tp=(a1-astrt)/ainc
610          indlu=int(tp)+1
611          value=(indlu-1)*ainc+astrt
612          aintrp=(a1-value)/ainc
613          tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
614          TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
615          TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
616          TLCL=AMIN1(TLCL,TMIX)
617          TVLCL=TLCL*(1.+0.608*QMIX)
618          ZLCL = ZMIX+(TLCL-TMIX)/GDRY
619          NK = LC-1
620          DO
621            NK = NK+1
622            KLCL=NK
623            IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
624              EXIT
625            ENDIF
626          ENDDO   
627          IF(NK.GT.KL)THEN
628            RETURN 
629          ENDIF
630          K=KLCL-1
631          DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
632!     
633!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
634!     
635          TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
636          QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
637          TVEN=TENV*(1.+0.608*QENV)
638!     
639!...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
640!...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
641!...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
642!...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
643!...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
644!...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
645!...SUCCESS AT GRID LENGTHS NEAR 25 km.  FOR DIFFERENT GRID-LENGTHS,
646!...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
647!...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
648!     
649          IF(ZLCL.LT.2.E3)THEN
650            WKLCL=0.02*ZLCL/2.E3
651          ELSE
652            WKLCL=0.02
653          ENDIF
654          WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
655          IF(WKL.LT.0.0001)THEN
656            DTLCL=0.
657          ELSE
658            DTLCL=4.64*WKL**0.33
659          ENDIF
660!
661!...for ETA model, give parcel an extra temperature perturbation based
662!...the threshold RH for condensation (U00)...
663!
664!...for now, just assume U00=0.75...
665!...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
666!         U00 = 0.75
667!         IF(U00.lt.1.)THEN
668!           QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
669!           RHLCL = QENV/QSLCL
670!           DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
671!           IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
672!             DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
673!           ELSEIF(RHLCL.GT.0.95)THEN
674!             DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
675!           ELSE
676               DTRH = 0.
677!           ENDIF
678!         ENDIF   
679!         IF(ISHALL.EQ.1)IPRNT=.TRUE.
680!         IPRNT=.TRUE.
681!         IF(TLCL+DTLCL.GT.TENV)GOTO 45
682!
683trigger:  IF(TLCL+DTLCL+DTRH.LT.TENV)THEN   
684!
685! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
686!
687            CYCLE usl
688!
689          ELSE                            ! Parcel is buoyant, determine updraft
690!     
691!...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
692!...EQUIVALENT POTENTIAL TEMPERATURE
693!...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
694!     
695            CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
696!
697!...modify calculation of initial parcel vertical velocity...jsk 11/26/97
698!
699            DTTOT = DTLCL+DTRH
700            IF(DTTOT.GT.1.E-4)THEN
701              GDT=2.*G*DTTOT*500./TVEN
702              WLCL=1.+0.5*SQRT(GDT)
703              WLCL = AMIN1(WLCL,3.)
704            ELSE
705              WLCL=1.
706            ENDIF
707            PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
708            WTW=WLCL*WLCL
709!
710            TVLCL=TLCL*(1.+0.608*QMIX)
711            RHOLCL=PLCL/(R*TVLCL)
712!       
713            LCL=KLCL
714            LET=LCL
715! make RAD a function of background vertical velocity...
716            IF(WKL.LT.0.)THEN
717              RAD = 1000.
718            ELSEIF(WKL.GT.0.1)THEN
719              RAD = 2000.
720            ELSE
721              RAD = 1000.+1000*WKL/0.1
722            ENDIF
723!     
724!*******************************************************************
725!                                                                  *
726!                 COMPUTE UPDRAFT PROPERTIES                       *
727!                                                                  *
728!*******************************************************************
729!     
730!     
731!...
732!...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
733!     
734            WU(K)=WLCL
735            AU0=0.01*DXSQ
736            UMF(K)=RHOLCL*AU0
737            VMFLCL=UMF(K)
738            UPOLD=VMFLCL
739            UPNEW=UPOLD
740!     
741!...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
742!...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
743!...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
744!...PRODUCTION...
745!     
746            RATIO2(K)=0.
747            UER(K)=0.
748            ABE=0.
749            TRPPT=0.
750            TU(K)=TLCL
751            TVU(K)=TVLCL
752            QU(K)=QMIX
753            EQFRC(K)=1.
754            QLIQ(K)=0.
755            QICE(K)=0.
756            QLQOUT(K)=0.
757            QICOUT(K)=0.
758            DETLQ(K)=0.
759            DETIC(K)=0.
760            PPTLIQ(K)=0.
761            PPTICE(K)=0.
762            IFLAG=0
763!     
764!...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
765!...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
766!...FREEZING IS SPECIFIED TO BEGIN.  WITHIN THE GLACIATION
767!...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
768!...PREVIOUS MODEL LEVEL...
769!     
770            TTEMP=TTFRZ
771!     
772!...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
773!...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
774!...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
775!     
776!     
777            EE1=1.
778            UD1=0.
779            REI = 0.
780            DILBE = 0.
781updraft:    DO NK=K,KL-1
782              NK1=NK+1
783              RATIO2(NK1)=RATIO2(NK)
784              FRC1=0.
785              TU(NK1)=T0(NK1)
786              THETEU(NK1)=THETEU(NK)
787              QU(NK1)=QU(NK)
788              QLIQ(NK1)=QLIQ(NK)
789              QICE(NK1)=QICE(NK)
790              call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1),        &
791                     qice(nk1),qnewlq,qnewic,XLV1,XLV0)
792!
793!
794!...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
795!...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
796!...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
797!...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
798!...LIQUID WATER IS FROZEN AT EACH LEVEL...
799!
800              IF(TU(NK1).LE.TTFRZ)THEN
801                IF(TU(NK1).GT.TBFRZ)THEN
802                  IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
803                  FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
804                ELSE
805                  FRC1=1.
806                  IFLAG=1
807                ENDIF
808                TTEMP=TU(NK1)
809!
810!  DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
811!...IS BELOW TTFRZ...
812!
813                QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
814                QNEWIC=QNEWIC+QNEWLQ*FRC1
815                QNEWLQ=QNEWLQ-QNEWLQ*FRC1
816                QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
817                QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
818                CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ,         &
819                          QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
820              ENDIF
821              TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
822!
823!  CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
824!
825              IF(NK.EQ.K)THEN
826                BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
827                BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
828                DZZ=Z0(NK1)-ZLCL
829              ELSE
830                BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
831                BOTERM=2.*DZA(NK)*G*BE/1.5
832                DZZ=DZA(NK)
833              ENDIF
834              ENTERM=2.*REI*WTW/UPOLD
835
836              CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,      &
837                        RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
838!
839!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
840!...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
841!
842              IF(WTW.LT.1.E-3)THEN
843                EXIT
844              ELSE
845                WU(NK1)=SQRT(WTW)
846              ENDIF
847!...Calculate value of THETA-E in environment to entrain into updraft...
848!
849              CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
850!
851!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
852!
853              REI=VMFLCL*DP(NK1)*0.03/RAD
854              TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
855              IF(NK.EQ.K)THEN
856                DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
857              ELSE
858                DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
859              ENDIF
860              IF(DILBE.GT.0.)ABE=ABE+DILBE*G
861!
862!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
863!...ENTRAINMENT (0.5*REI) IS IMPOSED...
864!
865              IF(TVQU(NK1).LE.TV0(NK1))THEN    ! Entrain/Detrain IF BLOCK
866                EE2=0.5
867                UD2=1.
868                EQFRC(NK1)=0.
869              ELSE
870                LET=NK1
871                TTMP=TVQU(NK1)
872!
873!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
874!
875                F1=0.95
876                F2=1.-F1
877                THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
878                QTMP=F1*Q0(NK1)+F2*QU(NK1)
879                TMPLIQ=F2*QLIQ(NK1)
880                TMPICE=F2*QICE(NK1)
881                call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice,        &
882                           qnewlq,qnewic,XLV1,XLV0)
883                TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
884                IF(TU95.GT.TV0(NK1))THEN
885                  EE2=1.
886                  UD2=0.
887                  EQFRC(NK1)=1.0
888                ELSE
889                  F1=0.10
890                  F2=1.-F1
891                  THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
892                  QTMP=F1*Q0(NK1)+F2*QU(NK1)
893                  TMPLIQ=F2*QLIQ(NK1)
894                  TMPICE=F2*QICE(NK1)
895                  call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice,        &
896                               qnewlq,qnewic,XLV1,XLV0)
897                  TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
898                  TVDIFF = ABS(TU10-TVQU(NK1))
899                  IF(TVDIFF.LT.1.e-3)THEN
900                    EE2=1.
901                    UD2=0.
902                    EQFRC(NK1)=1.0
903                  ELSE
904                    EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
905                    EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
906                    EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
907                    IF(EQFRC(NK1).EQ.1)THEN
908                      EE2=1.
909                      UD2=0.
910                    ELSEIF(EQFRC(NK1).EQ.0.)THEN
911                      EE2=0.
912                      UD2=1.
913                    ELSE
914!
915!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
916!   FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
917!
918                      CALL PROF5(EQFRC(NK1),EE2,UD2)
919                    ENDIF
920                  ENDIF
921                ENDIF
922              ENDIF                            ! End of Entrain/Detrain IF BLOCK
923!
924!
925!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
926!   VALUES IN THE LAYER...
927!
928              EE2 = AMAX1(EE2,0.5)
929              UD2 = 1.5*UD2
930              UER(NK1)=0.5*REI*(EE1+EE2)
931              UDR(NK1)=0.5*REI*(UD1+UD2)
932!
933!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
934!   UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
935!
936              IF(UMF(NK)-UDR(NK1).LT.10.)THEN
937!
938!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
939!   FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
940!   First, correct ABE calculation if needed...
941!
942                IF(DILBE.GT.0.)THEN
943                  ABE=ABE-DILBE*G
944                ENDIF
945                LET=NK
946!               WRITE(98,1015)P0(NK1)/100.
947                EXIT
948              ELSE
949                EE1=EE2
950                UD1=UD2
951                UPOLD=UMF(NK)-UDR(NK1)
952                UPNEW=UPOLD+UER(NK1)
953                UMF(NK1)=UPNEW
954                DILFRC(NK1) = UPNEW/UPOLD
955!
956!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
957!...ICE IN THE DETRAINING UPDRAFT MASS...
958!
959                DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
960                DETIC(NK1)=QICE(NK1)*UDR(NK1)
961                QDT(NK1)=QU(NK1)
962                QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
963                THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
964                QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
965                QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
966!
967!...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
968!...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
969!...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
970!...CURRENT MODEL LEVEL...
971!
972                PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
973                PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
974!
975                TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
976                IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
977              ENDIF
978!
979            END DO updraft
980!
981!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
982!   TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
983!   THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BE
984!   THE LET AND CLOUD TOP...
985!     
986!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOC
987!   FIRST BECOMES NEGATIVE...
988!     
989            LTOP=NK
990            CLDHGT(LC)=Z0(LTOP)-ZLCL
991!
992!...Instead of using the same minimum cloud height (for deep convection)
993!...everywhere, try specifying minimum cloud depth as a function of TLCL...
994!
995!
996!
997            IF(TLCL.GT.293.)THEN
998              CHMIN = 4.E3
999            ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
1000              CHMIN = 2.E3 + 100.*(TLCL-273.)
1001            ELSEIF(TLCL.LT.273.)THEN
1002              CHMIN = 2.E3
1003            ENDIF
1004
1005!     
1006!...If cloud top height is less than the specified minimum for deep
1007!...convection, save value to consider this level as source for
1008!...shallow convection, go back up to check next level...
1009!     
1010!...Try specifying minimum cloud depth as a function of TLCL...
1011!
1012!
1013!...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
1014!
1015!...            1.) if there is no CAPE, or
1016!...            2.) cloud top is at model level just above LCL, or
1017!...            3.) cloud top is within updraft source layer, or
1018!...            4.) cloud-top detrainment layer begins within
1019!...                updraft source layer.
1020!
1021            IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN  ! No Convection Allowed
1022              CLDHGT(LC)=0.
1023              DO NK=K,LTOP
1024                UMF(NK)=0.
1025                UDR(NK)=0.
1026                UER(NK)=0.
1027                DETLQ(NK)=0.
1028                DETIC(NK)=0.
1029                PPTLIQ(NK)=0.
1030                PPTICE(NK)=0.
1031              ENDDO
1032!       
1033            ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN      ! Deep Convection allowed
1034              ISHALL=0
1035              EXIT usl
1036            ELSE
1037!
1038!...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
1039              ISHALL = 1
1040              IF(NU.EQ.NUCHM)THEN
1041                EXIT usl               ! Shallow Convection from this layer
1042              ELSE
1043! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
1044                DO NK=K,LTOP
1045                  UMF(NK)=0.
1046                  UDR(NK)=0.
1047                  UER(NK)=0.
1048                  DETLQ(NK)=0.
1049                  DETIC(NK)=0.
1050                  PPTLIQ(NK)=0.
1051                  PPTICE(NK)=0.
1052                ENDDO
1053              ENDIF
1054            ENDIF
1055          ENDIF trigger
1056        END DO usl
1057    IF(ISHALL.EQ.1)THEN
1058      KSTART=MAX0(KPBL,KLCL)
1059      LET=KSTART
1060    endif
1061!     
1062!...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1063!   THIS LEVEL...
1064!     
1065    IF(LET.EQ.LTOP)THEN
1066      UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1067      DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1068      DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1069      UER(LTOP)=0.
1070      UMF(LTOP)=0.
1071    ELSE
1072!     
1073!   BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1074!     
1075      DPTT=0.
1076      DO NJ=LET+1,LTOP
1077        DPTT=DPTT+DP(NJ)
1078      ENDDO
1079      DUMFDP=UMF(LET)/DPTT
1080!     
1081!...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1082!   RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1083!     
1084      DO NK=LET+1,LTOP
1085!
1086!...entrainment is allowed at every level except for LTOP, so disallow
1087!...entrainment at LTOP and adjust entrainment rates between LET and LTOP
1088!...so the the dilution factor due to entyrianment is not changed but
1089!...the actual entrainment rate will change due due forced total
1090!...detrainment in this layer...
1091!
1092        IF(NK.EQ.LTOP)THEN
1093          UDR(NK) = UMF(NK-1)
1094          UER(NK) = 0.
1095          DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
1096          DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
1097        ELSE
1098          UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
1099          UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
1100          UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
1101          DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
1102          DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
1103        ENDIF
1104        IF(NK.GE.LET+2)THEN
1105          TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1106          PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1107          PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1108          TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1109        ENDIF
1110      ENDDO
1111    ENDIF
1112!     
1113! Initialize some arrays below cloud base and above cloud top...
1114!
1115    DO NK=1,K
1116      IF(NK.GE.LC)THEN
1117        IF(NK.EQ.LC)THEN
1118          UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1119          UER(NK)=VMFLCL*DP(NK)/DPTHMX
1120        ELSEIF(NK.LE.KPBL)THEN
1121          UER(NK)=VMFLCL*DP(NK)/DPTHMX
1122          UMF(NK)=UMF(NK-1)+UER(NK)
1123        ELSE
1124          UMF(NK)=VMFLCL
1125          UER(NK)=0.
1126        ENDIF
1127        TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1128        QU(NK)=QMIX
1129        WU(NK)=WLCL
1130      ELSE
1131        TU(NK)=0.
1132        QU(NK)=0.
1133        UMF(NK)=0.
1134        WU(NK)=0.
1135        UER(NK)=0.
1136      ENDIF
1137      UDR(NK)=0.
1138      QDT(NK)=0.
1139      QLIQ(NK)=0.
1140      QICE(NK)=0.
1141      QLQOUT(NK)=0.
1142      QICOUT(NK)=0.
1143      PPTLIQ(NK)=0.
1144      PPTICE(NK)=0.
1145      DETLQ(NK)=0.
1146      DETIC(NK)=0.
1147      RATIO2(NK)=0.
1148      CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
1149      EQFRC(NK)=1.0
1150    ENDDO
1151!     
1152      LTOP1=LTOP+1
1153      LTOPM1=LTOP-1
1154!     
1155!...DEFINE VARIABLES ABOVE CLOUD TOP...
1156!     
1157      DO NK=LTOP1,KX
1158        UMF(NK)=0.
1159        UDR(NK)=0.
1160        UER(NK)=0.
1161        QDT(NK)=0.
1162        QLIQ(NK)=0.
1163        QICE(NK)=0.
1164        QLQOUT(NK)=0.
1165        QICOUT(NK)=0.
1166        DETLQ(NK)=0.
1167        DETIC(NK)=0.
1168        PPTLIQ(NK)=0.
1169        PPTICE(NK)=0.
1170        IF(NK.GT.LTOP1)THEN
1171          TU(NK)=0.
1172          QU(NK)=0.
1173          WU(NK)=0.
1174        ENDIF
1175        THTA0(NK)=0.
1176        THTAU(NK)=0.
1177        EMS(NK)=0.
1178        EMSD(NK)=0.
1179        TG(NK)=T0(NK)
1180        QG(NK)=Q0(NK)
1181        QLG(NK)=0.
1182        QIG(NK)=0.
1183        QRG(NK)=0.
1184        QSG(NK)=0.
1185        OMG(NK)=0.
1186      ENDDO
1187        OMG(KX+1)=0.
1188        DO NK=1,LTOP
1189          EMS(NK)=DP(NK)*DXSQ/G
1190          EMSD(NK)=1./EMS(NK)
1191!     
1192!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCH
1193!     
1194          EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1195          THTAU(NK)=TU(NK)*EXN(NK)
1196          EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1197          THTA0(NK)=T0(NK)*EXN(NK)
1198          DDILFRC(NK) = 1./DILFRC(NK)
1199          OMG(NK)=0.
1200        ENDDO
1201!     IF (XTIME.LT.10.)THEN
1202!      WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1203!    * TMIX-T00,PMIX,QMIX,ABE
1204!      WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1205!    * WLCL,CLDHGT
1206!     ENDIF
1207!     
1208!...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1209!...AND MIDTROPOSPHERE IS USED.
1210!     
1211        WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1212        WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1213        WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1214        VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1215!...for ETA model, DX is a function of location...
1216!       TIMEC=DX(I,J)/VCONV
1217        TIMEC=DX/VCONV
1218        TADVEC=TIMEC
1219        TIMEC=AMAX1(1800.,TIMEC)
1220        TIMEC=AMIN1(3600.,TIMEC)
1221        IF(ISHALL.EQ.1)TIMEC=2400.
1222        NIC=NINT(TIMEC/DT)
1223        TIMEC=FLOAT(NIC)*DT
1224!     
1225!...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1226!     
1227        IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1228          SHSIGN=1.
1229        ELSE
1230          SHSIGN=-1.
1231        ENDIF
1232        VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))*   &
1233            (V0(LTOP)-V0(KLCL))
1234        VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1235        PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1236        PEF=AMAX1(PEF,.2)
1237        PEF=AMIN1(PEF,.9)
1238!     
1239!...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1240!     
1241        CBH=(ZLCL-Z0(1))*3.281E-3
1242        IF(CBH.LT.3.)THEN
1243          RCBH=.02
1244        ELSE
1245          RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(-            &
1246               1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1247        ENDIF
1248        IF(CBH.GT.25)RCBH=2.4
1249        PEFCBH=1./(1.+RCBH)
1250        PEFCBH=AMIN1(PEFCBH,.9)
1251!     
1252!... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1253!     
1254        PEFF=.5*(PEF+PEFCBH)
1255        PEFF2 = PEFF                                ! JSK MODS
1256       IF(IPRNT)THEN 
1257         WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1258!       call flush(98)   
1259       endif     
1260!        WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1261!*****************************************************************
1262!                                                                *
1263!                  COMPUTE DOWNDRAFT PROPERTIES                  *
1264!                                                                *
1265!*****************************************************************
1266!     
1267!     
1268       TDER=0.
1269 devap:IF(ISHALL.EQ.1)THEN
1270         LFS = 1
1271       ELSE
1272!
1273!...start downdraft about 150 mb above cloud base...
1274!
1275!        KSTART=MAX0(KPBL,KLCL)
1276!        KSTART=KPBL                                  ! Changed 7/23/99
1277         KSTART=KPBL+1                                ! Changed 7/23/99
1278         KLFS = LET-1
1279         DO NK = KSTART+1,KL
1280           DPPP = P0(KSTART)-P0(NK)
1281!          IF(DPPP.GT.200.E2)THEN
1282           IF(DPPP.GT.150.E2)THEN
1283             KLFS = NK
1284             EXIT
1285           ENDIF
1286         ENDDO
1287         KLFS = MIN0(KLFS,LET-1)
1288         LFS = KLFS
1289!
1290!...if LFS is not at least 50 mb above cloud base (implying that the
1291!...level of equil temp, LET, is just above cloud base) do not allow a
1292!...downdraft...
1293!
1294        IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
1295          THETED(LFS) = THETEE(LFS)
1296          QD(LFS) = Q0(LFS)
1297!
1298!...call tpmix2dd to find wet-bulb temp, qv...
1299!
1300          call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
1301          THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
1302!     
1303!...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
1304!     
1305          TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
1306          RDD=P0(LFS)/(R*TVD(LFS))
1307          A1=(1.-PEFF)*AU0
1308          DMF(LFS)=-A1*RDD
1309          DER(LFS)=DMF(LFS)
1310          DDR(LFS)=0.
1311          RHBAR = RH(LFS)*DP(LFS)
1312          DPTT = DP(LFS)
1313          DO ND = LFS-1,KSTART,-1
1314            ND1 = ND+1
1315            DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
1316            DDR(ND)=0.
1317            DMF(ND)=DMF(ND1)+DER(ND)
1318            THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1319            QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)   
1320            DPTT = DPTT+DP(ND)
1321            RHBAR = RHBAR+RH(ND)*DP(ND)
1322          ENDDO
1323          RHBAR = RHBAR/DPTT
1324          DMFFRC = 2.*(1.-RHBAR)
1325          DPDD = 0.
1326!...Calculate melting effect
1327!... first, compute total frozen precipitation generated...
1328!
1329          pptmlt = 0.
1330          DO NK = KLCL,LTOP
1331            PPTMLT = PPTMLT+PPTICE(NK)
1332          ENDDO
1333          if(lc.lt.ml)then
1334!...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
1335!...if DMFFRC=1.  Otherwise, for small DMFFRC, DTMELT gets too large!
1336!...12/14/98 jsk...
1337            DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
1338          else
1339            DTMELT = 0.
1340          endif
1341          LDT = MIN0(LFS-1,KSTART-1)
1342!
1343          call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
1344!
1345          tz(kstart) = tz(kstart)-dtmelt
1346          ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
1347          QSS=0.622*ES/(P0(KSTART)-ES)
1348          THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))*    &
1349                EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
1350!.... 
1351          LDT = MIN0(LFS-1,KSTART-1)
1352          DO ND = LDT,1,-1
1353            DPDD = DPDD+DP(ND)
1354            THETED(ND) = THETED(KSTART)
1355            QD(ND)     = QD(KSTART)       
1356!
1357!...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
1358!
1359            call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
1360            qsd(nd) = qss
1361!
1362!...specify RH decrease of 20%/km in downdraft...
1363!
1364            RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
1365!
1366!...adjust downdraft TEMP, Q to specified RH:
1367!
1368            IF(RHH.LT.1.)THEN
1369              DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1370              RL=XLV0-XLV1*TZ(ND)
1371              DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
1372              T1RH=TZ(ND)+DTMP
1373              ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1374              QSRH=0.622*ES/(P0(ND)-ES)
1375!
1376!...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1377!...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1378!
1379              IF(QSRH.LT.QD(ND))THEN
1380                QSRH=QD(ND)
1381                T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
1382              ENDIF
1383              TZ(ND)=T1RH
1384              QSS=QSRH
1385              QSD(ND) = QSS
1386            ENDIF         
1387            TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
1388            IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
1389              LDB=ND
1390              EXIT
1391            ENDIF
1392          ENDDO
1393          IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN   ! minimum Downdraft depth!
1394            DO ND=LDT,LDB,-1
1395              ND1 = ND+1
1396              DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
1397              DER(ND) = 0.
1398              DMF(ND) = DMF(ND1)+DDR(ND)
1399              TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
1400              QD(ND)=QSD(nd)
1401              THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1402            ENDDO
1403          ENDIF
1404        ENDIF
1405      ENDIF devap
1406!
1407!...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1408!...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1409!
1410d_mf:   IF(TDER.LT.1.)THEN
1411!           WRITE(98,3004)I,J
1412!3004       FORMAT(' ','No Downdraft!;  I=',I3,2X,'J=',I3,'ISHALL =',I2)
1413          PPTFLX=TRPPT
1414          CPR=TRPPT
1415          TDER=0.
1416          CNDTNF=0.
1417          UPDINC=1.
1418          LDB=LFS
1419          DO NDK=1,LTOP
1420            DMF(NDK)=0.
1421            DER(NDK)=0.
1422            DDR(NDK)=0.
1423            THTAD(NDK)=0.
1424            WD(NDK)=0.
1425            TZ(NDK)=0.
1426            QD(NDK)=0.
1427          ENDDO
1428          AINCM2=100.
1429        ELSE
1430          DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
1431          UPDINC=1.
1432          IF(TDER*DDINC.GT.TRPPT)THEN
1433            DDINC = TRPPT/TDER
1434          ENDIF
1435          TDER = TDER*DDINC
1436          DO NK=LDB,LFS
1437            DMF(NK)=DMF(NK)*DDINC
1438            DER(NK)=DER(NK)*DDINC
1439            DDR(NK)=DDR(NK)*DDINC
1440          ENDDO
1441         CPR=TRPPT
1442         PPTFLX = TRPPT-TDER
1443         PEFF=PPTFLX/TRPPT
1444         IF(IPRNT)THEN
1445           write(98,*)'PRECIP EFFICIENCY =',PEFF
1446!          call flush(98)   
1447         ENDIF
1448!
1449!
1450!...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1451!   DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1452!   FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1453!     
1454!         DO NK=LC,LFS
1455!           UMF(NK)=UMF(NK)*UPDINC
1456!           UDR(NK)=UDR(NK)*UPDINC
1457!           UER(NK)=UER(NK)*UPDINC
1458!           PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1459!           PPTICE(NK)=PPTICE(NK)*UPDINC
1460!           DETLQ(NK)=DETLQ(NK)*UPDINC
1461!           DETIC(NK)=DETIC(NK)*UPDINC
1462!         ENDDO
1463!     
1464!...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1465!...DOWNDRAFT...
1466!     
1467         IF(LDB.GT.1)THEN
1468           DO NK=1,LDB-1
1469             DMF(NK)=0.
1470             DER(NK)=0.
1471             DDR(NK)=0.
1472             WD(NK)=0.
1473             TZ(NK)=0.
1474             QD(NK)=0.
1475             THTAD(NK)=0.
1476           ENDDO
1477         ENDIF
1478         DO NK=LFS+1,KX
1479           DMF(NK)=0.
1480           DER(NK)=0.
1481           DDR(NK)=0.
1482           WD(NK)=0.
1483           TZ(NK)=0.
1484           QD(NK)=0.
1485           THTAD(NK)=0.
1486         ENDDO
1487         DO NK=LDT+1,LFS-1
1488           TZ(NK)=0.
1489           QD(NK)=0.
1490           THTAD(NK)=0.
1491         ENDDO
1492       ENDIF d_mf
1493!
1494!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFL
1495!   INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILAB
1496!   IN THAT LAYER INITIALLY...
1497!     
1498       AINCMX=1000.
1499       LMAX=MAX0(KLCL,LFS)
1500       DO NK=LC,LMAX
1501         IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
1502           AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
1503           AINCMX=AMIN1(AINCMX,AINCM1)
1504         ENDIF
1505       ENDDO
1506       AINC=1.
1507       IF(AINCMX.LT.AINC)AINC=AINCMX
1508!     
1509!...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
1510!...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
1511!...CLOSURE...
1512!     
1513       TDER2=TDER
1514       PPTFL2=PPTFLX
1515       DO NK=1,LTOP
1516         DETLQ2(NK)=DETLQ(NK)
1517         DETIC2(NK)=DETIC(NK)
1518         UDR2(NK)=UDR(NK)
1519         UER2(NK)=UER(NK)
1520         DDR2(NK)=DDR(NK)
1521         DER2(NK)=DER(NK)
1522         UMF2(NK)=UMF(NK)
1523         DMF2(NK)=DMF(NK)
1524       ENDDO
1525       FABE=1.
1526       STAB=0.95
1527       NOITR=0
1528       ISTOP=0
1529!
1530        IF(ISHALL.EQ.1)THEN                              ! First for shallow convection
1531!
1532! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
1533! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
1534! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
1535!
1536!...find the maximum TKE value between LC and KLCL...
1537!         TKEMAX = 0.
1538          TKEMAX = 5.
1539!          DO 173 K = LC,KLCL
1540!            NK = KX-K+1
1541!            TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
1542! 173      CONTINUE
1543!          TKEMAX = AMIN1(TKEMAX,10.)
1544!          TKEMAX = AMAX1(TKEMAX,5.)
1545!c         TKEMAX = 10.
1546!c...3_24_99...DPMIN was changed for shallow convection so that it is the
1547!c...          the same as for deep convection (5.E3).  Since this doubles
1548!c...          (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
1549!c...          lation of EVAC...
1550!c         EVAC  = TKEMAX*0.1
1551          EVAC  = 0.5*TKEMAX*0.1
1552!         AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
1553!          AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
1554          AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
1555          TDER=TDER2*AINC
1556          PPTFLX=PPTFL2*AINC
1557          DO NK=1,LTOP
1558            UMF(NK)=UMF2(NK)*AINC
1559            DMF(NK)=DMF2(NK)*AINC
1560            DETLQ(NK)=DETLQ2(NK)*AINC
1561            DETIC(NK)=DETIC2(NK)*AINC
1562            UDR(NK)=UDR2(NK)*AINC
1563            UER(NK)=UER2(NK)*AINC
1564            DER(NK)=DER2(NK)*AINC
1565            DDR(NK)=DDR2(NK)*AINC
1566          ENDDO
1567        ENDIF                                           ! Otherwise for deep convection
1568! use iterative procedure to find mass fluxes...
1569iter:     DO NCOUNT=1,10
1570!     
1571!*****************************************************************
1572!                                                                *
1573!           COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE     *
1574!                                                                *
1575!*****************************************************************
1576!     
1577!...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1578!...SATISFY MASS CONTINUITY...
1579!     
1580            DTT=TIMEC
1581            DO NK=1,LTOP
1582              DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1583              IF(NK.GT.1)THEN
1584                OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1585                ABSOMG = ABS(OMG(NK))
1586                ABSOMGTC = ABSOMG*TIMEC
1587                FRDP = 0.75*DP(NK-1)
1588                IF(ABSOMGTC.GT.FRDP)THEN
1589                  DTT1 = FRDP/ABSOMG
1590                  DTT=AMIN1(DTT,DTT1)
1591                ENDIF
1592              ENDIF
1593            ENDDO
1594            DO NK=1,LTOP
1595              THPA(NK)=THTA0(NK)
1596              QPA(NK)=Q0(NK)
1597              NSTEP=NINT(TIMEC/DTT+1)
1598              DTIME=TIMEC/FLOAT(NSTEP)
1599              FXM(NK)=OMG(NK)*DXSQ/G
1600            ENDDO
1601!     
1602!...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1603!     
1604        DO NTC=1,NSTEP
1605!     
1606!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1607!...SIGN OF OMEGA...
1608!     
1609            DO  NK=1,LTOP
1610              THFXIN(NK)=0.
1611              THFXOUT(NK)=0.
1612              QFXIN(NK)=0.
1613              QFXOUT(NK)=0.
1614            ENDDO
1615            DO NK=2,LTOP
1616              IF(OMG(NK).LE.0.)THEN
1617                THFXIN(NK)=-FXM(NK)*THPA(NK-1)
1618                QFXIN(NK)=-FXM(NK)*QPA(NK-1)
1619                THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
1620                QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
1621              ELSE
1622                THFXOUT(NK)=FXM(NK)*THPA(NK)
1623                QFXOUT(NK)=FXM(NK)*QPA(NK)
1624                THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
1625                QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
1626              ENDIF
1627            ENDDO
1628!     
1629!...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
1630!     
1631            DO NK=1,LTOP
1632              THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)*      &
1633                       THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))*  &
1634                       DTIME*EMSD(NK)
1635              QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)-    &
1636                      QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1637            ENDDO   
1638          ENDDO   
1639          DO NK=1,LTOP
1640            THTAG(NK)=THPA(NK)
1641            QG(NK)=QPA(NK)
1642          ENDDO
1643!     
1644!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE;  IF SO, BORRO
1645!...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
1646!     
1647        DO NK=1,LTOP
1648          IF(QG(NK).LT.0.)THEN
1649            IF(NK.EQ.1)THEN                             ! JSK MODS
1650!              PRINT *,' PROBLEM WITH KF SCHEME:  ' ! JSK MODS
1651!              PRINT *,'QG = 0 AT THE SURFACE!!!!!!!'    ! JSK MODS
1652              CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
1653            ENDIF                                       ! JSK MODS
1654            NK1=NK+1
1655            IF(NK.EQ.LTOP)THEN
1656              NK1=KLCL
1657            ENDIF
1658            TMA=QG(NK1)*EMS(NK1)
1659            TMB=QG(NK-1)*EMS(NK-1)
1660            TMM=(QG(NK)-1.E-9)*EMS(NK  )
1661            BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1662            ACOEFF=BCOEFF*TMA/TMB
1663            TMB=TMB*(1.-BCOEFF)
1664            TMA=TMA*(1.-ACOEFF)
1665            IF(NK.EQ.LTOP)THEN
1666              QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1667!              IF(ABS(QVDIFF).GT.1.)THEN
1668!             PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ',     &
1669!                      QVDIFF,                                                &
1670!                     '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ',     &
1671!                     'VALUES IN KAIN-FRITSCH'
1672!              ENDIF
1673            ENDIF
1674            QG(NK)=1.E-9
1675            QG(NK1)=TMA*EMSD(NK1)
1676            QG(NK-1)=TMB*EMSD(NK-1)
1677          ENDIF
1678        ENDDO
1679        TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1680        IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1681!       WRITE(99,*)'ERROR:  MASS DOES NOT BALANCE IN KF SCHEME;            &
1682!      TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1683!      TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1684          ISTOP=1
1685          IPRNT=.TRUE.
1686          EXIT iter
1687        ENDIF
1688!     
1689!...CONVERT THETA TO T...
1690!     
1691        DO NK=1,LTOP
1692          EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1693          TG(NK)=THTAG(NK)/EXN(NK)
1694          TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1695        ENDDO
1696        IF(ISHALL.EQ.1)THEN
1697          EXIT iter
1698        ENDIF
1699!     
1700!*******************************************************************
1701!                                                                  *
1702!     COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY.    *
1703!                                                                  *
1704!*******************************************************************
1705!     
1706!...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1707!     
1708!        THMIX=0.
1709          TMIX=0.
1710          QMIX=0.
1711!
1712!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
1713!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
1714!...LAYERS...
1715!
1716          DO NK=LC,KPBL
1717            TMIX=TMIX+DP(NK)*TG(NK)
1718            QMIX=QMIX+DP(NK)*QG(NK) 
1719          ENDDO
1720          TMIX=TMIX/DPTHMX
1721          QMIX=QMIX/DPTHMX
1722          ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
1723          QSS=0.622*ES/(PMIX-ES)
1724!     
1725!...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1726!     
1727          IF(QMIX.GT.QSS)THEN
1728            RL=XLV0-XLV1*TMIX
1729            CPM=CP*(1.+0.887*QMIX)
1730            DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1731            DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
1732            TMIX=TMIX+RL/CP*DQ
1733            QMIX=QMIX-DQ
1734            TLCL=TMIX
1735          ELSE
1736            QMIX=AMAX1(QMIX,0.)
1737            EMIX=QMIX*PMIX/(0.622+QMIX)
1738            astrt=1.e-3
1739            binc=0.075
1740            a1=emix/aliq
1741            tp=(a1-astrt)/binc
1742            indlu=int(tp)+1
1743            value=(indlu-1)*binc+astrt
1744            aintrp=(a1-value)/binc
1745            tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
1746            TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1747            TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
1748            TLCL=AMIN1(TLCL,TMIX)
1749          ENDIF
1750          TVLCL=TLCL*(1.+0.608*QMIX)
1751          ZLCL = ZMIX+(TLCL-TMIX)/GDRY
1752          DO NK = LC,KL
1753            KLCL=NK
1754            IF(ZLCL.LE.Z0(NK))THEN
1755              EXIT
1756            ENDIF
1757          ENDDO
1758          K=KLCL-1
1759          DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
1760!     
1761!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1762!     
1763          TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1764          QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1765          TVEN=TENV*(1.+0.608*QENV)
1766          PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1767          THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))*             &
1768                  EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
1769!     
1770!...COMPUTE ADJUSTED ABE(ABEG).
1771!     
1772          ABEG=0.
1773          DO NK=K,LTOPM1
1774            NK1=NK+1
1775            THETEU(NK1) = THETEU(NK)
1776!
1777            call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
1778!
1779            TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
1780            IF(NK.EQ.K)THEN
1781              DZZ=Z0(KLCL)-ZLCL
1782              DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
1783            ELSE
1784              DZZ=DZA(NK)
1785              DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
1786            ENDIF
1787            IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
1788!
1789!...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
1790!
1791            CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1792            THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
1793          ENDDO
1794!     
1795!...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1796!...THE PERIOD TIMEC...
1797!     
1798          IF(NOITR.EQ.1)THEN
1799!         write(98,*)' '
1800!         write(98,*)'TAU, I, J, =',NTSD,I,J
1801!         WRITE(98,1060)FABE
1802!          GOTO 265
1803          EXIT iter
1804          ENDIF
1805          DABE=AMAX1(ABE-ABEG,0.1*ABE)
1806          FABE=ABEG/ABE
1807          IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
1808!          WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
1809!     *GRID POINT; NO CONVECTION ALLOWED!'
1810            RETURN 
1811          ENDIF
1812          IF(NCOUNT.NE.1)THEN
1813            IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
1814              NOITR=1
1815              AINC=AINCOLD
1816              CYCLE iter
1817            ENDIF
1818            DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1819            IF(DFDA.GT.0.)THEN
1820              NOITR=1
1821              AINC=AINCOLD
1822              CYCLE iter
1823            ENDIF
1824          ENDIF
1825          AINCOLD=AINC
1826          FABEOLD=FABE
1827          IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1828!           write(98,*)' '
1829!           write(98,*)'TAU, I, J, =',NTSD,I,J
1830!           WRITE(98,1055)FABE
1831!            GOTO 265
1832            EXIT
1833          ENDIF
1834          IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
1835            EXIT iter
1836          ELSE
1837            IF(NCOUNT.GT.10)THEN
1838!             write(98,*)' '
1839!             write(98,*)'TAU, I, J, =',NTSD,I,J
1840!             WRITE(98,1060)FABE
1841!             GOTO 265
1842              EXIT
1843            ENDIF
1844!     
1845!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTI
1846!...MASS FLUX BY THE FACTOR AINC:
1847!     
1848            IF(FABE.EQ.0.)THEN
1849              AINC=AINC*0.5
1850            ELSE
1851              IF(DABE.LT.1.e-4)THEN
1852                NOITR=1
1853                AINC=AINCOLD
1854                CYCLE iter
1855              ELSE
1856                AINC=AINC*STAB*ABE/DABE
1857              ENDIF
1858            ENDIF
1859!           AINC=AMIN1(AINCMX,AINC)
1860            AINC=AMIN1(AINCMX,AINC)
1861!...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
1862!...WILL BE MINIMAL SO JUST IGNORE IT...              ! JSK MODS
1863            IF(AINC.LT.0.05)then
1864              RETURN                          ! JSK MODS
1865            ENDIF
1866!            AINC=AMAX1(AINC,0.05)                        ! JSK MODS
1867            TDER=TDER2*AINC
1868            PPTFLX=PPTFL2*AINC
1869!           IF (XTIME.LT.10.)THEN
1870!           WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
1871!          *              FABEOLD,AINCOLD
1872!           ENDIF
1873            DO NK=1,LTOP
1874              UMF(NK)=UMF2(NK)*AINC
1875              DMF(NK)=DMF2(NK)*AINC
1876              DETLQ(NK)=DETLQ2(NK)*AINC
1877              DETIC(NK)=DETIC2(NK)*AINC
1878              UDR(NK)=UDR2(NK)*AINC
1879              UER(NK)=UER2(NK)*AINC
1880              DER(NK)=DER2(NK)*AINC
1881              DDR(NK)=DDR2(NK)*AINC
1882            ENDDO
1883!     
1884!...GO BACK UP FOR ANOTHER ITERATION...
1885!     
1886          ENDIF
1887        ENDDO iter
1888!     
1889!...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1890!     
1891!...FRC2 IS THE FRACTION OF TOTAL CONDENSATE      !  PPT FB MODS
1892!...GENERATED THAT GOES INTO PRECIPITIATION       !  PPT FB MODS
1893!
1894!  Redistribute hydormeteors according to the final mass-flux values:
1895!
1896        IF(CPR.GT.0.)THEN
1897          FRC2=PPTFLX/(CPR*AINC)                    !  PPT FB MODS
1898        ELSE
1899           FRC2=0.
1900        ENDIF
1901        DO NK=1,LTOP
1902          QLPA(NK)=QL0(NK)
1903          QIPA(NK)=QI0(NK)
1904          QRPA(NK)=QR0(NK)
1905          QSPA(NK)=QS0(NK)
1906          RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2   !  PPT FB MODS
1907          SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2   !  PPT FB MODS
1908        ENDDO
1909        DO NTC=1,NSTEP
1910!     
1911!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAY
1912!...BASED ON THE SIGN OF OMEGA...
1913!     
1914          DO NK=1,LTOP
1915            QLFXIN(NK)=0.
1916            QLFXOUT(NK)=0.
1917            QIFXIN(NK)=0.
1918            QIFXOUT(NK)=0.
1919            QRFXIN(NK)=0.
1920            QRFXOUT(NK)=0.
1921            QSFXIN(NK)=0.
1922            QSFXOUT(NK)=0.
1923          ENDDO   
1924          DO NK=2,LTOP
1925            IF(OMG(NK).LE.0.)THEN
1926              QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1927              QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1928              QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1929              QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1930              QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1931              QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1932              QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1933              QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1934            ELSE
1935              QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1936              QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1937              QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1938              QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1939              QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1940              QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1941              QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1942              QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1943            ENDIF
1944          ENDDO   
1945!     
1946!...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1947!     
1948          DO NK=1,LTOP
1949            QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
1950            QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
1951            QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK)         !  PPT FB MODS
1952            QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK)         !  PPT FB MODS
1953          ENDDO     
1954        ENDDO
1955        DO NK=1,LTOP
1956          QLG(NK)=QLPA(NK)
1957          QIG(NK)=QIPA(NK)
1958          QRG(NK)=QRPA(NK)
1959          QSG(NK)=QSPA(NK)
1960        ENDDO   
1961!
1962!...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1963!...GRID POINT...
1964!     
1965!     IF (XTIME.LT.10.)THEN
1966!     WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1967!     ENDIF
1968       IF(IPRNT)THEN 
1969         WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1970!        call flush(98)   
1971       endif 
1972!     
1973!...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1974!     
1975!297   IF(IPRNT)then
1976       IF(IPRNT)then
1977!    if(I.eq.16 .and. J.eq.41)then
1978!      IF(ISTOP.EQ.1)THEN
1979         write(98,*)
1980!        write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
1981         write(98,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100.,       &
1982                     TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
1983         write(98,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL,       &
1984                      DTRH,TENV   
1985         WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,       &
1986         TMIX-T00,PMIX,QMIX,ABE
1987         WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,  &
1988         WLCL,CLDHGT(LC)
1989         WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1990         write(98,*)'PRECIP EFFICIENCY =',PEFF
1991      WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1992!      ENDIF
1993!!!!! HERE !!!!!!!
1994           WRITE(98,1070)'  P  ','   DP ',' DT K/D ',' DR K/D ','   OMG  ',        &
1995          ' DOMGDP ','   UMF  ','   UER  ','   UDR  ','   DMF  ','   DER  '        &
1996          ,'   DDR  ','   EMS  ','    W0  ','  DETLQ ',' DETIC '
1997           write(98,*)'just before DO 300...'
1998!          call flush(98)
1999           DO NK=1,LTOP
2000             K=LTOP-NK+1
2001             DTT=(TG(K)-T0(K))*86400./TIMEC
2002             RL=XLV0-XLV1*TG(K)
2003             DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
2004             UDFRC=UDR(K)*TIMEC*EMSD(K)
2005             UEFRC=UER(K)*TIMEC*EMSD(K)
2006             DDFRC=DDR(K)*TIMEC*EMSD(K)
2007             DEFRC=-DER(K)*TIMEC*EMSD(K)
2008             WRITE(98,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4,       &
2009             UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11,           &
2010             W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)*                   &
2011             TIMEC*EMSD(K)*1.E3
2012           ENDDO
2013           WRITE(98,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG',             &
2014                  'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
2015           DO NK=1,KL
2016             K=KX-NK+1
2017             DTT=TG(K)-T0(K)
2018             TUC=TU(K)-T00
2019             IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
2020             TDC=TZ(K)-T00
2021             IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
2022             IF(T0(K).LT.T00)THEN
2023               ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2024             ELSE
2025               ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2026             ENDIF 
2027             QGS=ES*0.622/(P0(K)-ES)
2028             RH0=Q0(K)/QES(K)
2029             RHG=QG(K)/QGS
2030             WRITE(98,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC,            &
2031             TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)*                   &
2032             1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000.,                &
2033             QSG(K)*1000.,RH0,RHG
2034           ENDDO
2035!     
2036!...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
2037!...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2038!     
2039!         IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
2040
2041!         IF(ISHALL.NE.1)THEN
2042!            write(98,4421)i,j,iyr,imo,idy,ihr,imn
2043!           write(98)i,j,iyr,imo,idy,ihr,imn,kl
2044! 4421       format(7i4)
2045!            write(98,4422)kl
2046! 4422       format(i6)
2047            DO 310 NK = 1,KL
2048              k = kl - nk + 1
2049              write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000.,       &
2050                       u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
2051!             write(98) p0,t0,q0,u0,v0,w0,dp,tke
2052!           WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
2053!    *               U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
2054 310        CONTINUE
2055            IF(ISTOP.EQ.1)THEN
2056              CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
2057            ENDIF
2058!         ENDIF
2059  4455  format(8f11.3)
2060       ENDIF
2061        CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
2062        RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ     !  PPT FB MODS
2063!        RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ               !  PPT FB MODS
2064!         RNC=0.1*TIMEC*PPTFLX/DXSQ
2065        RNC=RAINCV(I,J)*NIC
2066       IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
2067
2068!     WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2069!     
2070!  EVALUATE MOISTURE BUDGET...
2071!     
2072
2073        QINIT=0.
2074        QFNL=0.
2075        DPT=0.
2076        DO 315 NK=1,LTOP
2077          DPT=DPT+DP(NK)
2078          QINIT=QINIT+Q0(NK)*EMS(NK)
2079          QFNL=QFNL+QG(NK)*EMS(NK)
2080          QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
2081  315   CONTINUE
2082        QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)       !  PPT FB MODS
2083!        QFNL=QFNL+PPTFLX*TIMEC                 !  PPT FB MODS
2084        ERR2=(QFNL-QINIT)*100./QINIT
2085       IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
2086      IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
2087!       write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
2088!       WRITE(99,1110)QINIT,QFNL,ERR2
2089        IPRNT=.TRUE.
2090        ISTOP=1
2091            write(98,4422)kl
2092 4422       format(i6)
2093            DO 311 NK = 1,KL
2094              k = kl - nk + 1
2095!             write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000.,       &
2096!                      u0(k),v0(k),W0AVG1D(K),dp(k)
2097!             write(98) p0,t0,q0,u0,v0,w0,dp,tke
2098!           WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,          &
2099!                    U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2100            WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,          &
2101                     U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2102 311        CONTINUE
2103!           call flush(98)
2104
2105!        GOTO 297
2106!         STOP 'QVERR'
2107      ENDIF
2108 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2109 4456  format(8f12.3)
2110        IF(PPTFLX.GT.0.)THEN
2111          RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
2112        ELSE
2113          RELERR=0.
2114        ENDIF
2115     IF(IPRNT)THEN
2116        WRITE(98,1120)RELERR
2117        WRITE(98,*)'TDER, CPR, TRPPT =',              &
2118          TDER,CPR*AINC,TRPPT*AINC
2119     ENDIF
2120!     
2121!...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
2122!     
2123!...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
2124!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2125!     
2126        IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
2127        NCA(I,J)=FLOAT(NIC)
2128        IF(ISHALL.EQ.1)THEN
2129          TIMEC = 2400.
2130          NCA(I,J) = FLOAT(NTST)
2131          NSHALL = NSHALL+1
2132        ENDIF
2133        DO K=1,KX
2134!         IF(IMOIST(INEST).NE.2)THEN
2135!
2136!...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
2137!...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
2138!...NOTE:  THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
2139!...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
2140!...OF QG...
2141!
2142!           RLC=XLV0-XLV1*TG(K)
2143!           RLS=XLS0-XLS1*TG(K)
2144!           CPM=CP*(1.+0.887*QG(K))
2145!           TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
2146!           QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
2147!           DQLDT(I,J,NK)=0.
2148!           DQIDT(I,J,NK)=0.
2149!           DQRDT(I,J,NK)=0.
2150!           DQSDT(I,J,NK)=0.
2151!         ELSE
2152!
2153!...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2154!
2155          IF(.NOT. F_QI .and. warm_rain)THEN
2156
2157            CPM=CP*(1.+0.887*QG(K))
2158            TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2159            DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2160            DQIDT(K)=0.
2161            DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2162            DQSDT(K)=0.
2163          ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
2164!
2165!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
2166!...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2167!
2168            CPM=CP*(1.+0.887*QG(K))
2169            IF(K.LE.ML)THEN
2170              TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2171            ELSEIF(K.GT.ML)THEN
2172              TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2173            ENDIF
2174            DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2175            DQIDT(K)=0.
2176            DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2177            DQSDT(K)=0.
2178          ELSEIF(F_QI) THEN
2179!
2180!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDEN
2181!...OF HYDROMETEORS DIRECTLY...
2182!
2183            DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2184            DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2185            DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2186            IF (F_QS) THEN
2187               DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2188            ELSE
2189               DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2190            ENDIF
2191          ELSE
2192!              PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
2193              CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
2194          ENDIF
2195          DTDT(K)=(TG(K)-T0(K))/TIMEC
2196          DQDT(K)=(QG(K)-Q0(K))/TIMEC
2197        ENDDO
2198        RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ     !  PPT FB MODS
2199!        RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ               !  PPT FB MODS
2200!         RNC=0.1*TIMEC*PPTFLX/DXSQ
2201        RNC=RAINCV(I,J)*NIC
2202 909     FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
2203!      write (98,909)I,J,RNC
2204!      write (6,909)I,J,RNC
2205!      WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
2206!     *            NCCNT
2207!      call flush(98)
22081000  FORMAT(' ',10A8)
22091005  FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
22101010  FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
22111015   FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
22121025   FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M',                         &
2213        ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=',   &
2214        I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1,          &
2215        ' CAPE=',0PF7.1)
22161030   FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =',   &
2217      E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =',                  &
2218      F8.1)
22191035  FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL='       &
2220      ,F6.3,'VWS=',F5.2)
2221!1055  FORMAT('*** DEGREE OF STABILIZATION =',F5.3,                  &
2222!      ', NO MORE MASS FLUX IS ALLOWED!')
2223!1060     FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED    &
2224!      &DEGREE OF STABILIZATION!  FABE= ',F6.4)
2225 1070 FORMAT (16A8)
2226 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2227 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=',           &
2228              2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2229 1085 FORMAT (A3,16A7,2A8)
2230 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2231 1095 FORMAT(' ','  PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
22321105   FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
2233       E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
22341110   FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5,       &
2235       ' TOTAL WATER CHANGE =',F8.2,'%')
2236! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
22371120   FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2238!
2239!-----------------------------------------------------------------------
2240!--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
2241!-----------------------------------------------------------------------
2242!
2243      CUTOP(I,J)=REAL(LTOP)
2244      CUBOT(I,J)=REAL(LCL)
2245!
2246!-----------------------------------------------------------------------
2247   END SUBROUTINE  KF_eta_PARA
2248!********************************************************************
2249! ***********************************************************************
2250   SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
2251!
2252! Lookup table variables:
2253!     INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2254!     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2255!     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2256!     REAL, SAVE, DIMENSION(1:200) :: ALU
2257!     REAL, SAVE :: RDPR,RDTHK,PLUTOP
2258! End of Lookup table variables:
2259!-----------------------------------------------------------------------
2260   IMPLICIT NONE
2261!-----------------------------------------------------------------------
2262   REAL,         INTENT(IN   )   :: P,THES,XLV1,XLV0
2263   REAL,         INTENT(OUT  )   :: QNEWLQ,QNEWIC
2264   REAL,         INTENT(INOUT)   :: TU,QU,QLIQ,QICE
2265   REAL    ::    TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11,          &
2266                 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
2267   INTEGER ::    IPTB,ITHTB
2268!-----------------------------------------------------------------------
2269
2270!c******** LOOKUP TABLE VARIABLES... ****************************
2271!      parameter(kfnt=250,kfnp=220)
2272!c
2273!      COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
2274!     *              alu(200),rdpr,rdthk,plutop
2275!C***************************************************************
2276!c
2277!c***********************************************************************
2278!c     scaling pressure and tt table index                         
2279!c***********************************************************************
2280!c
2281      tp=(p-plutop)*rdpr
2282      qq=tp-aint(tp)
2283      iptb=int(tp)+1
2284
2285!
2286!***********************************************************************
2287!              base and scaling factor for the                           
2288!***********************************************************************
2289!
2290!  scaling the and tt table index                                       
2291      bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2292      tth=(thes-bth)*rdthk
2293      pp   =tth-aint(tth)
2294      ithtb=int(tth)+1
2295       IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
2296         write(98,*)'**** OUT OF BOUNDS *********'
2297!        call flush(98)
2298       ENDIF
2299!
2300      t00=ttab(ithtb  ,iptb  )
2301      t10=ttab(ithtb+1,iptb  )
2302      t01=ttab(ithtb  ,iptb+1)
2303      t11=ttab(ithtb+1,iptb+1)
2304!
2305      q00=qstab(ithtb  ,iptb  )
2306      q10=qstab(ithtb+1,iptb  )
2307      q01=qstab(ithtb  ,iptb+1)
2308      q11=qstab(ithtb+1,iptb+1)
2309!
2310!***********************************************************************
2311!              parcel temperature                                       
2312!***********************************************************************
2313!
2314      temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2315!
2316      qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2317!
2318      DQ=QS-QU
2319      IF(DQ.LE.0.)THEN
2320        QNEW=QU-QS
2321        QU=QS
2322      ELSE
2323!
2324!   IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2325!   ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
2326!
2327        QNEW=0.
2328        QTOT=QLIQ+QICE
2329!
2330!   IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
2331!   WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
2332!   RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
2333!   DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
2334!   ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
2335!
2336!...subsaturated values only occur in calculations involving various mixtures of
2337!...updraft and environmental air for estimation of entrainment and detrainment.
2338!...For these purposes, assume that reasonable estimates can be given using
2339!...liquid water saturation calculations only - i.e., ignore the effect of the
2340!...ice phase in this process only...will not affect conservative properties...
2341!
2342        IF(QTOT.GE.DQ)THEN
2343          qliq=qliq-dq*qliq/(qtot+1.e-10)
2344          qice=qice-dq*qice/(qtot+1.e-10)
2345          QU=QS
2346        ELSE
2347          RLL=XLV0-XLV1*TEMP
2348          CPP=1004.5*(1.+0.89*QU)
2349          IF(QTOT.LT.1.E-10)THEN
2350!
2351!...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2352            TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
2353          ELSE
2354!
2355!...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
2356!   THE TEMPERATURE IS GIVEN BY:
2357!
2358            TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
2359            QU=QU+QTOT
2360            QTOT=0.
2361            QLIQ=0.
2362            QICE=0.
2363          ENDIF
2364        ENDIF
2365      ENDIF
2366      TU=TEMP
2367      qnewlq=qnew
2368      qnewic=0.
2369!
2370   END SUBROUTINE TPMIX2
2371!******************************************************************************
2372      SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
2373!-----------------------------------------------------------------------
2374   IMPLICIT NONE
2375!-----------------------------------------------------------------------
2376   REAL,         INTENT(IN   )   :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
2377   REAL,         INTENT(INOUT)   :: TU,THTEU,QU,QICE
2378   REAL    ::    RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
2379!-----------------------------------------------------------------------
2380!
2381!...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
2382!...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
2383!...TTFRZ TO TBFRZ...
2384!...FOR COLDER TERMPERATURES, FREEZE ALL LIQUID WATER...
2385!...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
2386!...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
2387!
2388      RLC=2.5E6-2369.276*(TU-273.16)
2389      RLS=2833922.-259.532*(TU-273.16)
2390      RLF=RLS-RLC
2391      CPP=1004.5*(1.+0.89*QU)
2392!
2393!  A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
2394!  FOR SATURATION VAPOR PRESSURE...
2395!
2396      A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
2397      DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
2398      TU = TU+DTFRZ
2399     
2400      ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2401      QS = ES*0.622/(P-ES)
2402!
2403!...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
2404!...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
2405!...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
2406!...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
2407!...TEMPERATURE TO THE SATURATION VALUE...
2408!
2409      DQEVAP = QS-QU
2410      QICE = QICE-DQEVAP
2411      QU = QU+DQEVAP
2412      PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
2413      THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
2414!
2415   END SUBROUTINE DTFRZNEW
2416! --------------------------------------------------------------------------------
2417
2418      SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ,           &
2419                          QNEWIC,QLQOUT,QICOUT,G)
2420
2421!-----------------------------------------------------------------------
2422   IMPLICIT NONE
2423!-----------------------------------------------------------------------
2424!  9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2425!  BY OGURA AND CHO (1973).  LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2426!  CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2427!  CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2428!  RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2429
2430      REAL, INTENT(IN   )   :: G
2431      REAL, INTENT(IN   )   :: DZ,BOTERM,ENTERM,RATE
2432      REAL, INTENT(INOUT)   :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2433      REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2434
2435!
2436!  9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2437!  BY OGURA AND CHO (1973).  LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2438!  CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-   
2439!  CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL       
2440!  RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).                   
2441      QTOT=QLIQ+QICE                                                   
2442      QNEW=QNEWLQ+QNEWIC                                               
2443!                                                                       
2444!  ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
2445!  BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2446!  LEVELS...                                                           
2447!                                                                       
2448      QEST=0.5*(QTOT+QNEW)                                             
2449      G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5                             
2450      IF(G1.LT.0.0)G1=0.                                               
2451      WAVG=0.5*(SQRT(WTW)+SQRT(G1))                                     
2452      CONV=RATE*DZ/WAVG                                                 
2453!                                                                       
2454!  RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2455!  THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2456!  IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2457!  SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...       
2458!                                                                       
2459      RATIO3=QNEWLQ/(QNEW+1.E-8)                                       
2460!     OLDQ=QTOT                                                         
2461      QTOT=QTOT+0.6*QNEW                                               
2462      OLDQ=QTOT                                                         
2463      RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)                           
2464      QTOT=QTOT*EXP(-CONV)                                             
2465!                                                                       
2466!  DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT 
2467!  PARCEL AT THIS LEVEL...                                             
2468!                                                                       
2469      DQ=OLDQ-QTOT                                                     
2470      QLQOUT=RATIO4*DQ                                                 
2471      QICOUT=(1.-RATIO4)*DQ                                             
2472!                                                                       
2473!  ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2474!  LATE VERTICAL VELOCITY                                               
2475!                                                                       
2476      PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)                                   
2477      WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5                         
2478      IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
2479!                                                                       
2480!  DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2481!  DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...                 
2482!                                                                       
2483      QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW                                 
2484      QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW                       
2485      QNEWLQ=0.                                                         
2486      QNEWIC=0.                                                         
2487
2488   END SUBROUTINE CONDLOAD
2489
2490! ----------------------------------------------------------------------
2491   SUBROUTINE PROF5(EQ,EE,UD)                                       
2492!
2493!***********************************************************************
2494!*****    GAUSSIAN TYPE MIXING PROFILE....******************************
2495!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2496!  THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN 
2497!  DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
2498!  "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
2499!  ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
2500!  MATHEMATICS SERIES.  JUNE, 1964., MAY, 1968.                         
2501!                                     JACK KAIN                         
2502!                                     7/6/89                           
2503!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2504!-----------------------------------------------------------------------
2505   IMPLICIT NONE
2506!-----------------------------------------------------------------------
2507   REAL,         INTENT(IN   )   :: EQ
2508   REAL,         INTENT(INOUT)   :: EE,UD
2509   REAL ::       SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2510
2511      DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676,       &
2512           0.9372980,0.33267,0.166666667,0.202765151/                       
2513      X=(EQ-0.5)/SIGMA                                                 
2514      Y=6.*EQ-3.                                                       
2515      EY=EXP(Y*Y/(-2))                                                 
2516      E45=EXP(-4.5)                                                     
2517      T2=1./(1.+P*ABS(Y))                                               
2518      T1=0.500498                                                       
2519      C1=A1*T1+A2*T1*T1+A3*T1*T1*T1                                     
2520      C2=A1*T2+A2*T2*T2+A3*T2*T2*T2                                     
2521      IF(Y.GE.0.)THEN                                                   
2522        EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2523        UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.-    &
2524           EQ)                                                         
2525      ELSE                                                             
2526        EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.       
2527        UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ*   &
2528           EQ/2.-EQ)                                                   
2529      ENDIF                                                             
2530      EE=EE/FE                                                         
2531      UD=UD/FE                                                         
2532
2533   END SUBROUTINE PROF5
2534
2535! ------------------------------------------------------------------------
2536   SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
2537!
2538! Lookup table variables:
2539!     INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2540!     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2541!     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2542!     REAL, SAVE, DIMENSION(1:200) :: ALU
2543!     REAL, SAVE :: RDPR,RDTHK,PLUTOP
2544! End of Lookup table variables:
2545!-----------------------------------------------------------------------
2546   IMPLICIT NONE
2547!-----------------------------------------------------------------------
2548   REAL,         INTENT(IN   )   :: P,THES
2549   REAL,         INTENT(INOUT)   :: TS,QS
2550   INTEGER,      INTENT(IN   )   :: i,j     ! avail for debugging
2551   REAL    ::    TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
2552   INTEGER ::    IPTB,ITHTB
2553   CHARACTER*256 :: MESS
2554!-----------------------------------------------------------------------
2555
2556!
2557!******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
2558!     parameter(kfnt=250,kfnp=220)
2559!
2560!     COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),        &
2561!                   alu(200),rdpr,rdthk,plutop
2562!***************************************************************
2563!
2564!***********************************************************************
2565!     scaling pressure and tt table index                         
2566!***********************************************************************
2567!
2568      tp=(p-plutop)*rdpr
2569      qq=tp-aint(tp)
2570      iptb=int(tp)+1
2571!
2572!***********************************************************************
2573!              base and scaling factor for the                           
2574!***********************************************************************
2575!
2576!  scaling the and tt table index                                       
2577      bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2578      tth=(thes-bth)*rdthk
2579      pp   =tth-aint(tth)
2580      ithtb=int(tth)+1
2581!
2582      t00=ttab(ithtb  ,iptb  )
2583      t10=ttab(ithtb+1,iptb  )
2584      t01=ttab(ithtb  ,iptb+1)
2585      t11=ttab(ithtb+1,iptb+1)
2586!
2587      q00=qstab(ithtb  ,iptb  )
2588      q10=qstab(ithtb+1,iptb  )
2589      q01=qstab(ithtb  ,iptb+1)
2590      q11=qstab(ithtb+1,iptb+1)
2591!
2592!***********************************************************************
2593!              parcel temperature and saturation mixing ratio                                       
2594!***********************************************************************
2595!
2596      ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2597!
2598      qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2599!
2600   END SUBROUTINE TPMIX2DD
2601
2602! -----------------------------------------------------------------------
2603  SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)                       
2604!
2605!-----------------------------------------------------------------------
2606   IMPLICIT NONE
2607!-----------------------------------------------------------------------
2608   REAL,         INTENT(IN   )   :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
2609   REAL,         INTENT(INOUT)   :: THT1
2610   REAL    ::    EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT,      &
2611                 T00,P00,C1,C2,C3,C4,C5
2612   INTEGER ::    INDLU
2613!-----------------------------------------------------------------------
2614      DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,   &
2615           0.278296,1.0723E-3/                                         
2616!                                                                       
2617!  CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...         
2618!                                                                       
2619! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
2620!
2621      EE=Q1*P1/(0.622+Q1)                                             
2622!     TLOG=ALOG(EE/ALIQ)                                             
2623! ...calculate LOG term using lookup table...
2624!
2625      astrt=1.e-3
2626      ainc=0.075
2627      a1=ee/aliq
2628      tp=(a1-astrt)/ainc
2629      indlu=int(tp)+1
2630      value=(indlu-1)*ainc+astrt
2631      aintrp=(a1-value)/ainc
2632      tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
2633!
2634      TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)                               
2635      TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2636      THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))                         
2637      THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))                     
2638!
2639  END SUBROUTINE ENVIRTHT                                                             
2640! ***********************************************************************
2641!====================================================================
2642   SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN,      &
2643                     RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS,         &
2644                     SVP1,SVP2,SVP3,SVPT0,                          &
2645                     P_FIRST_SCALAR,restart,allowed_to_read,        &
2646                     ids, ide, jds, jde, kds, kde,                  &
2647                     ims, ime, jms, jme, kms, kme,                  &
2648                     its, ite, jts, jte, kts, kte                   )
2649!--------------------------------------------------------------------
2650   IMPLICIT NONE
2651!--------------------------------------------------------------------
2652   LOGICAL , INTENT(IN)           ::  restart,allowed_to_read
2653   INTEGER , INTENT(IN)           ::  ids, ide, jds, jde, kds, kde, &
2654                                      ims, ime, jms, jme, kms, kme, &
2655                                      its, ite, jts, jte, kts, kte
2656   INTEGER , INTENT(IN)           ::  P_QI,P_QS,P_FIRST_SCALAR
2657
2658   REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::       &
2659                                                          RTHCUTEN, &
2660                                                          RQVCUTEN, &
2661                                                          RQCCUTEN, &
2662                                                          RQRCUTEN, &
2663                                                          RQICUTEN, &
2664                                                          RQSCUTEN
2665
2666   REAL ,   DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2667
2668   REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2669
2670   INTEGER :: i, j, k, itf, jtf, ktf
2671   REAL, INTENT(IN)    :: SVP1,SVP2,SVP3,SVPT0
2672
2673   jtf=min0(jte,jde-1)
2674   ktf=min0(kte,kde-1)
2675   itf=min0(ite,ide-1)
2676
2677   IF(.not.restart)THEN
2678
2679      DO j=jts,jtf
2680      DO k=kts,ktf
2681      DO i=its,itf
2682         RTHCUTEN(i,k,j)=0.
2683         RQVCUTEN(i,k,j)=0.
2684         RQCCUTEN(i,k,j)=0.
2685         RQRCUTEN(i,k,j)=0.
2686      ENDDO
2687      ENDDO
2688      ENDDO
2689
2690      IF (P_QI .ge. P_FIRST_SCALAR) THEN
2691         DO j=jts,jtf
2692         DO k=kts,ktf
2693         DO i=its,itf
2694            RQICUTEN(i,k,j)=0.
2695         ENDDO
2696         ENDDO
2697         ENDDO
2698      ENDIF
2699
2700      IF (P_QS .ge. P_FIRST_SCALAR) THEN
2701         DO j=jts,jtf
2702         DO k=kts,ktf
2703         DO i=its,itf
2704            RQSCUTEN(i,k,j)=0.
2705         ENDDO
2706         ENDDO
2707         ENDDO
2708      ENDIF
2709
2710      DO j=jts,jtf
2711      DO i=its,itf
2712         NCA(i,j)=-100.
2713      ENDDO
2714      ENDDO
2715
2716      DO j=jts,jtf
2717      DO k=kts,ktf
2718      DO i=its,itf
2719         W0AVG(i,k,j)=0.
2720      ENDDO
2721      ENDDO
2722      ENDDO
2723
2724   endif
2725 
2726   CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
2727
2728   END SUBROUTINE kf_eta_init
2729
2730!-------------------------------------------------------
2731
2732      subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
2733!
2734!  This subroutine is a lookup table.
2735!  Given a series of series of saturation equivalent potential
2736!  temperatures, the temperature is calculated.
2737!
2738!--------------------------------------------------------------------
2739   IMPLICIT NONE
2740!--------------------------------------------------------------------
2741! Lookup table variables
2742!     INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
2743!     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2744!     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2745!     REAL, SAVE, DIMENSION(1:200) :: ALU
2746!     REAL, SAVE :: RDPR,RDTHK,PLUTOP
2747! End of Lookup table variables
2748
2749     INTEGER :: KP,IT,ITCNT,I
2750     REAL :: DTH,TMIN,TOLER,PBOT,DPR,                               &
2751             TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
2752             ASTRT,AINC,A1,THTGS
2753!    REAL    :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
2754     REAL    :: ALIQ,BLIQ,CLIQ,DLIQ
2755     REAL, INTENT(IN)    :: SVP1,SVP2,SVP3,SVPT0
2756!
2757! equivalent potential temperature increment
2758      data dth/1./
2759! minimum starting temp
2760      data tmin/150./
2761! tolerance for accuracy of temperature
2762      data toler/0.001/
2763! top pressure (pascals)
2764      plutop=5000.0
2765! bottom pressure (pascals)
2766      pbot=110000.0
2767
2768      ALIQ = SVP1*1000.
2769      BLIQ = SVP2
2770      CLIQ = SVP2*SVPT0
2771      DLIQ = SVP3
2772
2773!
2774! compute parameters
2775!
2776! 1._over_(sat. equiv. theta increment)
2777      rdthk=1./dth
2778! pressure increment
2779!
2780      DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
2781!      dpr=(pbot-plutop)/REAL(kfnp-1)
2782! 1._over_(pressure increment)
2783      rdpr=1./dpr
2784! compute the spread of thes
2785!     thespd=dth*(kfnt-1)
2786!
2787! calculate the starting sat. equiv. theta
2788!
2789      temp=tmin
2790      p=plutop-dpr
2791      do kp=1,kfnp
2792        p=p+dpr
2793        es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
2794        qs=0.622*es/(p-es)
2795        pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2796        the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs*        &
2797               (1.+0.81*qs))
2798      enddo   
2799!
2800! compute temperatures for each sat. equiv. potential temp.
2801!
2802      p=plutop-dpr
2803      do kp=1,kfnp
2804        thes=the0k(kp)-dth
2805        p=p+dpr
2806        do it=1,kfnt
2807! define sat. equiv. pot. temp.
2808          thes=thes+dth
2809! iterate to find temperature
2810! find initial guess
2811          if(it.eq.1) then
2812            tgues=tmin
2813          else
2814            tgues=ttab(it-1,kp)
2815          endif
2816          es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
2817          qs=0.622*es/(p-es)
2818          pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2819          thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs*      &
2820               (1.+0.81*qs))
2821          f0=thgues-thes
2822          t1=tgues-0.5*f0
2823          t0=tgues
2824          itcnt=0
2825! iteration loop
2826          do itcnt=1,11
2827            es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
2828            qs=0.622*es/(p-es)
2829            pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2830            thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
2831            f1=thtgs-thes
2832            if(abs(f1).lt.toler)then
2833              exit
2834            endif
2835!           itcnt=itcnt+1
2836            dt=f1*(t1-t0)/(f1-f0)
2837            t0=t1
2838            f0=f1
2839            t1=t1-dt
2840          enddo
2841          ttab(it,kp)=t1
2842          qstab(it,kp)=qs
2843        enddo
2844      enddo   
2845!
2846! lookup table for tlog(emix/aliq)
2847!
2848! set up intial values for lookup tables
2849!
2850       astrt=1.e-3
2851       ainc=0.075
2852!
2853       a1=astrt-ainc
2854       do i=1,200
2855         a1=a1+ainc
2856         alu(i)=alog(a1)
2857       enddo   
2858!
2859   END SUBROUTINE KF_LUTAB
2860
2861END MODULE module_cu_kfeta
Note: See TracBrowser for help on using the repository browser.