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

Legend:
Removed from v.190  
changed lines
  Added in v.195

  ViewVC Help
Powered by ViewVC 1.1.21