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

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

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

revision 187 by guez, Mon Mar 21 18:01:02 2016 UTC revision 205 by guez, Tue Jun 21 15:16:03 2016 UTC
# Line 4  module cv30_yield_m Line 4  module cv30_yield_m
4    
5  contains  contains
6    
7    SUBROUTINE cv30_yield(nloc, ncum, nd, na, icb, inb, delt, t, rr, u, v, gz, &    SUBROUTINE cv30_yield(icb, inb, delt, t, rr, u, v, gz, p, ph, h, hp, lv, &
8         p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, wt, water, &         cpn, th, ep, clw, m, tp, mp, qp, up, vp, wt, water, evap, b, ment, &
9         evap, b, ment, qent, uent, vent, nent, elij, sig, tv, tvp, iflag, &         qent, uent, vent, nent, elij, sig, tv, tvp, iflag, precip, VPrecip, &
10         precip, VPrecip, ft, fr, fu, fv, upwd, dnwd, dnwd0, ma, mike, tls, &         ft, fr, fu, fv, upwd, dnwd, ma, mike, tls, tps, qcondc)
11         tps, qcondc, wd)  
12        ! Tendencies, precipitation, variables of interface with other
13        ! processes, etc.
14    
15      use conema3_m, only: iflag_clw      use conema3_m, only: iflag_clw
16      use cv30_param_m, only: delta, minorig, nl, sigd      use cv30_param_m, only: minorig, nl, sigd
17      use cvthermo, only: cl, cpd, cpv, grav, rowl, rrd, rrv      use cv_thermo_m, only: rowl
18        USE dimphy, ONLY: klev, klon
19        use SUPHEC_M, only: rg, rcpd, rcw, rcpv, rd, rv
20    
21      ! inputs:      ! inputs:
     integer, intent(in):: ncum, nd, na, nloc  
     integer, intent(in):: icb(nloc), inb(nloc)  
     real, intent(in):: delt  
     real t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)  
     real sig(nloc, nd)  
     real gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)  
     real th(nloc, na), p(nloc, nd), tp(nloc, na)  
     real lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)  
     real m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)  
     real vp(nloc, na), wt(nloc, nd)  
     real water(nloc, na), evap(nloc, na), b(nloc, na)  
     real ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)  
     !ym real vent(nloc, na, na), nent(nloc, na), elij(nloc, na, na)  
     real vent(nloc, na, na), elij(nloc, na, na)  
     integer nent(nloc, na)  
     real tv(nloc, nd), tvp(nloc, nd)  
22    
23      ! input/output:      integer, intent(in):: icb(:)
24      integer iflag(nloc)  
25        integer, intent(in):: inb(:) ! (ncum)
26        ! first model level above the level of neutral buoyancy of the
27        ! parcel (1 <= inb <= nl - 1)
28    
29        real, intent(in):: delt
30        real, intent(in):: t(klon, klev), rr(klon, klev)
31        real, intent(in):: u(klon, klev), v(klon, klev)
32        real gz(klon, klev)
33        real p(klon, klev)
34        real ph(klon, klev + 1), h(klon, klev), hp(klon, klev)
35        real, intent(in):: lv(:, :) ! (klon, klev)
36    
37        real, intent(in):: cpn(:, :) ! (ncum, nl)
38        ! specific heat capacity at constant pressure of humid air, in J K-1 kg-1
39    
40        real, intent(in):: th(:, :) ! (ncum, nl)
41        real ep(klon, klev), clw(klon, klev)
42        real m(klon, klev)
43        real tp(klon, klev)
44    
45        real, intent(in):: mp(:, :) ! (ncum, nl) Mass flux of the
46        ! unsaturated downdraft, defined positive downward, in kg m-2
47        ! s-1. M_p in Emanuel (1991 928).
48    
49        real, intent(in):: qp(:, :), up(:, :) ! (klon, klev)
50        real, intent(in):: vp(:, 2:) ! (ncum, 2:nl)
51        real, intent(in):: wt(:, :) ! (ncum, nl - 1)
52        real, intent(in):: water(:, :), evap(:, :) ! (ncum, nl)
53        real, intent(in):: b(:, :) ! (ncum, nl - 1)
54        real ment(klon, klev, klev), qent(klon, klev, klev), uent(klon, klev, klev)
55        real vent(klon, klev, klev)
56        integer nent(klon, klev)
57        real elij(klon, klev, klev)
58        real sig(klon, klev)
59        real tv(klon, klev), tvp(klon, klev)
60    
61      ! outputs:      ! outputs:
62      real precip(nloc)      integer, intent(out):: iflag(:) ! (ncum)
63      real VPrecip(nloc, nd+1)      real precip(klon)
64      real ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)      real VPrecip(klon, klev + 1)
65      real upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)      real ft(klon, klev), fr(klon, klev), fu(klon, klev), fv(klon, klev)
66      real dnwd0(nloc, nd), mike(nloc, nd)      real upwd(klon, klev), dnwd(klon, klev)
67      real tls(nloc, nd), tps(nloc, nd)      real ma(klon, klev)
68      real qcondc(nloc, nd) ! cld      real mike(klon, klev)
69      real wd(nloc) ! gust      real tls(klon, klev), tps(klon, klev)
70        real qcondc(klon, klev)
71      ! local variables:  
72      integer i, k, il, n, j, num1      ! Local:
73      real rat, awat, delti      real, parameter:: delta = 0.01 ! interface cloud parameterization
74        integer ncum
75        integer i, k, il, n, j
76        real awat, delti
77      real ax, bx, cx, dx      real ax, bx, cx, dx
78      real cpinv, rdcp, dpinv      real cpinv, rdcp, dpinv
79      real lvcp(nloc, na)      real lvcp(klon, klev)
80      real am(nloc), work(nloc), ad(nloc), amp1(nloc)      real am(klon), work(klon), ad(klon), amp1(klon)
81  !!! real up1(nloc), dn1(nloc)      real up1(klon, klev, klev), dn1(klon, klev, klev)
82      real up1(nloc, nd, nd), dn1(nloc, nd, nd)      real asum(klon), bsum(klon), csum(klon), dsum(klon)
83      real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)      real qcond(klon, klev), nqcond(klon, klev), wa(klon, klev)
84      real qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld      real siga(klon, klev), sax(klon, klev), mac(klon, klev)
     real siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld  
85    
86      !-------------------------------------------------------------      !-------------------------------------------------------------
87    
88        ncum = size(icb)
89        iflag = 0
90    
91      ! initialization:      ! initialization:
92    
93      delti = 1.0/delt      delti = 1.0 / delt
94    
95      do il=1, ncum      do il = 1, ncum
96         precip(il)=0.0         precip(il) = 0.0
97         wd(il)=0.0 ! gust         VPrecip(il, klev + 1) = 0.
        VPrecip(il, nd+1)=0.  
98      enddo      enddo
99    
100      do i=1, nd      do i = 1, klev
101         do il=1, ncum         do il = 1, ncum
102            VPrecip(il, i)=0.0            VPrecip(il, i) = 0.0
103            ft(il, i)=0.0            ft(il, i) = 0.0
104            fr(il, i)=0.0            fr(il, i) = 0.0
105            fu(il, i)=0.0            fu(il, i) = 0.0
106            fv(il, i)=0.0            fv(il, i) = 0.0
107            qcondc(il, i)=0.0 ! cld            qcondc(il, i) = 0.0
108            qcond(il, i)=0.0 ! cld            qcond(il, i) = 0.0
109            nqcond(il, i)=0.0 ! cld            nqcond(il, i) = 0.0
110         enddo         enddo
111      enddo      enddo
112    
113      do i=1, nl      do i = 1, nl
114         do il=1, ncum         do il = 1, ncum
115            lvcp(il, i)=lv(il, i)/cpn(il, i)            lvcp(il, i) = lv(il, i) / cpn(il, i)
116         enddo         enddo
117      enddo      enddo
118    
119      ! *** calculate surface precipitation in mm/day ***      ! calculate surface precipitation in mm / day
120    
121      do il=1, ncum      do il = 1, ncum
122         if(ep(il, inb(il)) >= 0.0001)then         if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd &
123            precip(il)=wt(il, 1)*sigd*water(il, 1)*86400.*1000./(rowl*grav)              * water(il, 1) * 86400. * 1000. / (rowl * rg)
        endif  
124      enddo      enddo
125    
126      ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===      ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s ===
127    
128      ! MAF rajout pour lessivage      ! MAF rajout pour lessivage
129      do k=1, nl      do k = 1, nl - 1
130         do il=1, ncum         do il = 1, ncum
131            if (k <= inb(il)) then            if (k <= inb(il)) VPrecip(il, k) = wt(il, k) * sigd * water(il, k) &
132               VPrecip(il, k) = wt(il, k)*sigd*water(il, k)/grav                 / rg
           endif  
133         end do         end do
134      end do      end do
135    
136      ! *** calculate tendencies of lowest level potential temperature ***      ! calculate tendencies of lowest level potential temperature
137      ! *** and mixing ratio ***      ! and mixing ratio
138    
139      do il=1, ncum      do il = 1, ncum
140         work(il)=1.0/(ph(il, 1)-ph(il, 2))         work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
141         am(il)=0.0         am(il) = 0.0
142      enddo      enddo
143    
144      do k=2, nl      do k = 2, nl
145         do il=1, ncum         do il = 1, ncum
146            if (k <= inb(il)) then            if (k <= inb(il)) am(il) = am(il) + m(il, k)
              am(il)=am(il)+m(il, k)  
           endif  
147         enddo         enddo
148      enddo      enddo
149    
150      do il=1, ncum      do il = 1, ncum
151           if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1
        ! convect3 if((0.1*dpinv*am) >= delti)iflag(il)=4  
        if((0.01*grav*work(il)*am(il)) >= delti)iflag(il)=1!consist vect  
        ft(il, 1)=0.01*grav*work(il)*am(il)*(t(il, 2)-t(il, 1) &  
             +(gz(il, 2)-gz(il, 1))/cpn(il, 1))  
   
        ft(il, 1)=ft(il, 1)-0.5*lvcp(il, 1)*sigd*(evap(il, 1)+evap(il, 2))  
152    
153         ft(il, 1)=ft(il, 1)-0.009*grav*sigd*mp(il, 2) &         ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) &
154              *t(il, 1)*b(il, 1)*work(il)              + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) - 0.5 * lvcp(il, 1) &
155                * sigd * (evap(il, 1) + evap(il, 2)) - 0.009 * rg * sigd &
156                * mp(il, 2) * t(il, 1) * b(il, 1) * work(il) + 0.01 * sigd &
157                * wt(il, 1) * (rcw - rcpd) * water(il, 2) * (t(il, 2) - t(il, 1)) &
158                * work(il) / cpn(il, 1)
159    
160         ft(il, 1)=ft(il, 1)+0.01*sigd*wt(il, 1)*(cl-cpd)*water(il, 2)*(t(il, 2) &         ! jyg Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
161              -t(il, 1))*work(il)/cpn(il, 1)         ! (sb: pour l'instant, on ne fait que le chgt concernant rg, pas evap)
162           fr(il, 1) = 0.01 * rg * mp(il, 2) * (qp(il, 2) - rr(il, 1)) &
163                * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
164           ! + tard : + sigd * evap(il, 1)
165    
166         !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)         fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) &
167         ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)              * work(il)
        fr(il, 1)=0.01*grav*mp(il, 2)*(rp(il, 2)-rr(il, 1))*work(il) &  
             +sigd*0.5*(evap(il, 1)+evap(il, 2))  
        !+tard : +sigd*evap(il, 1)  
168    
169         fr(il, 1)=fr(il, 1)+0.01*grav*am(il)*(rr(il, 2)-rr(il, 1))*work(il)         fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
170                * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1)))
171         fu(il, 1)=fu(il, 1)+0.01*grav*work(il)*(mp(il, 2)*(up(il, 2)-u(il, 1)) &         fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
172              +am(il)*(u(il, 2)-u(il, 1)))              * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
173         fv(il, 1)=fv(il, 1)+0.01*grav*work(il)*(mp(il, 2)*(vp(il, 2)-v(il, 1)) &      enddo
             +am(il)*(v(il, 2)-v(il, 1)))  
     enddo ! il  
174    
175      do j=2, nl      do j = 2, nl
176         do il=1, ncum         do il = 1, ncum
177            if (j <= inb(il)) then            if (j <= inb(il)) then
178               fr(il, 1)=fr(il, 1) &               fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
179                    +0.01*grav*work(il)*ment(il, j, 1)*(qent(il, j, 1)-rr(il, 1))                    * (qent(il, j, 1) - rr(il, 1))
180               fu(il, 1)=fu(il, 1) &               fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
181                    +0.01*grav*work(il)*ment(il, j, 1)*(uent(il, j, 1)-u(il, 1))                    * (uent(il, j, 1) - u(il, 1))
182               fv(il, 1)=fv(il, 1) &               fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
183                    +0.01*grav*work(il)*ment(il, j, 1)*(vent(il, j, 1)-v(il, 1))                    * (vent(il, j, 1) - v(il, 1))
184            endif ! j            endif
185         enddo         enddo
186      enddo      enddo
187    
188      ! *** calculate tendencies of potential temperature and mixing ratio ***      ! calculate tendencies of potential temperature and mixing ratio
189      ! *** at levels above the lowest level ***      ! at levels above the lowest level
190    
191      ! *** first find the net saturated updraft and downdraft mass fluxes ***      ! first find the net saturated updraft and downdraft mass fluxes
192      ! *** through each level ***      ! through each level
193    
194      do i=2, nl+1 ! newvecto: mettre nl au lieu nl+1?      loop_i: do i = 2, nl - 1
195           if (any(inb >= i)) then
196              amp1(:ncum) = 0.
197              ad(:ncum) = 0.
198    
199              do k = i + 1, nl + 1
200                 do il = 1, ncum
201                    if (i <= inb(il) .and. k <= (inb(il) + 1)) then
202                       amp1(il) = amp1(il) + m(il, k)
203                    endif
204                 end do
205              end do
206    
207         num1=0            do k = 1, i
208         do il=1, ncum               do j = i + 1, nl + 1
209            if(i <= inb(il))num1=num1+1                  do il = 1, ncum
210         enddo                     if (i <= inb(il) .and. j <= (inb(il) + 1)) then
211         if(num1 <= 0) cycle                        amp1(il) = amp1(il) + ment(il, k, j)
212                       endif
213                    end do
214                 end do
215              end do
216    
217         amp1(:ncum) = 0.            do k = 1, i - 1
218         ad(:ncum) = 0.               do j = i, nl + 1 ! newvecto: nl au lieu nl + 1?
219                    do il = 1, ncum
220                       if (i <= inb(il) .and. j <= inb(il)) then
221                          ad(il) = ad(il) + ment(il, j, k)
222                       endif
223                    end do
224                 end do
225              end do
226    
227              do il = 1, ncum
228                 if (i <= inb(il)) then
229                    dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
230                    cpinv = 1.0 / cpn(il, i)
231    
232                    if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 1
233    
234         do k=i+1, nl+1                  ft(il, i) = 0.01 * rg * dpinv * (amp1(il) * (t(il, i + 1) &
235            do il=1, ncum                       - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) &
236               if (i <= inb(il) .and. k <= (inb(il)+1)) then                       - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) &
237                  amp1(il)=amp1(il)+m(il, k)                       - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) &
238                         * (evap(il, i) + evap(il, i + 1)) - 0.009 * rg * sigd &
239                         * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il, i) &
240                         * t(il, i - 1) * cpn(il, i - 1) * cpinv * b(il, i - 1)) &
241                         * dpinv + 0.01 * rg * dpinv * ment(il, i, i) &
242                         * (hp(il, i) - h(il, i) + t(il, i) * (rcpv - rcpd) &
243                         * (rr(il, i) - qent(il, i, i))) * cpinv + 0.01 * sigd &
244                         * wt(il, i) * (rcw - rcpd) * water(il, i + 1) &
245                         * (t(il, i + 1) - t(il, i)) * dpinv * cpinv
246                    fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) &
247                         - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
248                    fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) &
249                         * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
250                         - u(il, i - 1)))
251                    fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) &
252                         * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
253                         - v(il, i - 1)))
254               endif               endif
255            end do            end do
        end do  
256    
257         do k=1, i            do k = 1, i - 1
258            do j=i+1, nl+1               do il = 1, ncum
259               do il=1, ncum                  if (i <= inb(il)) then
260                  if (i <= inb(il) .and. j <= (inb(il)+1)) then                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
261                     amp1(il)=amp1(il)+ment(il, k, j)                     cpinv = 1.0 / cpn(il, i)
262                  endif  
263                       awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
264                       awat = amax1(awat, 0.0)
265    
266                       fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
267                            * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
268                       fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
269                            * ment(il, k, i) * (uent(il, k, i) - u(il, i))
270                       fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
271                            * ment(il, k, i) * (vent(il, k, i) - v(il, i))
272    
273                       ! (saturated updrafts resulting from mixing)
274                       qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat)
275                       nqcond(il, i) = nqcond(il, i) + 1.
276                    endif ! i
277               end do               end do
278            end do            end do
        end do  
279    
280         do k=1, i-1            do k = i, nl + 1
281            do j=i, nl+1 ! newvecto: nl au lieu nl+1?               do il = 1, ncum
282               do il=1, ncum                  if (i <= inb(il) .and. k <= inb(il)) then
283                  if (i <= inb(il) .and. j <= inb(il)) then                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
284                     ad(il)=ad(il)+ment(il, j, k)                     cpinv = 1.0 / cpn(il, i)
285    
286                       fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
287                            * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
288                       fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
289                            * ment(il, k, i) * (uent(il, k, i) - u(il, i))
290                       fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
291                            * ment(il, k, i) * (vent(il, k, i) - v(il, i))
292                  endif                  endif
293               end do               end do
294            end do            end do
        end do  
295    
296         do il=1, ncum            do il = 1, ncum
           if (i <= inb(il)) then  
              dpinv=1.0/(ph(il, i)-ph(il, i+1))  
              cpinv=1.0/cpn(il, i)  
   
              if((0.01*grav*dpinv*amp1(il)) >= delti)iflag(il)=1 ! vecto  
   
              ft(il, i)=0.01*grav*dpinv*(amp1(il)*(t(il, i+1)-t(il, i) &  
                   +(gz(il, i+1)-gz(il, i))*cpinv) &  
                   -ad(il)*(t(il, i)-t(il, i-1)+(gz(il, i)-gz(il, i-1))*cpinv)) &  
                   -0.5*sigd*lvcp(il, i)*(evap(il, i)+evap(il, i+1))  
              rat=cpn(il, i-1)*cpinv  
              ft(il, i)=ft(il, i)-0.009*grav*sigd*(mp(il, i+1)*t(il, i)*b(il, i) &  
                   -mp(il, i)*t(il, i-1)*rat*b(il, i-1))*dpinv  
              ft(il, i)=ft(il, i)+0.01*grav*dpinv*ment(il, i, i)*(hp(il, i)-h(il, i) &  
                   +t(il, i)*(cpv-cpd)*(rr(il, i)-qent(il, i, i)))*cpinv  
   
              ft(il, i)=ft(il, i)+0.01*sigd*wt(il, i)*(cl-cpd)*water(il, i+1) &  
                   *(t(il, i+1)-t(il, i))*dpinv*cpinv  
   
              fr(il, i)=0.01*grav*dpinv*(amp1(il)*(rr(il, i+1)-rr(il, i)) &  
                   -ad(il)*(rr(il, i)-rr(il, i-1)))  
              fu(il, i)=fu(il, i)+0.01*grav*dpinv*(amp1(il)*(u(il, i+1)-u(il, i)) &  
                   -ad(il)*(u(il, i)-u(il, i-1)))  
              fv(il, i)=fv(il, i)+0.01*grav*dpinv*(amp1(il)*(v(il, i+1)-v(il, i)) &  
                   -ad(il)*(v(il, i)-v(il, i-1)))  
           endif ! i  
        end do  
   
        do k=1, i-1  
           do il=1, ncum  
297               if (i <= inb(il)) then               if (i <= inb(il)) then
298                  dpinv=1.0/(ph(il, i)-ph(il, i+1))                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
299                  cpinv=1.0/cpn(il, i)                  cpinv = 1.0 / cpn(il, i)
   
                 awat=elij(il, k, i)-(1.-ep(il, i))*clw(il, i)  
                 awat=amax1(awat, 0.0)  
   
                 fr(il, i)=fr(il, i) &  
                      +0.01*grav*dpinv*ment(il, k, i)*(qent(il, k, i)-awat-rr(il, i))  
                 fu(il, i)=fu(il, i) &  
                      +0.01*grav*dpinv*ment(il, k, i)*(uent(il, k, i)-u(il, i))  
                 fv(il, i)=fv(il, i) &  
                      +0.01*grav*dpinv*ment(il, k, i)*(vent(il, k, i)-v(il, i))  
   
                 ! (saturated updrafts resulting from mixing) ! cld  
                 qcond(il, i)=qcond(il, i)+(elij(il, k, i)-awat) ! cld  
                 nqcond(il, i)=nqcond(il, i)+1. ! cld  
              endif ! i  
           end do  
        end do  
300    
301         do k=i, nl+1                  ! sb: on ne fait pas encore la correction permettant de mieux
302            do il=1, ncum                  ! conserver l'eau:
303               if (i <= inb(il) .and. k <= inb(il)) then                  fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) &
304                  dpinv=1.0/(ph(il, i)-ph(il, i+1))                       + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) &
305                  cpinv=1.0/cpn(il, i)                       * (qp(il, i + 1) - rr(il, i)) - mp(il, i) * (qp(il, i) &
306                         - rr(il, i - 1))) * dpinv
307                  fr(il, i)=fr(il, i) &  
308                       +0.01*grav*dpinv*ment(il, k, i)*(qent(il, k, i)-rr(il, i))                  fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) &
309                  fu(il, i)=fu(il, i) &                       * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
310                       +0.01*grav*dpinv*ment(il, k, i)*(uent(il, k, i)-u(il, i))                       - u(il, i - 1))) * dpinv
311                  fv(il, i)=fv(il, i) &                  fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) &
312                       +0.01*grav*dpinv*ment(il, k, i)*(vent(il, k, i)-v(il, i))                       * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) &
313               endif ! i and k                       - v(il, i - 1))) * dpinv
314                 endif
315            end do            end do
        end do  
316    
317         do il=1, ncum            ! sb: interface with the cloud parameterization:
           if (i <= inb(il)) then  
              dpinv=1.0/(ph(il, i)-ph(il, i+1))  
              cpinv=1.0/cpn(il, i)  
   
              ! sb: on ne fait pas encore la correction permettant de mieux  
              ! conserver l'eau:  
              fr(il, i)=fr(il, i)+0.5*sigd*(evap(il, i)+evap(il, i+1)) &  
                   +0.01*grav*(mp(il, i+1)*(rp(il, i+1)-rr(il, i))-mp(il, i) &  
                   *(rp(il, i)-rr(il, i-1)))*dpinv  
   
              fu(il, i)=fu(il, i)+0.01*grav*(mp(il, i+1)*(up(il, i+1)-u(il, i)) &  
                   -mp(il, i)*(up(il, i)-u(il, i-1)))*dpinv  
              fv(il, i)=fv(il, i)+0.01*grav*(mp(il, i+1)*(vp(il, i+1)-v(il, i)) &  
                   -mp(il, i)*(vp(il, i)-v(il, i-1)))*dpinv  
   
           endif ! i  
        end do  
318    
319         ! sb: interface with the cloud parameterization: ! cld            do k = i + 1, nl
320                 do il = 1, ncum
321                    if (k <= inb(il) .and. i <= inb(il)) then
322                       ! (saturated downdrafts resulting from mixing)
323                       qcond(il, i) = qcond(il, i) + elij(il, k, i)
324                       nqcond(il, i) = nqcond(il, i) + 1.
325                    endif
326                 enddo
327              enddo
328    
329         do k=i+1, nl            ! (particular case: no detraining level is found)
330            do il=1, ncum            do il = 1, ncum
331               if (k <= inb(il) .and. i <= inb(il)) then ! cld               if (i <= inb(il) .and. nent(il, i) == 0) then
332                  ! (saturated downdrafts resulting from mixing) ! cld                  qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i)
333                  qcond(il, i)=qcond(il, i)+elij(il, k, i) ! cld                  nqcond(il, i) = nqcond(il, i) + 1.
334                  nqcond(il, i)=nqcond(il, i)+1. ! cld               endif
335               endif ! cld            enddo
           enddo ! cld  
        enddo ! cld  
   
        ! (particular case: no detraining level is found) ! cld  
        do il=1, ncum ! cld  
           if (i <= inb(il) .and. nent(il, i) == 0) then ! cld  
              qcond(il, i)=qcond(il, i)+(1.-ep(il, i))*clw(il, i) ! cld  
              nqcond(il, i)=nqcond(il, i)+1. ! cld  
           endif ! cld  
        enddo ! cld  
   
        do il=1, ncum ! cld  
           if (i <= inb(il) .and. nqcond(il, i) /= 0.) then ! cld  
              qcond(il, i)=qcond(il, i)/nqcond(il, i) ! cld  
           endif ! cld  
        enddo  
336    
337      end do            do il = 1, ncum
338                 if (i <= inb(il) .and. nqcond(il, i) /= 0.) then
339                    qcond(il, i) = qcond(il, i) / nqcond(il, i)
340                 endif
341              enddo
342           end if
343        end do loop_i
344    
345      ! *** move the detrainment at level inb down to level inb-1 ***      ! move the detrainment at level inb down to level inb - 1
346      ! *** in such a way as to preserve the vertically ***      ! in such a way as to preserve the vertically
347      ! *** integrated enthalpy and water tendencies ***      ! integrated enthalpy and water tendencies
348    
349      do il=1, ncum      do il = 1, ncum
350           ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) &
351         ax=0.1*ment(il, inb(il), inb(il))*(hp(il, inb(il))-h(il, inb(il)) &              - h(il, inb(il)) + t(il, inb(il)) * (rcpv - rcpd) &
352              +t(il, inb(il))*(cpv-cpd) &              * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) &
353              *(rr(il, inb(il))-qent(il, inb(il), inb(il)))) &              / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
354              /(cpn(il, inb(il))*(ph(il, inb(il))-ph(il, inb(il)+1)))         ft(il, inb(il)) = ft(il, inb(il)) - ax
355         ft(il, inb(il))=ft(il, inb(il))-ax         ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) &
356         ft(il, inb(il)-1)=ft(il, inb(il)-1)+ax*cpn(il, inb(il)) &              * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) &
357              *(ph(il, inb(il))-ph(il, inb(il)+1))/(cpn(il, inb(il)-1) &              * (ph(il, inb(il) - 1) - ph(il, inb(il))))
358              *(ph(il, inb(il)-1)-ph(il, inb(il))))  
359           bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) &
360         bx=0.1*ment(il, inb(il), inb(il))*(qent(il, inb(il), inb(il)) &              - rr(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
361              -rr(il, inb(il)))/(ph(il, inb(il))-ph(il, inb(il)+1))         fr(il, inb(il)) = fr(il, inb(il)) - bx
362         fr(il, inb(il))=fr(il, inb(il))-bx         fr(il, inb(il) - 1) = fr(il, inb(il) - 1) &
363         fr(il, inb(il)-1)=fr(il, inb(il)-1) &              + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
364              +bx*(ph(il, inb(il))-ph(il, inb(il)+1)) &              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
365              /(ph(il, inb(il)-1)-ph(il, inb(il)))  
366           cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) &
367         cx=0.1*ment(il, inb(il), inb(il))*(uent(il, inb(il), inb(il)) &              - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
368              -u(il, inb(il)))/(ph(il, inb(il))-ph(il, inb(il)+1))         fu(il, inb(il)) = fu(il, inb(il)) - cx
369         fu(il, inb(il))=fu(il, inb(il))-cx         fu(il, inb(il) - 1) = fu(il, inb(il) - 1) &
370         fu(il, inb(il)-1)=fu(il, inb(il)-1) &              + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
371              +cx*(ph(il, inb(il))-ph(il, inb(il)+1)) &              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
372              /(ph(il, inb(il)-1)-ph(il, inb(il)))  
373           dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) &
374         dx=0.1*ment(il, inb(il), inb(il))*(vent(il, inb(il), inb(il)) &              - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
375              -v(il, inb(il)))/(ph(il, inb(il))-ph(il, inb(il)+1))         fv(il, inb(il)) = fv(il, inb(il)) - dx
376         fv(il, inb(il))=fv(il, inb(il))-dx         fv(il, inb(il) - 1) = fv(il, inb(il) - 1) &
377         fv(il, inb(il)-1)=fv(il, inb(il)-1) &              + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
378              +dx*(ph(il, inb(il))-ph(il, inb(il)+1)) &              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
             /(ph(il, inb(il)-1)-ph(il, inb(il)))  
379    
380      end do      end do
381    
382      ! *** homoginize tendencies below cloud base ***      ! homoginize tendencies below cloud base
383    
384      do il=1, ncum      do il = 1, ncum
385         asum(il)=0.0         asum(il) = 0.0
386         bsum(il)=0.0         bsum(il) = 0.0
387         csum(il)=0.0         csum(il) = 0.0
388         dsum(il)=0.0         dsum(il) = 0.0
389      enddo      enddo
390    
391      do i=1, nl      do i = 1, nl
392         do il=1, ncum         do il = 1, ncum
393            if (i <= (icb(il)-1)) then            if (i <= (icb(il) - 1)) then
394               asum(il)=asum(il)+ft(il, i)*(ph(il, i)-ph(il, i+1))               asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
395               bsum(il)=bsum(il)+fr(il, i)*(lv(il, i)+(cl-cpd)*(t(il, i)-t(il, 1))) &               bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (rcw - rcpd) &
396                    *(ph(il, i)-ph(il, i+1))                    * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
397               csum(il)=csum(il)+(lv(il, i)+(cl-cpd)*(t(il, i)-t(il, 1))) &               csum(il) = csum(il) + (lv(il, i) + (rcw - rcpd) * (t(il, i) &
398                    *(ph(il, i)-ph(il, i+1))                    - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
399               dsum(il)=dsum(il)+t(il, i)*(ph(il, i)-ph(il, i+1))/th(il, i)               dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) &
400                      / th(il, i)
401            endif            endif
402         enddo         enddo
403      enddo      enddo
404    
405      do i=1, nl      do i = 1, nl
406         do il=1, ncum         do il = 1, ncum
407            if (i <= (icb(il)-1)) then            if (i <= (icb(il) - 1)) then
408               ft(il, i)=asum(il)*t(il, i)/(th(il, i)*dsum(il))               ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
409               fr(il, i)=bsum(il)/csum(il)               fr(il, i) = bsum(il) / csum(il)
410            endif            endif
411         enddo         enddo
412      enddo      enddo
413    
414      ! *** reset counter and return ***      ! reset counter and return
415    
416      do il=1, ncum      do il = 1, ncum
417         sig(il, nd)=2.0         sig(il, klev) = 2.0
418      enddo      enddo
419    
420      do i=1, nd      do i = 1, klev
421         do il=1, ncum         do il = 1, ncum
422            upwd(il, i)=0.0            upwd(il, i) = 0.0
423            dnwd(il, i)=0.0            dnwd(il, i) = 0.0
424         enddo         enddo
425      enddo      enddo
426    
427      do i=1, nl      do i = 1, nl
428         do il=1, ncum         do il = 1, ncum
           dnwd0(il, i)=-mp(il, i)  
        enddo  
     enddo  
     do i=nl+1, nd  
        do il=1, ncum  
           dnwd0(il, i)=0.  
        enddo  
     enddo  
   
     do i=1, nl  
        do il=1, ncum  
429            if (i >= icb(il) .and. i <= inb(il)) then            if (i >= icb(il) .and. i <= inb(il)) then
430               upwd(il, i)=0.0               upwd(il, i) = 0.0
431               dnwd(il, i)=0.0               dnwd(il, i) = 0.0
432            endif            endif
433         enddo         enddo
434      enddo      enddo
435    
436      do i=1, nl      do i = 1, nl
437         do k=1, nl         do k = 1, nl
438            do il=1, ncum            do il = 1, ncum
439               up1(il, k, i)=0.0               up1(il, k, i) = 0.0
440               dn1(il, k, i)=0.0               dn1(il, k, i) = 0.0
441            enddo            enddo
442         enddo         enddo
443      enddo      enddo
444    
445      do i=1, nl      do i = 1, nl
446         do k=i, nl         do k = i, nl
447            do n=1, i-1            do n = 1, i - 1
448               do il=1, ncum               do il = 1, ncum
449                  if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then                  if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
450                     up1(il, k, i)=up1(il, k, i)+ment(il, n, k)                     up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
451                     dn1(il, k, i)=dn1(il, k, i)-ment(il, k, n)                     dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
452                  endif                  endif
453               enddo               enddo
454            enddo            enddo
455         enddo         enddo
456      enddo      enddo
457    
458      do i=2, nl      do i = 2, nl
459         do k=i, nl         do k = i, nl
460            do il=1, ncum            do il = 1, ncum
461               if (i <= inb(il).and.k <= inb(il)) then               if (i <= inb(il).and.k <= inb(il)) then
462                  upwd(il, i)=upwd(il, i)+m(il, k)+up1(il, k, i)                  upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
463                  dnwd(il, i)=dnwd(il, i)+dn1(il, k, i)                  dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
464               endif               endif
465            enddo            enddo
466         enddo         enddo
467      enddo      enddo
468    
469      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      ! D\'etermination de la variation de flux ascendant entre
470      ! determination de la variation de flux ascendant entre      ! deux niveaux non dilu\'es mike
     ! deux niveau non dilue mike  
     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
471    
472      do i=1, nl      do i = 1, nl
473         do il=1, ncum         do il = 1, ncum
474            mike(il, i)=m(il, i)            mike(il, i) = m(il, i)
475         enddo         enddo
476      enddo      enddo
477    
478      do i=nl+1, nd      do i = nl + 1, klev
479         do il=1, ncum         do il = 1, ncum
480            mike(il, i)=0.            mike(il, i) = 0.
481         enddo         enddo
482      enddo      enddo
483    
484      do i=1, nd      do i = 1, klev
485         do il=1, ncum         do il = 1, ncum
486            ma(il, i)=0            ma(il, i) = 0
487         enddo         enddo
488      enddo      enddo
489    
490      do i=1, nl      do i = 1, nl
491         do j=i, nl         do j = i, nl
492            do il=1, ncum            do il = 1, ncum
493               ma(il, i)=ma(il, i)+m(il, j)               ma(il, i) = ma(il, i) + m(il, j)
494            enddo            enddo
495         enddo         enddo
496      enddo      enddo
497    
498      do i=nl+1, nd      do i = nl + 1, klev
499         do il=1, ncum         do il = 1, ncum
500            ma(il, i)=0.            ma(il, i) = 0.
501         enddo         enddo
502      enddo      enddo
503    
504      do i=1, nl      do i = 1, nl
505         do il=1, ncum         do il = 1, ncum
506            if (i <= (icb(il)-1)) then            if (i <= (icb(il) - 1)) then
507               ma(il, i)=0               ma(il, i) = 0
508            endif            endif
509         enddo         enddo
510      enddo      enddo
511    
512      !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      ! icb repr\'esente le niveau o\`u se trouve la base du nuage, et
513      ! icb represente de niveau ou se trouve la      ! inb le sommet du nuage
514      ! base du nuage, et inb le top du nuage  
515      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      do i = 1, klev
516           DO il = 1, ncum
517      do i=1, nd            rdcp = (rd * (1. - rr(il, i)) - rr(il, i) * rv) &
518         DO il=1, ncum                 / (rcpd * (1. - rr(il, i)) + rr(il, i) * rcpv)
519            rdcp=(rrd*(1.-rr(il, i))-rr(il, i)*rrv) &            tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
520                 /(cpd*(1.-rr(il, i))+rr(il, i)*cpv)            tps(il, i) = tp(il, i)
           tls(il, i)=t(il, i)*(1000.0/p(il, i))**rdcp  
           tps(il, i)=tp(il, i)  
521         end DO         end DO
522      enddo      enddo
523    
524      ! *** diagnose the in-cloud mixing ratio *** ! cld      ! Diagnose the in-cloud mixing ratio of condensed water
525      ! *** of condensed water *** ! cld  
526      ! ! cld      do i = 1, klev
527           do il = 1, ncum
528      do i=1, nd ! cld            mac(il, i) = 0.0
529         do il=1, ncum ! cld            wa(il, i) = 0.0
530            mac(il, i)=0.0 ! cld            siga(il, i) = 0.0
531            wa(il, i)=0.0 ! cld            sax(il, i) = 0.0
532            siga(il, i)=0.0 ! cld         enddo
533            sax(il, i)=0.0 ! cld      enddo
534         enddo ! cld  
535      enddo ! cld      do i = minorig, nl
536           do k = i + 1, nl + 1
537      do i=minorig, nl ! cld            do il = 1, ncum
538         do k=i+1, nl+1 ! cld               if (i <= inb(il) .and. k <= (inb(il) + 1)) then
539            do il=1, ncum ! cld                  mac(il, i) = mac(il, i) + m(il, k)
540               if (i <= inb(il) .and. k <= (inb(il)+1)) then ! cld               endif
541                  mac(il, i)=mac(il, i)+m(il, k) ! cld            enddo
542               endif ! cld         enddo
543            enddo ! cld      enddo
544         enddo ! cld  
545      enddo ! cld      do i = 1, nl
546           do j = 1, i
547      do i=1, nl ! cld            do il = 1, ncum
548         do j=1, i ! cld               if (i >= icb(il) .and. i <= (inb(il) - 1) &
549            do il=1, ncum ! cld                    .and. j >= icb(il)) then
550               if (i >= icb(il) .and. i <= (inb(il)-1) &                  sax(il, i) = sax(il, i) + rd * (tvp(il, j) - tv(il, j)) &
551                    .and. j >= icb(il)) then ! cld                       * (ph(il, j) - ph(il, j + 1)) / p(il, j)
552                  sax(il, i)=sax(il, i)+rrd*(tvp(il, j)-tv(il, j)) &               endif
553                       *(ph(il, j)-ph(il, j+1))/p(il, j) ! cld            enddo
554               endif ! cld         enddo
555            enddo ! cld      enddo
556         enddo ! cld  
557      enddo ! cld      do i = 1, nl
558           do il = 1, ncum
559      do i=1, nl ! cld            if (i >= icb(il) .and. i <= (inb(il) - 1) &
560         do il=1, ncum ! cld                 .and. sax(il, i) > 0.0) then
561            if (i >= icb(il) .and. i <= (inb(il)-1) &               wa(il, i) = sqrt(2. * sax(il, i))
562                 .and. sax(il, i) > 0.0) then ! cld            endif
563               wa(il, i)=sqrt(2.*sax(il, i)) ! cld         enddo
564            endif ! cld      enddo
565         enddo ! cld  
566      enddo ! cld      do i = 1, nl
567           do il = 1, ncum
568      do i=1, nl ! cld            if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rd &
569         do il=1, ncum ! cld                 * tvp(il, i) / p(il, i) / 100. / delta
570            if (wa(il, i) > 0.0) &            siga(il, i) = min(siga(il, i), 1.0)
571                 siga(il, i)=mac(il, i)/wa(il, i) &  
                *rrd*tvp(il, i)/p(il, i)/100./delta ! cld  
           siga(il, i) = min(siga(il, i), 1.0) ! cld  
           !IM cf. FH  
572            if (iflag_clw == 0) then            if (iflag_clw == 0) then
573               qcondc(il, i)=siga(il, i)*clw(il, i)*(1.-ep(il, i)) &               qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) &
574                    + (1.-siga(il, i))*qcond(il, i) ! cld                    + (1. - siga(il, i)) * qcond(il, i)
575            else if (iflag_clw == 1) then            else if (iflag_clw == 1) then
576               qcondc(il, i)=qcond(il, i) ! cld               qcondc(il, i) = qcond(il, i)
577            endif            endif
578           enddo
579         enddo ! cld      enddo
     enddo ! cld  
580    
581    end SUBROUTINE cv30_yield    end SUBROUTINE cv30_yield
582    
Legend:
Removed from v.187  
changed lines
  Added in v.205
  ViewVC Help
Powered by ViewVC 1.1.21