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

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

  ViewVC Help
Powered by ViewVC 1.1.21