/[lmdze]/trunk/Sources/phylmd/CV3_routines/cv3_undilute2.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV3_routines/cv3_undilute2.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 182 by guez, Fri Mar 11 18:47:26 2016 UTC revision 183 by guez, Wed Mar 16 14:42:58 2016 UTC
# Line 1  Line 1 
1    module cv3_undilute2_m
2    
3        SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk &    implicit none
                              ,tnk,qnk,gznk,t,qs,gz &  
                              ,p,h,tv,lv,pbase,buoybase,plcl &  
                              ,inb,tp,tvp,clw,hp,ep,sigp,buoy)  
       use conema3_m  
             use cv3_param_m  
             use cvthermo  
       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  
 !  
 ! Main differences convect3/convect4:  
 !     - icbs (input) is the first level above LCL (may differ from icb)  
 !     - many minor differences in the iterations  
 !     - condensed water not removed from tvp in convect3  
 !   - vertical profile of buoyancy computed here (use of buoybase)  
 !   - the determination of inb is different  
 !   - no inb1, only inb in output  
 !---------------------------------------------------------------------  
   
   
 ! inputs:  
       integer, intent(in):: ncum, nd, nloc  
       integer icb(nloc), icbs(nloc), nk(nloc)  
       real t(nloc,nd), qs(nloc,nd), gz(nloc,nd)  
       real p(nloc,nd)  
       real tnk(nloc), qnk(nloc), gznk(nloc)  
       real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)  
       real pbase(nloc), buoybase(nloc), plcl(nloc)  
   
 ! outputs:  
       integer inb(nloc)  
       real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)  
       real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)  
       real buoy(nloc,nd)  
   
 ! local variables:  
       integer i, k  
       real tg,qg,ahg,alv,s,tc,es,denom  
       real pden  
       real ah0(nloc)  
   
 !=====================================================================  
 ! --- SOME INITIALIZATIONS  
 !=====================================================================  
   
       do 170 k=1,nl  
       do 160 i=1,ncum  
        ep(i,k)=0.0  
        sigp(i,k)=spfac  
  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)-273.15))+gznk(i)  
  240  continue  
 !  
 !  
 !    ***  Find lifted parcel quantities above cloud base    ***  
 !  
 !  
       do 300 k=minorig+1,nl  
         do 290 i=1,ncum  
 ! ori        if(k.ge.(icb(i)+1))then  
           if(k.ge.(icbs(i)+1))then ! convect3  
             tg=t(i,k)  
             qg=qs(i,k)  
 !debug         alv=lv0-clmcpv*(t(i,k)-t0)  
             alv=lv0-clmcpv*(t(i,k)-273.15)  
 !  
 ! First iteration.  
 !  
 ! ori           s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))  
            s=cpd*(1.-qnk(i))+cl*qnk(i) &  
             +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3  
              s=1./s  
 ! ori           ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)  
            ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3  
              tg=tg+s*(ah0(i)-ahg)  
 ! ori           tg=max(tg,35.0)  
 !debug          tc=tg-t0  
              tc=tg-273.15  
              denom=243.5+tc  
            denom=MAX(denom,1.0) ! convect3  
 ! ori           if(tc.ge.0.0)then  
             es=6.112*exp(17.67*tc/denom)  
 ! ori           else  
 ! ori          es=exp(23.33086-6111.72784/tg+0.15215*log(tg))  
 ! ori           endif  
             qg=eps*es/(p(i,k)-es*(1.-eps))  
 !  
 ! Second iteration.  
 !  
 ! ori           s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))  
 ! ori           s=1./s  
 ! ori           ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)  
            ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3  
              tg=tg+s*(ah0(i)-ahg)  
 ! ori           tg=max(tg,35.0)  
 !debug          tc=tg-t0  
              tc=tg-273.15  
              denom=243.5+tc  
            denom=MAX(denom,1.0) ! convect3  
 ! ori           if(tc.ge.0.0)then  
             es=6.112*exp(17.67*tc/denom)  
 ! ori           else  
 ! ori          es=exp(23.33086-6111.72784/tg+0.15215*log(tg))  
 ! ori           endif  
             qg=eps*es/(p(i,k)-es*(1.-eps))  
 !  
 !debug          alv=lv0-clmcpv*(t(i,k)-t0)  
              alv=lv0-clmcpv*(t(i,k)-273.15)  
 !      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  
   
 ! ori c approximation here:  
 ! ori        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd  
   
 ! convect3: no approximation:  
            tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))  
   
                clw(i,k)=qnk(i)-qg  
                clw(i,k)=max(0.0,clw(i,k))  
 ! convect3: (qg utilise au lieu du vrai mixing ratio rg):  
                tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing  
             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)  
 !=====================================================================  
 !  
 ! ori      do 320 k=minorig+1,nl  
       do 320 k=1,nl ! convect3  
         do 310 i=1,ncum  
            pden=ptcrit-pbcrit  
            ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax  
            ep(i,k)=amax1(ep(i,k),0.0)  
            ep(i,k)=amin1(ep(i,k),epmax)  
            sigp(i,k)=spfac  
 ! ori          if(k.ge.(nk(i)+1))then  
 ! ori            tca=tp(i,k)-t0  
 ! ori            if(tca.ge.0.0)then  
 ! ori              elacrit=elcrit  
 ! ori            else  
 ! ori              elacrit=elcrit*(1.0-tca/tlcrit)  
 ! ori            endif  
 ! ori            elacrit=max(elacrit,0.0)  
 ! ori            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)  
 ! ori            ep(i,k)=max(ep(i,k),0.0 )  
 ! ori            ep(i,k)=min(ep(i,k),1.0 )  
 ! ori            sigp(i,k)=sigs  
 ! ori          endif  
  310    continue  
  320  continue  
 !  
 !=====================================================================  
 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL  
 ! --- VIRTUAL TEMPERATURE  
 !=====================================================================  
 !  
 ! dans convect3, tvp est calcule en une seule fois, et sans retirer  
 ! l'eau condensee (~> reversible CAPE)  
 !  
 ! ori      do 340 k=minorig+1,nl  
 ! ori        do 330 i=1,ncum  
 ! ori        if(k.ge.(icb(i)+1))then  
 ! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))  
 ! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'  
 ! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)  
 ! ori        endif  
 ! ori 330    continue  
 ! ori 340  continue  
   
 ! ori      do 350 i=1,ncum  
 ! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd  
 ! ori 350  continue  
   
       do 350 i=1,ncum       ! convect3  
        tp(i,nlp)=tp(i,nl)   ! convect3  
  350  continue              ! convect3  
 !  
 !=====================================================================  
 !  --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):  
 !=====================================================================  
   
 !-- this is for convect3 only:  
   
 ! first estimate of buoyancy:  
   
       do 500 i=1,ncum  
        do 501 k=1,nl  
         buoy(i,k)=tvp(i,k)-tv(i,k)  
  501   continue  
  500  continue  
   
 ! set buoyancy=buoybase for all levels below base  
 ! for safety, set buoy(icb)=buoybase  
   
       do 505 i=1,ncum  
        do 506 k=1,nl  
         if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then  
          buoy(i,k)=buoybase(i)  
         endif  
  506   continue  
        buoy(icb(i),k)=buoybase(i)  
  505  continue  
   
 !-- end convect3  
   
 !=====================================================================  
 !  --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S  
 !  --- LEVEL OF NEUTRAL BUOYANCY  
 !=====================================================================  
 !  
 !-- this is for convect3 only:  
   
       do 510 i=1,ncum  
        inb(i)=nl-1  
  510  continue  
   
       do 530 i=1,ncum  
        do 535 k=1,nl-1  
         if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then  
          inb(i)=MIN(inb(i),k)  
         endif  
  535   continue  
  530  continue  
   
 !-- end convect3  
   
 ! ori      do 510 i=1,ncum  
 ! ori        cape(i)=0.0  
 ! ori        capem(i)=0.0  
 ! ori        inb(i)=icb(i)+1  
 ! ori        inb1(i)=inb(i)  
 ! ori 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  
 !  
 ! ori      call zilch(byp,ncum)  
 ! ori      do 515 i=1,ncum  
 ! ori        lcape(i)=.true.  
 ! ori 515  continue  
 ! ori      do 530 k=minorig+1,nl-1  
 ! ori        do 520 i=1,ncum  
 ! ori          if(cape(i).lt.0.0)lcape(i)=.false.  
 ! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then  
 ! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)  
 ! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)  
 ! ori            cape(i)=cape(i)+by  
 ! ori            if(by.ge.0.0)inb1(i)=k+1  
 ! ori            if(cape(i).gt.0.0)then  
 ! ori              inb(i)=k+1  
 ! ori              capem(i)=cape(i)  
 ! ori            endif  
 ! ori          endif  
 ! ori 520    continue  
 ! ori 530  continue  
 ! ori      do 540 i=1,ncum  
 ! ori          cape(i)=capem(i)+byp(i)  
 ! ori          defrac=capem(i)-cape(i)  
 ! ori          defrac=max(defrac,0.001)  
 ! ori          frac(i)=-cape(i)/defrac  
 ! ori          frac(i)=min(frac(i),1.0)  
 ! ori          frac(i)=max(frac(i),0.0)  
 ! ori 540  continue  
 !  
 !=====================================================================  
 ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL  
 !=====================================================================  
 !  
 !ym      do i=1,ncum*nlp  
 !ym       hp(i,1)=h(i,1)  
 !ym      enddo  
   
       do k=1,nlp  
         do i=1,ncum  
         hp(i,k)=h(i,k)  
       enddo  
       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  
4    
5          return  contains
6          end  
7      SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
8           qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
9           ep, sigp, buoy)
10    
11        ! Purpose: find the rest of the lifted parcel temperatures;
12        ! compute the precipitation efficiencies and the fraction of
13        ! precipitation falling outside of cloud; find the level of
14        ! neutral buoyancy.
15    
16        ! Vertical profile of buoyancy computed here (use of buoybase).
17    
18        use conema3_m, only: epmax
19        use cv3_param_m, only: dtovsh, minorig, nl, nlp, pbcrit, ptcrit, spfac
20        use cvthermo, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
21    
22        ! inputs:
23        integer, intent(in):: nloc, ncum, nd
24        integer icb(nloc), icbs(nloc), nk(nloc)
25        ! icbs (input) is the first level above LCL (may differ from icb)
26        real tnk(nloc), qnk(nloc), gznk(nloc)
27        real t(nloc, nd), qs(nloc, nd), gz(nloc, nd)
28        real p(nloc, nd), h(nloc, nd)
29        real tv(nloc, nd), lv(nloc, nd)
30        real pbase(nloc), buoybase(nloc), plcl(nloc)
31    
32        ! outputs:
33        integer inb(nloc)
34        real tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
35        ! condensed water not removed from tvp
36        real hp(nloc, nd), ep(nloc, nd), sigp(nloc, nd)
37        real buoy(nloc, nd)
38    
39        ! Local:
40        integer i, k
41        real tg, qg, ahg, alv, s, tc, es, denom
42        real pden
43        real ah0(nloc)
44    
45        !---------------------------------------------------------------------
46    
47        ! SOME INITIALIZATIONS
48    
49        do k = 1, nl
50           do i = 1, ncum
51              ep(i, k) = 0.0
52              sigp(i, k) = spfac
53           end do
54        end do
55    
56        ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
57    
58        ! The procedure is to solve the equation.
59        ! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.
60    
61        ! Calculate certain parcel quantities, including static energy
62    
63        do i = 1, ncum
64           ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) &
65                + qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i)
66        end do
67    
68        ! Find lifted parcel quantities above cloud base
69    
70        do k = minorig + 1, nl
71           do i = 1, ncum
72              if (k >= (icbs(i) + 1)) then
73                 tg = t(i, k)
74                 qg = qs(i, k)
75                 alv = lv0 - clmcpv * (t(i, k) - 273.15)
76    
77                 ! First iteration.
78    
79                 s = cpd * (1. - qnk(i)) + cl * qnk(i) &
80                      + alv * alv * qg / (rrv * t(i, k) * t(i, k))
81                 s = 1. / s
82    
83                 ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
84                 tg = tg + s * (ah0(i) - ahg)
85    
86                 tc = tg - 273.15
87                 denom = 243.5 + tc
88                 denom = MAX(denom, 1.0)
89    
90                 es = 6.112 * exp(17.67 * tc / denom)
91    
92                 qg = eps * es / (p(i, k) - es * (1. - eps))
93    
94                 ! Second iteration.
95    
96                 ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
97                 tg = tg + s * (ah0(i) - ahg)
98    
99                 tc = tg - 273.15
100                 denom = 243.5 + tc
101                 denom = MAX(denom, 1.0)
102    
103                 es = 6.112 * exp(17.67 * tc / denom)
104    
105                 qg = eps * es / (p(i, k) - es * (1. - eps))
106    
107                 alv = lv0 - clmcpv * (t(i, k) - 273.15)
108    
109                 ! no approximation:
110                 tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) &
111                      / (cpd + (cl - cpd) * qnk(i))
112    
113                 clw(i, k) = qnk(i) - qg
114                 clw(i, k) = max(0.0, clw(i, k))
115                 ! qg utilise au lieu du vrai mixing ratio rg:
116                 tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing
117              endif
118           end do
119        end do
120    
121        ! SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
122        ! PRECIPITATION FALLING OUTSIDE OF CLOUD
123        ! THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
124        do k = 1, nl
125           do i = 1, ncum
126              pden = ptcrit - pbcrit
127              ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax
128              ep(i, k) = amax1(ep(i, k), 0.0)
129              ep(i, k) = amin1(ep(i, k), epmax)
130              sigp(i, k) = spfac
131           end do
132        end do
133    
134        ! CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
135        ! VIRTUAL TEMPERATURE
136    
137        ! tvp est calcule en une seule fois, et sans retirer
138        ! l'eau condensee (~> reversible CAPE)
139        do i = 1, ncum
140           tp(i, nlp) = tp(i, nl)
141        end do
142    
143        ! EFFECTIVE VERTICAL PROFILE OF BUOYANCY:
144    
145        ! first estimate of buoyancy:
146        do i = 1, ncum
147           do k = 1, nl
148              buoy(i, k) = tvp(i, k) - tv(i, k)
149           end do
150        end do
151    
152        ! set buoyancy = buoybase for all levels below base
153        ! for safety, set buoy(icb) = buoybase
154        do i = 1, ncum
155           do k = 1, nl
156              if ((k >= icb(i)) .and. (k <= nl) .and. (p(i, k) >= pbase(i))) then
157                 buoy(i, k) = buoybase(i)
158              endif
159           end do
160           buoy(icb(i), k) = buoybase(i)
161        end do
162    
163        ! FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
164        ! LEVEL OF NEUTRAL BUOYANCY
165    
166        do i = 1, ncum
167           inb(i) = nl - 1
168        end do
169    
170        do i = 1, ncum
171           do k = 1, nl - 1
172              if ((k >= icb(i)) .and. (buoy(i, k) < dtovsh)) then
173                 inb(i) = MIN(inb(i), k)
174              endif
175           end do
176        end do
177    
178        ! CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
179    
180        do k = 1, nlp
181           do i = 1, ncum
182              hp(i, k) = h(i, k)
183           enddo
184        enddo
185    
186        do k = minorig + 1, nl
187           do i = 1, ncum
188              if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, nk(i)) &
189                   + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
190           end do
191        end do
192    
193      end SUBROUTINE cv3_undilute2
194    
195    end module cv3_undilute2_m
Legend:
Removed from v.182  
changed lines
  Added in v.183
  ViewVC Help
Powered by ViewVC 1.1.21