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

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

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

revision 189 by guez, Tue Mar 29 15:20:23 2016 UTC revision 196 by guez, Mon May 23 13:50:39 2016 UTC
# Line 1  Line 1 
1    module cv30_undilute1_m
2    
3        SUBROUTINE cv30_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb &    implicit none
4                               ,tp,tvp,clw,icbs)  
5              use cv30_param_m  contains
6              use cvthermo  
7        implicit none    SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, nk1, icb1, tp1, &
8           tvp1, clw1, icbs1)
9    
10      ! UNDILUTE (ADIABATIC) UPDRAFT / 1st part      ! UNDILUTE (ADIABATIC) UPDRAFT / 1st part
11      ! (up through ICB for convect4, up through ICB + 1 for convect3)      ! (up through ICB1 + 1)
12      ! Calculates the lifted parcel virtual temperature at nk, the      ! Calculates the lifted parcel virtual temperature at nk1, the
13      ! actual temperature, and the adiabatic liquid water content.      ! actual temperature, and the adiabatic liquid water content.
14    
15  !----------------------------------------------------------------      ! Equivalent de TLIFT entre NK1 et ICB1+1 inclus
16  ! Equivalent de TLIFT entre NK et ICB+1 inclus  
17  !      ! Differences with convect4:
18  ! Differences with convect4:      ! - icbs1 is the first level above LCL (may differ from icb1)
19  !        - specify plcl in input      ! - in the iterations, used x(icbs1) instead x(icb1)
20  !       - icbs is the first level above LCL (may differ from icb)      ! - tvp1 is computed in only one time
21  !       - in the iterations, used x(icbs) instead x(icb)      ! - icbs1: first level above Plcl1 (IMIN de TLIFT) in output
22  !       - many minor differences in the iterations      ! - if icbs1=icb1, compute also tp1(icb1+1), tvp1(icb1+1) & clw1(icb1+1)
23  !        - tvp is computed in only one time  
24  !        - icbs: first level above Plcl (IMIN de TLIFT) in output      use cv30_param_m, only: minorig, nl
25  !       - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)      use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
26  !----------------------------------------------------------------      USE dimphy, ONLY: klev, klon
27    
28        ! inputs:
29  ! inputs:      integer, intent(in):: nk1(klon), icb1(klon)
30        integer, intent(in):: len, nd      real, intent(in):: t1(klon, klev)
31        integer nk(len), icb(len)      real, intent(in):: q1(klon, klev), qs1(klon, klev), gz1(klon, klev)
32        real, intent(in):: t(len,nd)      real, intent(in):: p1(klon, klev)
33        real, intent(in):: q(len,nd), qs(len,nd), gz(len,nd)      real, intent(in):: plcl1(klon)
34        real, intent(in):: p(len,nd)  
35        real plcl(len) ! convect3      ! outputs:
36        real tp1(klon, klev), tvp1(klon, klev), clw1(klon, klev)
37  ! outputs:  
38        real tp(len,nd), tvp(len,nd), clw(len,nd)      ! local variables:
39        integer i, k
40  ! local variables:      integer icbs1(klon), icbsmax2
41        integer i, k      real tg, qg, alv, s, ahg, tc, denom, es
42        integer icb1(len), icbs(len), icbsmax2 ! convect3      real ah0(klon), cpp(klon)
43        real tg, qg, alv, s, ahg, tc, denom, es      real tnk(klon), qnk(klon), gznk(klon), ticb(klon), gzicb(klon)
44        real ah0(len), cpp(len)      real qsicb(klon)
45        real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)      real cpinv(klon)
46        real qsicb(len) ! convect3  
47        real cpinv(len) ! convect3      !-------------------------------------------------------------------
48    
49  !-------------------------------------------------------------------      !  Calculates the lifted parcel virtual temperature at nk1,
50  ! --- Calculates the lifted parcel virtual temperature at nk,      !  the actual temperature, and the adiabatic
51  ! --- the actual temperature, and the adiabatic      !  liquid water content. The procedure is to solve the equation.
52  ! --- liquid water content. The procedure is to solve the equation.      ! cp*tp1+L*qp+phi=cp*tnk+L*qnk+gznk.
53  !     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.  
54  !-------------------------------------------------------------------      do i=1, klon
55           tnk(i)=t1(i, nk1(i))
56        do 320 i=1,len         qnk(i)=q1(i, nk1(i))
57          tnk(i)=t(i,nk(i))         gznk(i)=gz1(i, nk1(i))
58          qnk(i)=q(i,nk(i))      end do
59          gznk(i)=gz(i,nk(i))  
60  ! ori        ticb(i)=t(i,icb(i))      ! *** Calculate certain parcel quantities, including static energy ***
61  ! ori        gzicb(i)=gz(i,icb(i))  
62   320  continue      do i=1, klon
63  !         ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &
64  !   ***  Calculate certain parcel quantities, including static energy   ***              +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
65  !         cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
66        do 330 i=1,len         cpinv(i)=1./cpp(i)
67          ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &      end do
68                 +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)  
69          cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv      ! *** Calculate lifted parcel quantities below cloud base ***
70          cpinv(i)=1./cpp(i)  
71   330  continue      do i=1, klon
72  !         ! if icb1 is below LCL, start loop at ICB1+1:
73  !   ***   Calculate lifted parcel quantities below cloud base   ***         ! (icbs1 est le premier niveau au-dessus du LCL)
74  !         icbs1(i)=MIN(max(icb1(i), 2), nl)
75          do i=1,len                      !convect3         if (plcl1(i) < p1(i, icbs1(i))) then
76           icb1(i)=MAX(icb(i),2)          !convect3            icbs1(i)=MIN(icbs1(i)+1, nl)
77           icb1(i)=MIN(icb(i),nl)         !convect3         endif
78  ! if icb is below LCL, start loop at ICB+1:      enddo
79  ! (icbs est le premier niveau au-dessus du LCL)  
80           icbs(i)=icb1(i)                !convect3      do i=1, klon
81           if (plcl(i).lt.p(i,icb1(i))) then         ticb(i)=t1(i, icbs1(i))
82               icbs(i)=MIN(icbs(i)+1,nl)  !convect3         gzicb(i)=gz1(i, icbs1(i))
83           endif         qsicb(i)=qs1(i, icbs1(i))
84          enddo                           !convect3      enddo
85    
86          do i=1,len                      !convect3      ! Re-compute icbsmax (icbsmax2):
87           ticb(i)=t(i,icbs(i))           !convect3      icbsmax2=2
88           gzicb(i)=gz(i,icbs(i))         !convect3      do i=1, klon
89           qsicb(i)=qs(i,icbs(i))         !convect3         icbsmax2=max(icbsmax2, icbs1(i))
90          enddo                           !convect3      end do
91    
92  !      ! initialization outputs:
93  ! Re-compute icbsmax (icbsmax2):        !convect3  
94  !                                       !convect3      do k=1, icbsmax2
95        icbsmax2=2                        !convect3         do i=1, klon
96        do 310 i=1,len                    !convect3            tp1(i, k) = 0.0
97          icbsmax2=max(icbsmax2,icbs(i))  !convect3            tvp1(i, k) = 0.0
98   310  continue                          !convect3            clw1(i, k) = 0.0
99           enddo
100  ! initialization outputs:      enddo
101    
102        do k=1,icbsmax2     ! convect3      ! tp1 and tvp1 below cloud base:
103         do i=1,len         ! convect3  
104          tp(i,k)  = 0.0    ! convect3      do k=minorig, icbsmax2-1
105          tvp(i,k) = 0.0    ! convect3         do i=1, klon
106          clw(i,k) = 0.0    ! convect3            tp1(i, k)=tnk(i)-(gz1(i, k)-gznk(i))*cpinv(i)
107         enddo              ! convect3            tvp1(i, k)=tp1(i, k)*(1.+qnk(i)/eps-qnk(i))
108        enddo               ! convect3         end do
109        end do
110  ! tp and tvp below cloud base:  
111        ! *** Find lifted parcel quantities above cloud base ***
112          do 350 k=minorig,icbsmax2-1  
113            do 340 i=1,len      do i=1, klon
114             tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))*cpinv(i)         tg=ticb(i)
115             tvp(i,k)=tp(i,k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)         qg=qsicb(i)
116    340     continue         !debug alv=lv0-clmcpv*(ticb(i)-t0)
117    350   continue         alv=lv0-clmcpv*(ticb(i)-273.15)
118  !  
119  !    ***  Find lifted parcel quantities above cloud base    ***         ! First iteration.
120  !  
121          do 360 i=1,len         s=cpd*(1.-qnk(i))+cl*qnk(i) &
122           tg=ticb(i)              +alv*alv*qg/(rrv*ticb(i)*ticb(i))
123  ! ori         qg=qs(i,icb(i))         s=1./s
124           qg=qsicb(i) ! convect3  
125  !debug         alv=lv0-clmcpv*(ticb(i)-t0)         ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
126           alv=lv0-clmcpv*(ticb(i)-273.15)         tg=tg+s*(ah0(i)-ahg)
127  !  
128  ! First iteration.         !debug tc=tg-t0
129  !         tc=tg-273.15
130  ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))         denom=243.5+tc
131            s=cpd*(1.-qnk(i))+cl*qnk(i) &         denom=MAX(denom, 1.0)
132              +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3  
133            s=1./s         es=6.112*exp(17.67*tc/denom)
134  ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)         qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))
135            ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3  
136            tg=tg+s*(ah0(i)-ahg)         ! Second iteration.
137  ! ori          tg=max(tg,35.0)  
138  !debug          tc=tg-t0         ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
139            tc=tg-273.15         tg=tg+s*(ah0(i)-ahg)
140            denom=243.5+tc  
141            denom=MAX(denom,1.0) ! convect3         !debug tc=tg-t0
142  ! ori          if(tc.ge.0.0)then         tc=tg-273.15
143             es=6.112*exp(17.67*tc/denom)         denom=243.5+tc
144  ! ori          else         denom=MAX(denom, 1.0)
145  ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))  
146  ! ori          endif         es=6.112*exp(17.67*tc/denom)
147  ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))  
148            qg=eps*es/(p(i,icbs(i))-es*(1.-eps))         qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))
149  !  
150  ! Second iteration.         alv=lv0-clmcpv*(ticb(i)-273.15)
151  !  
152           ! no approximation:
153  ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))         tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) &
154  ! ori          s=1./s              /(cpd+(cl-cpd)*qnk(i))
155  ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)  
156            ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3         clw1(i, icbs1(i))=qnk(i)-qg
157            tg=tg+s*(ah0(i)-ahg)         clw1(i, icbs1(i))=max(0.0, clw1(i, icbs1(i)))
158  ! ori          tg=max(tg,35.0)  
159  !debug          tc=tg-t0         ! (qg utilise au lieu du vrai mixing ratio rg)
160            tc=tg-273.15         tvp1(i, icbs1(i))=tp1(i, icbs1(i))*(1.+qg/eps-qnk(i))
161            denom=243.5+tc  
162            denom=MAX(denom,1.0) ! convect3      end do
163  ! ori          if(tc.ge.0.0)then  
164             es=6.112*exp(17.67*tc/denom)      ! * icbs1 is the first level above the LCL:
165  ! ori          else      ! if plcl1<p1(icb1), then icbs1=icb1+1
166  ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))      ! if plcl1>p1(icb1), then icbs1=icb1
167  ! ori          end if  
168  ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))      ! * the routine above computes tvp1 from minorig to icbs1 (included).
169            qg=eps*es/(p(i,icbs(i))-es*(1.-eps))  
170        ! * to compute buoybase (in cv30_trigger.F), both tvp1(icb1) and
171           alv=lv0-clmcpv*(ticb(i)-273.15)      ! tvp1(icb1+1) must be known. This is the case if icbs1=icb1+1,
172        ! but not if icbs1=icb1.
173  ! ori c approximation here:  
174  ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)      ! * therefore, in the case icbs1=icb1, we compute tvp1 at level icb1+1
175  ! ori     &   -gz(i,icb(i))-alv*qg)/cpd      ! (tvp1 at other levels will be computed in cv30_undilute2.F)
176    
177  ! convect3: no approximation:      do i=1, klon
178           tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg) &         ticb(i)=t1(i, icb1(i)+1)
179                        /(cpd+(cl-cpd)*qnk(i))         gzicb(i)=gz1(i, icb1(i)+1)
180           qsicb(i)=qs1(i, icb1(i)+1)
181  ! ori         clw(i,icb(i))=qnk(i)-qg      enddo
182  ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))  
183           clw(i,icbs(i))=qnk(i)-qg      do i=1, klon
184           clw(i,icbs(i))=max(0.0,clw(i,icbs(i)))         tg=ticb(i)
185           qg=qsicb(i)
186  ! convect3: (qg utilise au lieu du vrai mixing ratio rg)         !debug alv=lv0-clmcpv*(ticb(i)-t0)
187           tvp(i,icbs(i))=tp(i,icbs(i))*(1.+qg/eps-qnk(i)) !whole thing         alv=lv0-clmcpv*(ticb(i)-273.15)
188    
189    360   continue         ! First iteration.
190  !  
191  ! ori      do 380 k=minorig,icbsmax2         s=cpd*(1.-qnk(i))+cl*qnk(i) &
192  ! ori       do 370 i=1,len              +alv*alv*qg/(rrv*ticb(i)*ticb(i))
193  ! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)         s=1./s
194  ! ori 370   continue  
195  ! ori 380  continue         ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
196  !         tg=tg+s*(ah0(i)-ahg)
197    
198  ! -- The following is only for convect3:         !debug tc=tg-t0
199  !         tc=tg-273.15
200  ! * icbs is the first level above the LCL:         denom=243.5+tc
201  !    if plcl<p(icb), then icbs=icb+1         denom=MAX(denom, 1.0)
202  !    if plcl>p(icb), then icbs=icb  
203  !         es=6.112*exp(17.67*tc/denom)
204  ! * the routine above computes tvp from minorig to icbs (included).  
205  !         qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))
206  ! * to compute buoybase (in cv30_trigger.F), both tvp(icb) and tvp(icb+1)  
207  !    must be known. This is the case if icbs=icb+1, but not if icbs=icb.         ! Second iteration.
208  !  
209  ! * therefore, in the case icbs=icb, we compute tvp at level icb+1         ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
210  !   (tvp at other levels will be computed in cv30_undilute2.F)         tg=tg+s*(ah0(i)-ahg)
211  !  
212           !debug tc=tg-t0
213          do i=1,len         tc=tg-273.15
214           ticb(i)=t(i,icb(i)+1)         denom=243.5+tc
215           gzicb(i)=gz(i,icb(i)+1)         denom=MAX(denom, 1.0)
216           qsicb(i)=qs(i,icb(i)+1)  
217          enddo         es=6.112*exp(17.67*tc/denom)
218    
219          do 460 i=1,len         qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))
220           tg=ticb(i)  
221           qg=qsicb(i) ! convect3         alv=lv0-clmcpv*(ticb(i)-273.15)
222  !debug         alv=lv0-clmcpv*(ticb(i)-t0)  
223           alv=lv0-clmcpv*(ticb(i)-273.15)         ! no approximation:
224  !         tp1(i, icb1(i)+1)=(ah0(i)-gz1(i, icb1(i)+1)-alv*qg) &
225  ! First iteration.              /(cpd+(cl-cpd)*qnk(i))
226  !  
227  ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))         clw1(i, icb1(i)+1)=qnk(i)-qg
228            s=cpd*(1.-qnk(i))+cl*qnk(i) &         clw1(i, icb1(i)+1)=max(0.0, clw1(i, icb1(i)+1))
             +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3  
           s=1./s  
 ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)  
           ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! 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  
 ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))  
           qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))  
 !  
 ! Second iteration.  
 !  
   
 ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))  
 ! ori          s=1./s  
 ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)  
           ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! 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          end if  
 ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))  
           qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))  
   
          alv=lv0-clmcpv*(ticb(i)-273.15)  
   
 ! ori c approximation here:  
 ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)  
 ! ori     &   -gz(i,icb(i))-alv*qg)/cpd  
   
 ! convect3: no approximation:  
          tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg) &  
                       /(cpd+(cl-cpd)*qnk(i))  
   
 ! ori         clw(i,icb(i))=qnk(i)-qg  
 ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))  
          clw(i,icb(i)+1)=qnk(i)-qg  
          clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1))  
229    
230  ! convect3: (qg utilise au lieu du vrai mixing ratio rg)         ! (qg utilise au lieu du vrai mixing ratio rg)
231           tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing         tvp1(i, icb1(i)+1)=tp1(i, icb1(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
232        end do
233    
234    460   continue    end SUBROUTINE cv30_undilute1
235    
236        return  end module cv30_undilute1_m
       end  

Legend:
Removed from v.189  
changed lines
  Added in v.196

  ViewVC Help
Powered by ViewVC 1.1.21