/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_unsat.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_unsat.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 192 by guez, Thu May 12 13:00:07 2016 UTC
# Line 4  module cv30_unsat_m Line 4  module cv30_unsat_m
4    
5  contains  contains
6    
7    SUBROUTINE cv30_unsat(ncum, icb, inb, t, rr, rs, gz, u, v, p, ph, th, tv, &    SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &
8         lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, up, vp, wt, &         ep, sigp, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, &
9         water, evap, b)         evap, b)
10    
11      use cv30_param_m, only: nl, sigd      use cv30_param_m, only: nl, sigd
12      use cvthermo, only: cpd, ginv, grav      use cv_thermo_m, only: cpd, ginv, grav
13      USE dimphy, ONLY: klon, klev      USE dimphy, ONLY: klon, klev
14    
15      ! inputs:      integer, intent(in):: icb(:) ! (ncum)
16      integer, intent(in):: ncum  
17      integer, intent(in):: icb(:), inb(:) ! (ncum)      integer, intent(in):: inb(:) ! (ncum)
18      real t(klon, klev), rr(klon, klev), rs(klon, klev)      ! first model level above the level of neutral buoyancy of the
19      real gz(klon, klev)      ! parcel (1 <= inb <= nl - 1)
20      real u(klon, klev), v(klon, klev)  
21      real p(klon, klev), ph(klon, klev + 1)      real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (klon, klev)
22      real th(klon, klev)      real, intent(in):: gz(:, :) ! (klon, klev)
23      real tv(klon, klev)      real, intent(in):: u(:, :), v(:, :) ! (klon, klev)
24      real lv(klon, klev)      real, intent(in):: p(klon, klev), ph(klon, klev + 1)
25      real cpn(klon, klev)      real, intent(in):: th(klon, klev)
26      real, intent(in):: ep(klon, klev), sigp(klon, klev), clw(klon, klev)      real, intent(in):: tv(klon, klev)
27      real m(klon, klev), ment(klon, klev, klev), elij(klon, klev, klev)      real, intent(in):: lv(klon, klev)
28        real, intent(in):: cpn(klon, klev)
29        real, intent(in):: ep(:, :), sigp(:, :), clw(:, :) ! (ncum, klev)
30        real, intent(in):: m(:, :) ! (ncum, klev)
31        real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)
32        real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)
33      real, intent(in):: delt      real, intent(in):: delt
34      real plcl(klon)      real, intent(in):: plcl(klon)
35    
36      ! outputs:      ! outputs:
37      real mp(klon, klev), rp(klon, klev), up(klon, klev), vp(klon, klev)      real, intent(out):: mp(klon, klev)
38        real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
39      real wt(klon, klev), water(klon, klev), evap(klon, klev)      real wt(klon, klev), water(klon, klev), evap(klon, klev)
40      real b(:, :) ! (ncum, klev)      real, intent(out):: b(:, :) ! (ncum, nl - 1)
41    
42      ! Local:      ! Local:
43        integer ncum
44      integer i, j, il, num1      integer i, j, il, num1
45      real tinv, delti      real tinv, delti
46      real awat, afac, afac1, afac2, bfac      real awat, afac, afac1, afac2, bfac
# Line 41  contains Line 48  contains
48      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
49      real ampmax      real ampmax
50      real lvcp(klon, klev)      real lvcp(klon, klev)
51      real wdtrain(ncum)      real wdtrain(size(icb))
52      logical lwork(ncum)      logical lwork(size(icb))
53    
54      !------------------------------------------------------      !------------------------------------------------------
55    
56        ncum = size(icb)
57      delti = 1. / delt      delti = 1. / delt
58      tinv = 1. / 3.      tinv = 1. / 3.
59      mp = 0.      mp = 0.
60        b = 0.
61    
62      do i = 1, nl      do i = 1, nl
63         do il = 1, ncum         do il = 1, ncum
64            mp(il, i) = 0.            qp(il, i) = q(il, i)
           rp(il, i) = rr(il, i)  
65            up(il, i) = u(il, i)            up(il, i) = u(il, i)
66            vp(il, i) = v(il, i)            vp(il, i) = v(il, i)
67            wt(il, i) = 0.001            wt(il, i) = 0.001
68            water(il, i) = 0.            water(il, i) = 0.
69            evap(il, i) = 0.            evap(il, i) = 0.
           b(il, i) = 0.  
70            lvcp(il, i) = lv(il, i) / cpn(il, i)            lvcp(il, i) = lv(il, i) / cpn(il, i)
71         enddo         enddo
72      enddo      enddo
# Line 94  contains Line 101  contains
101                  do il = 1, ncum                  do il = 1, ncum
102                     if (i <= inb(il) .and. lwork(il)) then                     if (i <= inb(il) .and. lwork(il)) then
103                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
104                        awat = amax1(awat, 0.)                        awat = max(awat, 0.)
105                        wdtrain(il) = wdtrain(il) + grav * awat &                        wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
                            * ment(il, j, i)  
106                     endif                     endif
107                  enddo                  enddo
108               end do               end do
109            endif            endif
110    
111            ! find rain water and evaporation using provisional            ! find rain water and evaporation using provisional
112            ! estimates of rp(i)and rp(i - 1)            ! estimates of qp(i) and qp(i - 1)
113    
114            do il = 1, ncum            do il = 1, ncum
115               if (i <= inb(il) .and. lwork(il)) then               if (i <= inb(il) .and. lwork(il)) then
116                  wt(il, i) = 45.                  wt(il, i) = 45.
117    
118                  if (i < inb(il)) then                  if (i < inb(il)) then
119                     rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) &                     qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) &
120                          - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)                          - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
121                     rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))                     qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
122                  endif                  endif
123                  rp(il, i) = amax1(rp(il, i), 0.)  
124                  rp(il, i) = amin1(rp(il, i), rs(il, i))                  qp(il, i) = max(qp(il, i), 0.)
125                  rp(il, inb(il)) = rr(il, inb(il))                  qp(il, i) = min(qp(il, i), qs(il, i))
126                    qp(il, inb(il)) = q(il, inb(il))
127    
128                  if (i == 1) then                  if (i == 1) then
129                     afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) &                     afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
130                          / (1e4 + 2000. * p(il, 1) * rs(il, 1))                          / (1e4 + 2000. * p(il, 1) * qs(il, 1))
131                  else                  else
132                     rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &                     qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) &
133                          - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)                          - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
134                     rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))                     qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
135                     rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))                     qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
136                     rp(il, i - 1) = amax1(rp(il, i - 1), 0.)                     qp(il, i - 1) = max(qp(il, i - 1), 0.)
137                     afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &                     afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &
138                          / (1e4 + 2000. * p(il, i) * rs(il, i))                          / (1e4 + 2000. * p(il, i) * qs(il, i))
139                     afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &                     afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) &
140                          / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))                          / (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1))
141                     afac = 0.5 * (afac1 + afac2)                     afac = 0.5 * (afac1 + afac2)
142                  endif                  endif
143                  if (i == inb(il))afac = 0.  
144                  afac = amax1(afac, 0.)                  if (i == inb(il)) afac = 0.
145                    afac = max(afac, 0.)
146                  bfac = 1. / (sigd * wt(il, i))                  bfac = 1. / (sigd * wt(il, i))
147    
148                  ! prise en compte de la variation progressive de sigt dans                  ! prise en compte de la variation progressive de sigt dans
# Line 168  contains Line 176  contains
176                  ! hydrostatic approximation                  ! hydrostatic approximation
177    
178                  if (i /= 1) then                  if (i /= 1) then
179                     tevap = amax1(0., evap(il, i))                     tevap = max(0., evap(il, i))
180                     delth = amax1(0.001, (th(il, i) - th(il, i - 1)))                     delth = max(0.001, (th(il, i) - th(il, i - 1)))
181                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
182                          * (p(il, i - 1) - p(il, i)) / delth                          * (p(il, i - 1) - p(il, i)) / delth
183    
184                     ! if hydrostatic assumption fails,                     ! If hydrostatic assumption fails, solve cubic
185                     ! solve cubic difference equation for downdraft theta                     ! difference equation for downdraft theta and mass
186                     ! and mass flux from two simultaneous differential eqns                     ! flux from two simultaneous differential equations
187    
188                     amfac = sigd * sigd * 70. * ph(il, i) &                     amfac = sigd * sigd * 70. * ph(il, i) &
189                          * (p(il, i - 1) - p(il, i)) &                          * (p(il, i - 1) - p(il, i)) &
190                          * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))                          * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
191                     amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &                     amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
192                          * mp(il, i))                          * mp(il, i))
193                     if (amp2 > (0.1 * amfac)) then  
194                       if (amp2 > 0.1 * amfac) then
195                        xf = 100. * sigd * sigd * sigd * (ph(il, i) &                        xf = 100. * sigd * sigd * sigd * (ph(il, i) &
196                             - ph(il, i + 1))                             - ph(il, i + 1))
197                        tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &                        tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
# Line 192  contains Line 201  contains
201                             * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &                             * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
202                             - p(il, i)) * xf * tevap                             - p(il, i)) * xf * tevap
203                        fac2 = 1.                        fac2 = 1.
204                        if (bf < 0.)fac2 = - 1.                        if (bf < 0.) fac2 = - 1.
205                        bf = abs(bf)                        bf = abs(bf)
206                        ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv                        ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
207    
208                        if (ur >= 0.) then                        if (ur >= 0.) then
209                           sru = sqrt(ur)                           sru = sqrt(ur)
210                           fac = 1.                           fac = 1.
211                           if ((0.5 * bf - sru) < 0.)fac = - 1.                           if ((0.5 * bf - sru) < 0.) fac = - 1.
212                           mp(il, i) = mp(il, i + 1) * tinv &                           mp(il, i) = mp(il, i + 1) * tinv &
213                                + (0.5 * bf + sru)**tinv &                                + (0.5 * bf + sru)**tinv &
214                                + fac * (abs(0.5 * bf - sru))**tinv                                + fac * (abs(0.5 * bf - sru))**tinv
# Line 208  contains Line 218  contains
218                           mp(il, i) = mp(il, i + 1) * tinv + 2. &                           mp(il, i) = mp(il, i + 1) * tinv + 2. &
219                                * sqrt(af * tinv) * cos(d * tinv)                                * sqrt(af * tinv) * cos(d * tinv)
220                        endif                        endif
221                        mp(il, i) = amax1(0., mp(il, i))  
222                          mp(il, i) = max(0., mp(il, i))
223    
224                        ! Il y a vraisemblablement une erreur dans la                        ! Il y a vraisemblablement une erreur dans la
225                        ! ligne 2 suivante: il faut diviser par (mp(il,                        ! ligne suivante : il faut diviser par (mp(il,
226                        ! i) * sigd * grav) et non par (mp(il, i) + sigd                        ! i) * sigd * grav) et non par (mp(il, i) + sigd
227                        ! * 0.1).  Et il faut bien revoir les facteurs                        ! * 0.1).  Et il faut bien revoir les facteurs
228                        ! 100.                        ! 100.
# Line 219  contains Line 230  contains
230                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
231                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
232                             / (lvcp(il, i) * sigd * th(il, i))                             / (lvcp(il, i) * sigd * th(il, i))
233                        b(il, i - 1) = amax1(b(il, i - 1), 0.)                        b(il, i - 1) = max(b(il, i - 1), 0.)
234                     endif                     endif
235    
236                     ! limit magnitude of mp(i) to meet cfl condition                     ! limit magnitude of mp(i) to meet cfl condition
   
237                     ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti                     ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
238                     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti                     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
239                     ampmax = amin1(ampmax, amp2)                     ampmax = min(ampmax, amp2)
240                     mp(il, i) = amin1(mp(il, i), ampmax)                     mp(il, i) = min(mp(il, i), ampmax)
241    
242                     ! force mp to decrease linearly to zero                     ! force mp to decrease linearly to zero
243                     ! between cloud base and the surface                     ! between cloud base and the surface
244                       if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
245                     if (p(il, i) > p(il, icb(il))) then                          * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
                       mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) &  
                            / (p(il, 1) - p(il, icb(il)))  
                    endif  
246                  endif ! i == 1                  endif ! i == 1
247    
248                  ! find mixing ratio of precipitating downdraft                  ! find mixing ratio of precipitating downdraft
249    
250                  if (i /= inb(il)) then                  if (i /= inb(il)) then
251                     rp(il, i) = rr(il, i)                     qp(il, i) = q(il, i)
252    
253                     if (mp(il, i) > mp(il, i + 1)) then                     if (mp(il, i) > mp(il, i + 1)) then
254                        rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &                        qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
255                             * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &                             * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
256                             * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &                             * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
257                             * (evap(il, i + 1) + evap(il, i))                             * (evap(il, i + 1) + evap(il, i))
258                        rp(il, i) = rp(il, i) / mp(il, i)                        qp(il, i) = qp(il, i) / mp(il, i)
259                        up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &                        up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
260                             * (mp(il, i) - mp(il, i + 1))                             * (mp(il, i) - mp(il, i + 1))
261                        up(il, i) = up(il, i) / mp(il, i)                        up(il, i) = up(il, i) / mp(il, i)
# Line 257  contains Line 264  contains
264                        vp(il, i) = vp(il, i) / mp(il, i)                        vp(il, i) = vp(il, i) / mp(il, i)
265                     else                     else
266                        if (mp(il, i + 1) > 1e-16) then                        if (mp(il, i + 1) > 1e-16) then
267                           rp(il, i) = rp(il, i + 1) &                           qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
268                                + 100. * ginv * 0.5 * sigd * (ph(il, i) &                                * (ph(il, i) - ph(il, i + 1)) &
269                                - ph(il, i + 1)) &                                * (evap(il, i + 1) + evap(il, i)) / mp(il, i + 1)
                               * (evap(il, i + 1) + evap(il, i)) &  
                               / mp(il, i + 1)  
270                           up(il, i) = up(il, i + 1)                           up(il, i) = up(il, i + 1)
271                           vp(il, i) = vp(il, i + 1)                           vp(il, i) = vp(il, i + 1)
272                        endif                        endif
273                     endif                     endif
274                     rp(il, i) = amin1(rp(il, i), rs(il, i))  
275                     rp(il, i) = amax1(rp(il, i), 0.)                     qp(il, i) = min(qp(il, i), qs(il, i))
276                       qp(il, i) = max(qp(il, i), 0.)
277                  endif                  endif
278               endif               endif
279            end do            end do
Legend:
Removed from v.187  
changed lines
  Added in v.192
  ViewVC Help
Powered by ViewVC 1.1.21