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