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