/[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 186 by guez, Mon Mar 21 15:36:26 2016 UTC revision 188 by guez, Tue Mar 22 16:31:39 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, rr, rs, gz, u, v, p, ph, th, tv, lv, &
8         ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, &         cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, up, vp, wt, &
9         up, vp, wt, water, evap, b)         water, evap, b)
10    
11      use cv30_param_m, only: nl, sigd      use cv30_param_m, only: nl, sigd
     use cvflag, only: cvflag_grav  
12      use cvthermo, only: cpd, ginv, grav      use cvthermo, only: cpd, ginv, grav
13        USE dimphy, ONLY: klon, klev
14    
15      ! inputs:      ! inputs:
     integer, intent(in):: nloc, ncum, nd, na  
16      integer, intent(in):: icb(:), inb(:) ! (ncum)      integer, intent(in):: icb(:), inb(:) ! (ncum)
17      real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)      real t(klon, klev), rr(klon, klev), rs(klon, klev)
18      real gz(nloc, na)      real gz(klon, klev)
19      real u(nloc, nd), v(nloc, nd)      real u(klon, klev), v(klon, klev)
20      real p(nloc, nd), ph(nloc, nd + 1)      real p(klon, klev), ph(klon, klev + 1)
21      real th(nloc, na)      real th(klon, klev)
22      real tv(nloc, na)      real tv(klon, klev)
23      real lv(nloc, na)      real lv(klon, klev)
24      real cpn(nloc, na)      real cpn(klon, klev)
25      real ep(nloc, na), sigp(nloc, na), clw(nloc, na)      real, intent(in):: ep(klon, klev), sigp(klon, klev), clw(klon, klev)
26      real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)      real m(klon, klev), ment(klon, klev, klev), elij(klon, klev, klev)
27      real, intent(in):: delt      real, intent(in):: delt
28      real plcl(nloc)      real plcl(klon)
29    
30      ! outputs:      ! outputs:
31      real mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)      real mp(klon, klev), rp(klon, klev), up(klon, klev), vp(klon, klev)
32      real wt(nloc, na), water(nloc, na), evap(nloc, na)      real wt(klon, klev), water(klon, klev), evap(klon, klev)
33      real b(:, :) ! (nloc, na)      real, intent(out):: b(:, :) ! (ncum, nl)
34    
35      ! Local:      ! Local:
36        integer ncum
37      integer i, j, il, num1      integer i, j, il, num1
38      real tinv, delti      real tinv, delti
39      real awat, afac, afac1, afac2, bfac      real awat, afac, afac1, afac2, bfac
40      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real pr1, pr2, sigt, b6, c6, revap, tevap, delth
41      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
42      real ampmax      real ampmax
43      real lvcp(nloc, na)      real lvcp(klon, klev)
44      real wdtrain(nloc)      real wdtrain(size(icb))
45      logical lwork(nloc)      logical lwork(size(icb))
46    
47      !------------------------------------------------------      !------------------------------------------------------
48    
49        ncum = size(icb)
50      delti = 1. / delt      delti = 1. / delt
51      tinv = 1. / 3.      tinv = 1. / 3.
52      mp = 0.      mp = 0.
# Line 66  contains Line 67  contains
67    
68      ! check whether ep(inb) = 0, if so, skip precipitating      ! check whether ep(inb) = 0, if so, skip precipitating
69      ! downdraft calculation      ! downdraft calculation
70        forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
71    
72      do il = 1, ncum      wdtrain = 0.
        lwork(il) = .TRUE.  
        if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE.  
     enddo  
   
     wdtrain(:ncum) = 0.  
73    
74      downdraft_loop: DO i = nl + 1, 1, - 1      downdraft_loop: DO i = nl - 1, 1, - 1
75         num1 = 0         num1 = 0
76    
77         do il = 1, ncum         do il = 1, ncum
# Line 89  contains Line 86  contains
86    
87            do il = 1, ncum            do il = 1, ncum
88               if (i <= inb(il) .and. lwork(il)) then               if (i <= inb(il) .and. lwork(il)) then
89                  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  
90               endif               endif
91            enddo            enddo
92    
# Line 102  contains Line 95  contains
95                  do il = 1, ncum                  do il = 1, ncum
96                     if (i <= inb(il) .and. lwork(il)) then                     if (i <= inb(il) .and. lwork(il)) then
97                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
98                        awat = amax1(awat, 0.)                        awat = max(awat, 0.)
99                        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  
100                     endif                     endif
101                  enddo                  enddo
102               end do               end do
103            endif            endif
104    
105            ! find rain water and evaporation using provisional            ! find rain water and evaporation using provisional
106            ! estimates of rp(i)and rp(i - 1)            ! estimates of rp(i) and rp(i - 1)
107    
108            do il = 1, ncum            do il = 1, ncum
109               if (i <= inb(il) .and. lwork(il)) then               if (i <= inb(il) .and. lwork(il)) then
# Line 126  contains Line 114  contains
114                          - 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)
115                     rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))                     rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))
116                  endif                  endif
117                  rp(il, i) = amax1(rp(il, i), 0.)  
118                  rp(il, i) = amin1(rp(il, i), rs(il, i))                  rp(il, i) = max(rp(il, i), 0.)
119                    rp(il, i) = min(rp(il, i), rs(il, i))
120                  rp(il, inb(il)) = rr(il, inb(il))                  rp(il, inb(il)) = rr(il, inb(il))
121    
122                  if (i == 1) then                  if (i == 1) then
# Line 137  contains Line 126  contains
126                     rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &                     rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &
127                          - 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)
128                     rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))                     rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))
129                     rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))                     rp(il, i - 1) = min(rp(il, i - 1), rs(il, i - 1))
130                     rp(il, i - 1) = amax1(rp(il, i - 1), 0.)                     rp(il, i - 1) = max(rp(il, i - 1), 0.)
131                     afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &                     afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &
132                          / (1e4 + 2000. * p(il, i) * rs(il, i))                          / (1e4 + 2000. * p(il, i) * rs(il, i))
133                     afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &                     afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &
134                          / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))                          / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))
135                     afac = 0.5 * (afac1 + afac2)                     afac = 0.5 * (afac1 + afac2)
136                  endif                  endif
137                  if (i == inb(il))afac = 0.  
138                  afac = amax1(afac, 0.)                  if (i == inb(il)) afac = 0.
139                    afac = max(afac, 0.)
140                  bfac = 1. / (sigd * wt(il, i))                  bfac = 1. / (sigd * wt(il, i))
141    
142                  ! prise en compte de la variation progressive de sigt dans                  ! prise en compte de la variation progressive de sigt dans
# Line 180  contains Line 170  contains
170                  ! hydrostatic approximation                  ! hydrostatic approximation
171    
172                  if (i /= 1) then                  if (i /= 1) then
173                     tevap = amax1(0., evap(il, i))                     tevap = max(0., evap(il, i))
174                     delth = amax1(0.001, (th(il, i) - th(il, i - 1)))                     delth = max(0.001, (th(il, i) - th(il, i - 1)))
175                     if (cvflag_grav) then                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
176                        mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &                          * (p(il, i - 1) - p(il, i)) / delth
177                             * (p(il, i - 1) - p(il, i)) / delth  
178                     else                     ! If hydrostatic assumption fails, solve cubic
179                        mp(il, i) = 10. * lvcp(il, i) * sigd * tevap &                     ! difference equation for downdraft theta and mass
180                             * (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  
181    
182                     amfac = sigd * sigd * 70. * ph(il, i) &                     amfac = sigd * sigd * 70. * ph(il, i) &
183                          * (p(il, i - 1) - p(il, i)) &                          * (p(il, i - 1) - p(il, i)) &
184                          * (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))
185                     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) &
186                          * mp(il, i))                          * mp(il, i))
187                     if (amp2 > (0.1 * amfac)) then  
188                       if (amp2 > 0.1 * amfac) then
189                        xf = 100. * sigd * sigd * sigd * (ph(il, i) &                        xf = 100. * sigd * sigd * sigd * (ph(il, i) &
190                             - ph(il, i + 1))                             - ph(il, i + 1))
191                        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 195  contains
195                             * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &                             * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
196                             - p(il, i)) * xf * tevap                             - p(il, i)) * xf * tevap
197                        fac2 = 1.                        fac2 = 1.
198                        if (bf < 0.)fac2 = - 1.                        if (bf < 0.) fac2 = - 1.
199                        bf = abs(bf)                        bf = abs(bf)
200                        ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv                        ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
201    
202                        if (ur >= 0.) then                        if (ur >= 0.) then
203                           sru = sqrt(ur)                           sru = sqrt(ur)
204                           fac = 1.                           fac = 1.
205                           if ((0.5 * bf - sru) < 0.)fac = - 1.                           if ((0.5 * bf - sru) < 0.) fac = - 1.
206                           mp(il, i) = mp(il, i + 1) * tinv &                           mp(il, i) = mp(il, i + 1) * tinv &
207                                + (0.5 * bf + sru)**tinv &                                + (0.5 * bf + sru)**tinv &
208                                + fac * (abs(0.5 * bf - sru))**tinv                                + fac * (abs(0.5 * bf - sru))**tinv
# Line 225  contains Line 212  contains
212                           mp(il, i) = mp(il, i + 1) * tinv + 2. &                           mp(il, i) = mp(il, i + 1) * tinv + 2. &
213                                * sqrt(af * tinv) * cos(d * tinv)                                * sqrt(af * tinv) * cos(d * tinv)
214                        endif                        endif
                       mp(il, i) = amax1(0., mp(il, i))  
215    
216                        if (cvflag_grav) then                        mp(il, i) = max(0., mp(il, i))
217                           ! Il y a vraisemblablement une erreur dans la  
218                           ! ligne 2 suivante: il faut diviser par (mp(il,                        ! Il y a vraisemblablement une erreur dans la
219                           ! i) * sigd * grav) et non par (mp(il, i) + sigd                        ! ligne suivante : il faut diviser par (mp(il,
220                           ! * 0.1).  Et il faut bien revoir les facteurs                        ! i) * sigd * grav) et non par (mp(il, i) + sigd
221                           ! 100.                        ! * 0.1).  Et il faut bien revoir les facteurs
222                           b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &                        ! 100.
223                                - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &                        b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
224                                - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
225                                / (lvcp(il, i) * sigd * th(il, i))                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
226                        else                             / (lvcp(il, i) * sigd * th(il, i))
227                           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.)  
228                     endif                     endif
229    
230                     ! limit magnitude of mp(i) to meet cfl condition                     ! limit magnitude of mp(i) to meet cfl condition
   
231                     ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti                     ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
232                     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti                     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
233                     ampmax = amin1(ampmax, amp2)                     ampmax = min(ampmax, amp2)
234                     mp(il, i) = amin1(mp(il, i), ampmax)                     mp(il, i) = min(mp(il, i), ampmax)
235    
236                     ! force mp to decrease linearly to zero                     ! force mp to decrease linearly to zero
237                     ! between cloud base and the surface                     ! between cloud base and the surface
238                       if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
239                     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  
240                  endif ! i == 1                  endif ! i == 1
241    
242                  ! find mixing ratio of precipitating downdraft                  ! find mixing ratio of precipitating downdraft
# Line 268  contains Line 245  contains
245                     rp(il, i) = rr(il, i)                     rp(il, i) = rr(il, i)
246    
247                     if (mp(il, i) > mp(il, i + 1)) then                     if (mp(il, i) > mp(il, i + 1)) then
248                        if (cvflag_grav) then                        rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
249                           rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &                             * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
250                                * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &                             * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
251                                * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &                             * (evap(il, i + 1) + evap(il, i))
                               * (evap(il, i + 1) + evap(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  
252                        rp(il, i) = rp(il, i) / mp(il, i)                        rp(il, i) = rp(il, i) / mp(il, i)
253                        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) &
254                             * (mp(il, i) - mp(il, i + 1))                             * (mp(il, i) - mp(il, i + 1))
# Line 288  contains Line 258  contains
258                        vp(il, i) = vp(il, i) / mp(il, i)                        vp(il, i) = vp(il, i) / mp(il, i)
259                     else                     else
260                        if (mp(il, i + 1) > 1e-16) then                        if (mp(il, i + 1) > 1e-16) then
261                           if (cvflag_grav) then                           rp(il, i) = rp(il, i + 1) + 100. * ginv * 0.5 * sigd &
262                              rp(il, i) = rp(il, i + 1) &                                * (ph(il, i) - ph(il, i + 1)) &
263                                   + 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  
264                           up(il, i) = up(il, i + 1)                           up(il, i) = up(il, i + 1)
265                           vp(il, i) = vp(il, i + 1)                           vp(il, i) = vp(il, i + 1)
266                        endif                        endif
267                     endif                     endif
268                     rp(il, i) = amin1(rp(il, i), rs(il, i))  
269                     rp(il, i) = amax1(rp(il, i), 0.)                     rp(il, i) = min(rp(il, i), rs(il, i))
270                       rp(il, i) = max(rp(il, i), 0.)
271                  endif                  endif
272               endif               endif
273            end do            end do
Legend:
Removed from v.186  
changed lines
  Added in v.188
  ViewVC Help
Powered by ViewVC 1.1.21