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

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

  ViewVC Help
Powered by ViewVC 1.1.21