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

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

  ViewVC Help
Powered by ViewVC 1.1.21