phylmd/CV_routines/cv_undilute2.f


      SUBROUTINE cv_undilute2(nloc,ncum,nd,icb,nk &
                             ,tnk,qnk,gznk,t,q,qs,gz &
                             ,p,dph,h,tv,lv &
                             ,inb,inb1,tp,tvp,clw,hp,ep,sigp,frac)
            use cvthermo
            use cv_param
      implicit none

!---------------------------------------------------------------------
! Purpose:
!     FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
!     &
!     COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
!     FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
!     &
!     FIND THE LEVEL OF NEUTRAL BUOYANCY
!---------------------------------------------------------------------


! inputs:
      integer, intent(in):: ncum, nd, nloc
      integer icb(nloc), nk(nloc)
      real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
      real p(nloc,nd), dph(nloc,nd)
      real tnk(nloc), qnk(nloc), gznk(nloc)
      real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)

! outputs:
      integer inb(nloc), inb1(nloc)
      real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
      real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
      real frac(nloc)

! local variables:
      integer i, k
      real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
      real by, defrac
      real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
      logical lcape(nloc)

!=====================================================================
! --- SOME INITIALIZATIONS
!=====================================================================

      do 170 k=1,nl
      do 160 i=1,ncum
       ep(i,k)=0.0
       sigp(i,k)=sigs
 160  continue
 170  continue

!=====================================================================
! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
!=====================================================================
!
! ---       The procedure is to solve the equation.
!              cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
!
!   ***  Calculate certain parcel quantities, including static energy   ***
!
!
      do 240 i=1,ncum
         ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &
               +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
 240  continue
!
!
!    ***  Find lifted parcel quantities above cloud base    ***
!
!
      do 300 k=minorig+1,nl
        do 290 i=1,ncum
          if(k.ge.(icb(i)+1))then
            tg=t(i,k)
            qg=qs(i,k)
            alv=lv0-clmcpv*(t(i,k)-t0)
!
! First iteration.
!
             s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
             s=1./s
             ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
             tg=tg+s*(ah0(i)-ahg)
             tg=max(tg,35.0)
             tc=tg-t0
             denom=243.5+tc
             if(tc.ge.0.0)then
            es=6.112*exp(17.67*tc/denom)
             else
            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
             endif
            qg=eps*es/(p(i,k)-es*(1.-eps))
!
! Second iteration.
!
             s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
             s=1./s
             ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
             tg=tg+s*(ah0(i)-ahg)
             tg=max(tg,35.0)
             tc=tg-t0
             denom=243.5+tc
             if(tc.ge.0.0)then
            es=6.112*exp(17.67*tc/denom)
             else
            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
             endif
            qg=eps*es/(p(i,k)-es*(1.-eps))
!
             alv=lv0-clmcpv*(t(i,k)-t0)
!      print*,'cpd dans convect2 ',cpd
!      print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
!      print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
!              if (.not.cpd.gt.1000.) then
!                  print*,'CPD=',cpd
!                  stop
!              endif
               clw(i,k)=qnk(i)-qg
               clw(i,k)=max(0.0,clw(i,k))
               rg=qg/(1.-qnk(i))
               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
            endif
  290     continue
  300   continue
!
!=====================================================================
! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
!=====================================================================
!
      do 320 k=minorig+1,nl
        do 310 i=1,ncum
          if(k.ge.(nk(i)+1))then
            tca=tp(i,k)-t0
            if(tca.ge.0.0)then
              elacrit=elcrit
            else
              elacrit=elcrit*(1.0-tca/tlcrit)
            endif
            elacrit=max(elacrit,0.0)
            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
            ep(i,k)=max(ep(i,k),0.0 )
            ep(i,k)=min(ep(i,k),1.0 )
            sigp(i,k)=sigs
          endif
 310    continue
 320  continue
!
!=====================================================================
! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
! --- VIRTUAL TEMPERATURE
!=====================================================================
!
      do 340 k=minorig+1,nl
        do 330 i=1,ncum
        if(k.ge.(icb(i)+1))then
          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
!         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
!         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
        endif
 330    continue
 340  continue
      do 350 i=1,ncum
       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
 350  continue
!
!=====================================================================
!  --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
!  --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
!  --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
!=====================================================================
!
      do 510 i=1,ncum
        cape(i)=0.0
        capem(i)=0.0
        inb(i)=icb(i)+1
        inb1(i)=inb(i)
 510  continue
!
! Originial Code
!
!     do 530 k=minorig+1,nl-1
!       do 520 i=1,ncum
!         if(k.ge.(icb(i)+1))then
!           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
!           byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
!           cape(i)=cape(i)+by
!           if(by.ge.0.0)inb1(i)=k+1
!           if(cape(i).gt.0.0)then
!             inb(i)=k+1
!             capem(i)=cape(i)
!           endif
!         endif
!520    continue
!530  continue
!     do 540 i=1,ncum
!         byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
!         cape(i)=capem(i)+byp
!         defrac=capem(i)-cape(i)
!         defrac=max(defrac,0.001)
!         frac(i)=-cape(i)/defrac
!         frac(i)=min(frac(i),1.0)
!         frac(i)=max(frac(i),0.0)
!540   continue
!
! K Emanuel fix
!
!     call zilch(byp,ncum)
!     do 530 k=minorig+1,nl-1
!       do 520 i=1,ncum
!         if(k.ge.(icb(i)+1))then
!           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
!           cape(i)=cape(i)+by
!           if(by.ge.0.0)inb1(i)=k+1
!           if(cape(i).gt.0.0)then
!             inb(i)=k+1
!             capem(i)=cape(i)
!             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
!           endif
!         endif
!520    continue
!530  continue
!     do 540 i=1,ncum
!         inb(i)=max(inb(i),inb1(i))
!         cape(i)=capem(i)+byp(i)
!         defrac=capem(i)-cape(i)
!         defrac=max(defrac,0.001)
!         frac(i)=-cape(i)/defrac
!         frac(i)=min(frac(i),1.0)
!         frac(i)=max(frac(i),0.0)
!540   continue
!
! J Teixeira fix
!
      call zilch(byp,ncum)
      do 515 i=1,ncum
        lcape(i)=.true.
 515  continue
      do 530 k=minorig+1,nl-1
        do 520 i=1,ncum
          if(cape(i).lt.0.0)lcape(i)=.false.
          if((k.ge.(icb(i)+1)).and.lcape(i))then
            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
            cape(i)=cape(i)+by
            if(by.ge.0.0)inb1(i)=k+1
            if(cape(i).gt.0.0)then
              inb(i)=k+1
              capem(i)=cape(i)
            endif
          endif
 520    continue
 530  continue
      do 540 i=1,ncum
          cape(i)=capem(i)+byp(i)
          defrac=capem(i)-cape(i)
          defrac=max(defrac,0.001)
          frac(i)=-cape(i)/defrac
          frac(i)=min(frac(i),1.0)
          frac(i)=max(frac(i),0.0)
 540  continue
!
!=====================================================================
! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
!=====================================================================
!
! initialization:
      do i=1,ncum*nlp
       hp(i,1)=h(i,1)
      enddo

      do 600 k=minorig+1,nl
        do 590 i=1,ncum
        if((k.ge.icb(i)).and.(k.le.inb(i)))then
          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
        endif
 590    continue
 600  continue
!
        return
        end
1	guez	52
2			SUBROUTINE cv_undilute2(nloc,ncum,nd,icb,nk &
3			,tnk,qnk,gznk,t,q,qs,gz &
4			,p,dph,h,tv,lv &
5			,inb,inb1,tp,tvp,clw,hp,ep,sigp,frac)
6			use cvthermo
7	guez	103	use cv_param
8	guez	52	implicit none
9
10			!---------------------------------------------------------------------
11			! Purpose:
12			! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
13			! &
14			! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
15			! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
16			! &
17			! FIND THE LEVEL OF NEUTRAL BUOYANCY
18			!---------------------------------------------------------------------
19
20
21			! inputs:
22	guez	97	integer, intent(in):: ncum, nd, nloc
23	guez	52	integer icb(nloc), nk(nloc)
24			real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
25			real p(nloc,nd), dph(nloc,nd)
26			real tnk(nloc), qnk(nloc), gznk(nloc)
27			real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
28
29			! outputs:
30			integer inb(nloc), inb1(nloc)
31			real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
32			real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
33			real frac(nloc)
34
35			! local variables:
36			integer i, k
37			real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
38			real by, defrac
39			real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
40			logical lcape(nloc)
41
42			!=====================================================================
43			! --- SOME INITIALIZATIONS
44			!=====================================================================
45
46			do 170 k=1,nl
47			do 160 i=1,ncum
48			ep(i,k)=0.0
49			sigp(i,k)=sigs
50			160 continue
51			170 continue
52
53			!=====================================================================
54			! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
55			!=====================================================================
56			!
57			! --- The procedure is to solve the equation.
58			! cptp+Lqp+phi=cptnk+Lqnk+gznk.
59			!
60			! * Calculate certain parcel quantities, including static energy *
61			!
62			!
63			do 240 i=1,ncum
64			ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i) &
65			+qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
66			240 continue
67			!
68			!
69			! * Find lifted parcel quantities above cloud base *
70			!
71			!
72			do 300 k=minorig+1,nl
73			do 290 i=1,ncum
74			if(k.ge.(icb(i)+1))then
75			tg=t(i,k)
76			qg=qs(i,k)
77			alv=lv0-clmcpv*(t(i,k)-t0)
78			!
79			! First iteration.
80			!
81			s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
82			s=1./s
83			ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
84			tg=tg+s*(ah0(i)-ahg)
85			tg=max(tg,35.0)
86			tc=tg-t0
87			denom=243.5+tc
88			if(tc.ge.0.0)then
89			es=6.112exp(17.67tc/denom)
90			else
91			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
92			endif
93			qg=epses/(p(i,k)-es(1.-eps))
94			!
95			! Second iteration.
96			!
97			s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
98			s=1./s
99			ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
100			tg=tg+s*(ah0(i)-ahg)
101			tg=max(tg,35.0)
102			tc=tg-t0
103			denom=243.5+tc
104			if(tc.ge.0.0)then
105			es=6.112exp(17.67tc/denom)
106			else
107			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
108			endif
109			qg=epses/(p(i,k)-es(1.-eps))
110			!
111			alv=lv0-clmcpv*(t(i,k)-t0)
112			! print*,'cpd dans convect2 ',cpd
113			! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
114			! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
115			tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
116			! if (.not.cpd.gt.1000.) then
117			! print*,'CPD=',cpd
118			! stop
119			! endif
120			clw(i,k)=qnk(i)-qg
121			clw(i,k)=max(0.0,clw(i,k))
122			rg=qg/(1.-qnk(i))
123			tvp(i,k)=tp(i,k)(1.+rgepsi)
124			endif
125			290 continue
126			300 continue
127			!
128			!=====================================================================
129			! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
130			! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
131			! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
132			!=====================================================================
133			!
134			do 320 k=minorig+1,nl
135			do 310 i=1,ncum
136			if(k.ge.(nk(i)+1))then
137			tca=tp(i,k)-t0
138			if(tca.ge.0.0)then
139			elacrit=elcrit
140			else
141			elacrit=elcrit*(1.0-tca/tlcrit)
142			endif
143			elacrit=max(elacrit,0.0)
144			ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
145			ep(i,k)=max(ep(i,k),0.0 )
146			ep(i,k)=min(ep(i,k),1.0 )
147			sigp(i,k)=sigs
148			endif
149			310 continue
150			320 continue
151			!
152			!=====================================================================
153			! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
154			! --- VIRTUAL TEMPERATURE
155			!=====================================================================
156			!
157			do 340 k=minorig+1,nl
158			do 330 i=1,ncum
159			if(k.ge.(icb(i)+1))then
160			tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
161			! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
162			! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
163			endif
164			330 continue
165			340 continue
166			do 350 i=1,ncum
167			tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
168			350 continue
169			!
170			!=====================================================================
171			! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
172			! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
173			! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
174			!=====================================================================
175			!
176			do 510 i=1,ncum
177			cape(i)=0.0
178			capem(i)=0.0
179			inb(i)=icb(i)+1
180			inb1(i)=inb(i)
181			510 continue
182			!
183			! Originial Code
184			!
185			! do 530 k=minorig+1,nl-1
186			! do 520 i=1,ncum
187			! if(k.ge.(icb(i)+1))then
188			! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
189			! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
190			! cape(i)=cape(i)+by
191			! if(by.ge.0.0)inb1(i)=k+1
192			! if(cape(i).gt.0.0)then
193			! inb(i)=k+1
194			! capem(i)=cape(i)
195			! endif
196			! endif
197			!520 continue
198			!530 continue
199			! do 540 i=1,ncum
200			! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
201			! cape(i)=capem(i)+byp
202			! defrac=capem(i)-cape(i)
203			! defrac=max(defrac,0.001)
204			! frac(i)=-cape(i)/defrac
205			! frac(i)=min(frac(i),1.0)
206			! frac(i)=max(frac(i),0.0)
207			!540 continue
208			!
209			! K Emanuel fix
210			!
211			! call zilch(byp,ncum)
212			! do 530 k=minorig+1,nl-1
213			! do 520 i=1,ncum
214			! if(k.ge.(icb(i)+1))then
215			! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
216			! cape(i)=cape(i)+by
217			! if(by.ge.0.0)inb1(i)=k+1
218			! if(cape(i).gt.0.0)then
219			! inb(i)=k+1
220			! capem(i)=cape(i)
221			! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
222			! endif
223			! endif
224			!520 continue
225			!530 continue
226			! do 540 i=1,ncum
227			! inb(i)=max(inb(i),inb1(i))
228			! cape(i)=capem(i)+byp(i)
229			! defrac=capem(i)-cape(i)
230			! defrac=max(defrac,0.001)
231			! frac(i)=-cape(i)/defrac
232			! frac(i)=min(frac(i),1.0)
233			! frac(i)=max(frac(i),0.0)
234			!540 continue
235			!
236			! J Teixeira fix
237			!
238			call zilch(byp,ncum)
239			do 515 i=1,ncum
240			lcape(i)=.true.
241			515 continue
242			do 530 k=minorig+1,nl-1
243			do 520 i=1,ncum
244			if(cape(i).lt.0.0)lcape(i)=.false.
245			if((k.ge.(icb(i)+1)).and.lcape(i))then
246			by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
247			byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
248			cape(i)=cape(i)+by
249			if(by.ge.0.0)inb1(i)=k+1
250			if(cape(i).gt.0.0)then
251			inb(i)=k+1
252			capem(i)=cape(i)
253			endif
254			endif
255			520 continue
256			530 continue
257			do 540 i=1,ncum
258			cape(i)=capem(i)+byp(i)
259			defrac=capem(i)-cape(i)
260			defrac=max(defrac,0.001)
261			frac(i)=-cape(i)/defrac
262			frac(i)=min(frac(i),1.0)
263			frac(i)=max(frac(i),0.0)
264			540 continue
265			!
266			!=====================================================================
267			! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
268			!=====================================================================
269			!
270			! initialization:
271			do i=1,ncum*nlp
272			hp(i,1)=h(i,1)
273			enddo
274
275			do 600 k=minorig+1,nl
276			do 590 i=1,ncum
277			if((k.ge.icb(i)).and.(k.le.inb(i)))then
278			hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
279			endif
280			590 continue
281			600 continue
282			!
283			return
284			end