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

Diff of /trunk/phylmd/CV30_routines/cv30_unsat.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 186 by guez, Mon Mar 21 15:36:26 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(nloc, ncum, nd, na, icb, inb, t, rr, rs, gz, u, v, p, &    SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &
8         ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, &         ep, sigp, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, &
9         up, vp, wt, water, evap, b)         evap, b)
10    
11      use cv30_param_m, only: nl, sigd      use cv30_param_m, only: nl, sigd
12      use cvflag, only: cvflag_grav      use cv_thermo_m, only: cpd, ginv, grav
13      use cvthermo, only: cpd, ginv, grav      USE dimphy, ONLY: klon, klev
14    
15      ! inputs:      integer, intent(in):: icb(:) ! (ncum)
16      integer, intent(in):: nloc, ncum, nd, na  
17      integer, intent(in):: icb(:), inb(:) ! (ncum)      integer, intent(in):: inb(:) ! (ncum)
18      real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)      ! first model level above the level of neutral buoyancy of the
19      real gz(nloc, na)      ! parcel (1 <= inb <= nl - 1)
20      real u(nloc, nd), v(nloc, nd)  
21      real p(nloc, nd), ph(nloc, nd + 1)      real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (klon, klev)
22      real th(nloc, na)      real, intent(in):: gz(:, :) ! (klon, klev)
23      real tv(nloc, na)      real, intent(in):: u(:, :), v(:, :) ! (klon, klev)
24      real lv(nloc, na)      real, intent(in):: p(klon, klev), ph(klon, klev + 1)
25      real cpn(nloc, na)      real, intent(in):: th(klon, klev)
26      real ep(nloc, na), sigp(nloc, na), clw(nloc, na)      real, intent(in):: tv(klon, klev)
27      real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)      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(nloc)      real, intent(in):: plcl(klon)
35    
36      ! outputs:      ! outputs:
37      real mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)      real, intent(out):: mp(klon, klev)
38      real wt(nloc, na), water(nloc, na), evap(nloc, na)      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
39      real b(:, :) ! (nloc, na)      real wt(klon, klev), water(klon, klev), evap(klon, klev)
40        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
47      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real pr1, pr2, sigt, b6, c6, revap, tevap, delth
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(nloc, na)      real lvcp(klon, klev)
51      real wdtrain(nloc)      real wdtrain(size(icb))
52      logical lwork(nloc)      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
73    
74      ! check whether ep(inb) = 0, if so, skip precipitating      ! check whether ep(inb) = 0, if so, skip precipitating
75      ! downdraft calculation      ! downdraft calculation
76        forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
77    
78      do il = 1, ncum      wdtrain = 0.
        lwork(il) = .TRUE.  
        if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE.  
     enddo  
   
     wdtrain(:ncum) = 0.  
79    
80      downdraft_loop: DO i = nl + 1, 1, - 1      downdraft_loop: DO i = nl - 1, 1, - 1
81         num1 = 0         num1 = 0
82    
83         do il = 1, ncum         do il = 1, ncum
# Line 89  contains Line 92  contains
92    
93            do il = 1, ncum            do il = 1, ncum
94               if (i <= inb(il) .and. lwork(il)) then               if (i <= inb(il) .and. lwork(il)) then
95                  if (cvflag_grav) then                  wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
                    wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)  
                 else  
                    wdtrain(il) = 10. * ep(il, i) * m(il, i) * clw(il, i)  
                 endif  
96               endif               endif
97            enddo            enddo
98    
# Line 102  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                        if (cvflag_grav) then                        wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
                          wdtrain(il) = wdtrain(il) + grav * awat &  
                               * ment(il, j, i)  
                       else  
                          wdtrain(il) = wdtrain(il) + 10. * awat * ment(il, j, i)  
                       endif  
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 180  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                     if (cvflag_grav) then                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
182                        mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &                          * (p(il, i - 1) - p(il, i)) / delth
183                             * (p(il, i - 1) - p(il, i)) / delth  
184                     else                     ! If hydrostatic assumption fails, solve cubic
185                        mp(il, i) = 10. * lvcp(il, i) * sigd * tevap &                     ! difference equation for downdraft theta and mass
186                             * (p(il, i - 1) - p(il, i)) / delth                     ! flux from two simultaneous differential equations
                    endif  
   
                    ! if hydrostatic assumption fails,  
                    ! solve cubic difference equation for downdraft theta  
                    ! and mass flux from two simultaneous differential eqns  
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 209  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 225  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
                       mp(il, i) = amax1(0., mp(il, i))  
221    
222                        if (cvflag_grav) then                        mp(il, i) = max(0., mp(il, i))
223                           ! Il y a vraisemblablement une erreur dans la  
224                           ! ligne 2 suivante: il faut diviser par (mp(il,                        ! Il y a vraisemblablement une erreur dans la
225                           ! i) * sigd * grav) et non par (mp(il, i) + sigd                        ! ligne suivante : il faut diviser par (mp(il,
226                           ! * 0.1).  Et il faut bien revoir les facteurs                        ! i) * sigd * grav) et non par (mp(il, i) + sigd
227                           ! 100.                        ! * 0.1).  Et il faut bien revoir les facteurs
228                           b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &                        ! 100.
229                                - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &                        b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
230                                - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
231                                / (lvcp(il, i) * sigd * th(il, i))                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
232                        else                             / (lvcp(il, i) * sigd * th(il, i))
233                           b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &                        b(il, i - 1) = max(b(il, i - 1), 0.)
                               - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &  
                               - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &  
                               / (lvcp(il, i) * sigd * th(il, i))  
                       endif  
                       b(il, i - 1) = amax1(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                        if (cvflag_grav) then                        qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
255                           rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &                             * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
256                                * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &                             * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
257                                * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &                             * (evap(il, i + 1) + evap(il, i))
258                                * (evap(il, i + 1) + evap(il, i))                        qp(il, i) = qp(il, i) / mp(il, i)
                       else  
                          rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &  
                               * (mp(il, i) - mp(il, i + 1)) + 5. * sigd &  
                               * (ph(il, i) - ph(il, i + 1)) &  
                               * (evap(il, i + 1) + evap(il, i))  
                       endif  
                       rp(il, i) = rp(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 288  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                           if (cvflag_grav) then                           qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
268                              rp(il, i) = rp(il, i + 1) &                                * (ph(il, i) - ph(il, i + 1)) &
269                                   + 100. * ginv * 0.5 * sigd * (ph(il, i) &                                * (evap(il, i + 1) + evap(il, i)) / mp(il, i + 1)
                                  - ph(il, i + 1)) &  
                                  * (evap(il, i + 1) + evap(il, i)) &  
                                  / mp(il, i + 1)  
                          else  
                             rp(il, i) = rp(il, i + 1) &  
                                  + 5. * sigd * (ph(il, i) - ph(il, i + 1)) &  
                                  * (evap(il, i + 1) + evap(il, i)) &  
                                  / mp(il, i + 1)  
                          endif  
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.186  
changed lines
  Added in v.192
  ViewVC Help
Powered by ViewVC 1.1.21