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

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

  ViewVC Help
Powered by ViewVC 1.1.21