/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_undilute2.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_undilute2.f

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

trunk/Sources/phylmd/CV3_routines/cv3_undilute2.f revision 178 by guez, Fri Mar 11 18:47:26 2016 UTC trunk/Sources/phylmd/CV30_routines/cv30_undilute2.f revision 198 by guez, Tue May 31 16:17:35 2016 UTC
# Line 1  Line 1 
1    module cv30_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 cv30_undilute2(icb, icbs, tnk, qnk, gznk, t, qs, gz, p, h, tv, &
8           lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, ep, buoy)
9    
10        ! Undilute (adiabatic) updraft, second part. Purpose: find the
11        ! rest of the lifted parcel temperatures; compute the
12        ! precipitation efficiencies and the fraction of precipitation
13        ! falling outside of cloud; find the level of neutral buoyancy.
14    
15        ! Vertical profile of buoyancy computed here (use of buoybase).
16    
17        use conema3_m, only: epmax
18        use cv30_param_m, only: minorig, nl
19        use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, rrv
20        USE dimphy, ONLY: klon, klev
21        use SUPHEC_M, only: rlvtt
22    
23        integer, intent(in):: icb(:), icbs(:) ! (ncum)
24        ! icbs is the first level above LCL (may differ from icb)
25    
26        real, intent(in):: tnk(:), qnk(:), gznk(:) ! (klon)
27        real, intent(in):: t(klon, klev), qs(klon, klev), gz(klon, klev)
28        real, intent(in):: p(klon, klev), h(klon, klev)
29        real, intent(in):: tv(klon, klev), lv(klon, klev)
30        real, intent(in):: pbase(:), buoybase(:), plcl(:) ! (ncum)
31    
32        ! outputs:
33        integer, intent(out):: inb(:) ! (ncum)
34        ! first model level above the level of neutral buoyancy of the
35        ! parcel (1 <= inb <= nl - 1)
36    
37        real tp(klon, klev), tvp(klon, klev), clw(klon, klev)
38        ! condensed water not removed from tvp
39        real hp(klon, klev), ep(klon, klev)
40        real buoy(klon, klev)
41    
42        ! Local:
43    
44        integer ncum
45    
46        real, parameter:: pbcrit = 150.
47        ! critical cloud depth (mbar) beneath which the precipitation
48        ! efficiency is assumed to be zero
49    
50        real, parameter:: ptcrit = 500.
51        ! cloud depth (mbar) above which the precipitation efficiency is
52        ! assumed to be unity
53    
54        real, parameter:: dtovsh = - 0.2 ! dT for overshoot
55    
56        integer i, k
57        real tg, qg, ahg, alv, s, tc, es, denom
58        real pden
59        real ah0(klon)
60    
61        !---------------------------------------------------------------------
62    
63        ncum = size(icb)
64    
65        ! SOME INITIALIZATIONS
66    
67        do k = 1, nl
68           do i = 1, ncum
69              ep(i, k) = 0.
70           end do
71        end do
72    
73        ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
74    
75        ! The procedure is to solve the equation.
76        ! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.
77    
78        ! Calculate certain parcel quantities, including static energy
79    
80        do i = 1, ncum
81           ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) &
82                + qnk(i) * (rlvtt - clmcpv * (tnk(i) - 273.15)) + gznk(i)
83        end do
84    
85        ! Find lifted parcel quantities above cloud base
86    
87        do k = minorig + 1, nl
88           do i = 1, ncum
89              if (k >= (icbs(i) + 1)) then
90                 tg = t(i, k)
91                 qg = qs(i, k)
92                 alv = rlvtt - clmcpv * (t(i, k) - 273.15)
93    
94                 ! First iteration.
95    
96                 s = cpd * (1. - qnk(i)) + cl * qnk(i) &
97                      + alv * alv * qg / (rrv * t(i, k) * t(i, k))
98                 s = 1. / s
99    
100                 ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
101                 tg = tg + s * (ah0(i) - ahg)
102    
103                 tc = tg - 273.15
104                 denom = 243.5 + tc
105                 denom = MAX(denom, 1.)
106    
107                 es = 6.112 * exp(17.67 * tc / denom)
108    
109                 qg = eps * es / (p(i, k) - es * (1. - eps))
110    
111                 ! Second iteration.
112    
113                 ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
114                 tg = tg + s * (ah0(i) - ahg)
115    
116                 tc = tg - 273.15
117                 denom = 243.5 + tc
118                 denom = MAX(denom, 1.)
119    
120                 es = 6.112 * exp(17.67 * tc / denom)
121    
122                 qg = eps * es / (p(i, k) - es * (1. - eps))
123    
124                 alv = rlvtt - clmcpv * (t(i, k) - 273.15)
125    
126                 ! no approximation:
127                 tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) &
128                      / (cpd + (cl - cpd) * qnk(i))
129    
130                 clw(i, k) = qnk(i) - qg
131                 clw(i, k) = max(0., clw(i, k))
132                 ! qg utilise au lieu du vrai mixing ratio rg:
133                 tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing
134              endif
135           end do
136        end do
137    
138        ! SET THE PRECIPITATION EFFICIENCIES
139        ! It MAY BE a FUNCTION OF TP(I), P(I) AND CLW(I)
140        do k = 1, nl
141           do i = 1, ncum
142              pden = ptcrit - pbcrit
143              ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax
144              ep(i, k) = max(ep(i, k), 0.)
145              ep(i, k) = min(ep(i, k), epmax)
146           end do
147        end do
148    
149        ! CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
150        ! VIRTUAL TEMPERATURE
151    
152        ! tvp est calcule en une seule fois, et sans retirer
153        ! l'eau condensee (~> reversible CAPE)
154        do i = 1, ncum
155           tp(i, nl + 1) = tp(i, nl)
156        end do
157    
158        ! EFFECTIVE VERTICAL PROFILE OF BUOYANCY:
159    
160        ! first estimate of buoyancy:
161        do i = 1, ncum
162           do k = 1, nl
163              buoy(i, k) = tvp(i, k) - tv(i, k)
164           end do
165        end do
166    
167        ! set buoyancy = buoybase for all levels below base
168        ! for safety, set buoy(icb) = buoybase
169        do i = 1, ncum
170           do k = 1, nl
171              if ((k >= icb(i)) .and. (k <= nl) .and. (p(i, k) >= pbase(i))) then
172                 buoy(i, k) = buoybase(i)
173              endif
174           end do
175           buoy(icb(i), k) = buoybase(i)
176        end do
177    
178        ! Compute inb:
179    
180        inb = nl - 1
181    
182        do i = 1, ncum
183           do k = 1, nl - 1
184              if ((k >= icb(i)) .and. (buoy(i, k) < dtovsh)) then
185                 inb(i) = MIN(inb(i), k)
186              endif
187           end do
188        end do
189    
190        ! CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
191    
192        do k = 1, nl + 1
193           do i = 1, ncum
194              hp(i, k) = h(i, k)
195           enddo
196        enddo
197    
198        do k = minorig + 1, nl
199           do i = 1, ncum
200              if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, minorig) &
201                   + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
202           end do
203        end do
204    
205      end SUBROUTINE cv30_undilute2
206    
207    end module cv30_undilute2_m
Legend:
Removed from v.178  
changed lines
  Added in v.198
  ViewVC Help
Powered by ViewVC 1.1.21