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

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

  ViewVC Help
Powered by ViewVC 1.1.21