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

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

  ViewVC Help
Powered by ViewVC 1.1.21