/[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 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(ncum, icb, inb, t, rr, rs, gz, u, v, p, ph, th, tv, &    SUBROUTINE cv30_unsat(icb, inb, t, rr, rs, gz, u, v, p, ph, th, tv, lv, &
8         lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, up, vp, wt, &         cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, up, vp, wt, &
9         water, evap, b)         water, evap, b)
10    
11      use cv30_param_m, only: nl, sigd      use cv30_param_m, only: nl, sigd
# Line 13  contains Line 13  contains
13      USE dimphy, ONLY: klon, klev      USE dimphy, ONLY: klon, klev
14    
15      ! inputs:      ! inputs:
     integer, intent(in):: ncum  
16      integer, intent(in):: icb(:), inb(:) ! (ncum)      integer, intent(in):: icb(:), inb(:) ! (ncum)
17      real t(klon, klev), rr(klon, klev), rs(klon, klev)      real t(klon, klev), rr(klon, klev), rs(klon, klev)
18      real gz(klon, klev)      real gz(klon, klev)
# Line 31  contains Line 30  contains
30      ! outputs:      ! outputs:
31      real mp(klon, klev), rp(klon, klev), up(klon, klev), vp(klon, klev)      real mp(klon, klev), rp(klon, klev), up(klon, klev), vp(klon, klev)
32      real wt(klon, klev), water(klon, klev), evap(klon, klev)      real wt(klon, klev), water(klon, klev), evap(klon, klev)
33      real b(:, :) ! (ncum, klev)      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
# Line 41  contains Line 41  contains
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(klon, klev)      real lvcp(klon, klev)
44      real wdtrain(ncum)      real wdtrain(size(icb))
45      logical lwork(ncum)      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 94  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                        wdtrain(il) = wdtrain(il) + grav * awat &                        wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
                            * ment(il, j, i)  
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 114  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 125  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 168  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                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
176                          * (p(il, i - 1) - p(il, i)) / delth                          * (p(il, i - 1) - p(il, i)) / delth
177    
178                     ! if hydrostatic assumption fails,                     ! If hydrostatic assumption fails, solve cubic
179                     ! solve cubic difference equation for downdraft theta                     ! difference equation for downdraft theta and mass
180                     ! and mass flux from two simultaneous differential eqns                     ! flux from two simultaneous differential equations
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 192  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 208  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
215                        mp(il, i) = amax1(0., mp(il, i))  
216                          mp(il, i) = max(0., mp(il, i))
217    
218                        ! Il y a vraisemblablement une erreur dans la                        ! Il y a vraisemblablement une erreur dans la
219                        ! ligne 2 suivante: il faut diviser par (mp(il,                        ! ligne suivante : il faut diviser par (mp(il,
220                        ! i) * sigd * grav) et non par (mp(il, i) + sigd                        ! i) * sigd * grav) et non par (mp(il, i) + sigd
221                        ! * 0.1).  Et il faut bien revoir les facteurs                        ! * 0.1).  Et il faut bien revoir les facteurs
222                        ! 100.                        ! 100.
# Line 219  contains Line 224  contains
224                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
225                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
226                             / (lvcp(il, i) * sigd * th(il, i))                             / (lvcp(il, i) * sigd * th(il, i))
227                        b(il, i - 1) = amax1(b(il, i - 1), 0.)                        b(il, i - 1) = max(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 257  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                           rp(il, i) = rp(il, i + 1) &                           rp(il, i) = rp(il, i + 1) + 100. * ginv * 0.5 * sigd &
262                                + 100. * ginv * 0.5 * sigd * (ph(il, i) &                                * (ph(il, i) - ph(il, i + 1)) &
263                                - 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)  
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.187  
changed lines
  Added in v.188
  ViewVC Help
Powered by ViewVC 1.1.21