/[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 187 by guez, Mon Mar 21 18:01:02 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(ncum, icb, inb, t, rr, rs, gz, u, v, p, ph, th, tv, &
8         ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, &         lv, 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:
16      integer, intent(in):: nloc, ncum, nd, na      integer, intent(in):: ncum
17      integer, intent(in):: icb(:), inb(:) ! (ncum)      integer, intent(in):: icb(:), inb(:) ! (ncum)
18      real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)      real t(klon, klev), rr(klon, klev), rs(klon, klev)
19      real gz(nloc, na)      real gz(klon, klev)
20      real u(nloc, nd), v(nloc, nd)      real u(klon, klev), v(klon, klev)
21      real p(nloc, nd), ph(nloc, nd + 1)      real p(klon, klev), ph(klon, klev + 1)
22      real th(nloc, na)      real th(klon, klev)
23      real tv(nloc, na)      real tv(klon, klev)
24      real lv(nloc, na)      real lv(klon, klev)
25      real cpn(nloc, na)      real cpn(klon, klev)
26      real ep(nloc, na), sigp(nloc, na), clw(nloc, na)      real, intent(in):: ep(klon, klev), sigp(klon, klev), clw(klon, klev)
27      real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)      real m(klon, klev), ment(klon, klev, klev), elij(klon, klev, klev)
28      real, intent(in):: delt      real, intent(in):: delt
29      real plcl(nloc)      real plcl(klon)
30    
31      ! outputs:      ! outputs:
32      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)
33      real wt(nloc, na), water(nloc, na), evap(nloc, na)      real wt(klon, klev), water(klon, klev), evap(klon, klev)
34      real b(:, :) ! (nloc, na)      real b(:, :) ! (ncum, klev)
35    
36      ! Local:      ! Local:
37      integer i, j, il, num1      integer i, j, il, num1
# Line 40  contains Line 40  contains
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(ncum)
45      logical lwork(nloc)      logical lwork(ncum)
46    
47      !------------------------------------------------------      !------------------------------------------------------
48    
# Line 66  contains Line 66  contains
66    
67      ! check whether ep(inb) = 0, if so, skip precipitating      ! check whether ep(inb) = 0, if so, skip precipitating
68      ! downdraft calculation      ! downdraft calculation
69        forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
70    
71      do il = 1, ncum      wdtrain = 0.
        lwork(il) = .TRUE.  
        if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE.  
     enddo  
   
     wdtrain(:ncum) = 0.  
72    
73      downdraft_loop: DO i = nl + 1, 1, - 1      downdraft_loop: DO i = nl - 1, 1, - 1
74         num1 = 0         num1 = 0
75    
76         do il = 1, ncum         do il = 1, ncum
# Line 89  contains Line 85  contains
85    
86            do il = 1, ncum            do il = 1, ncum
87               if (i <= inb(il) .and. lwork(il)) then               if (i <= inb(il) .and. lwork(il)) then
88                  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  
89               endif               endif
90            enddo            enddo
91    
# Line 103  contains Line 95  contains
95                     if (i <= inb(il) .and. lwork(il)) then                     if (i <= inb(il) .and. lwork(il)) then
96                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
97                        awat = amax1(awat, 0.)                        awat = amax1(awat, 0.)
98                        if (cvflag_grav) then                        wdtrain(il) = wdtrain(il) + grav * awat &
99                           wdtrain(il) = wdtrain(il) + grav * awat &                             * ment(il, j, i)
                               * 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
# Line 182  contains Line 170  contains
170                  if (i /= 1) then                  if (i /= 1) then
171                     tevap = amax1(0., evap(il, i))                     tevap = amax1(0., evap(il, i))
172                     delth = amax1(0.001, (th(il, i) - th(il, i - 1)))                     delth = amax1(0.001, (th(il, i) - th(il, i - 1)))
173                     if (cvflag_grav) then                     mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
174                        mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &                          * (p(il, i - 1) - p(il, i)) / delth
                            * (p(il, i - 1) - p(il, i)) / delth  
                    else  
                       mp(il, i) = 10. * lvcp(il, i) * sigd * tevap &  
                            * (p(il, i - 1) - p(il, i)) / delth  
                    endif  
175    
176                     ! if hydrostatic assumption fails,                     ! if hydrostatic assumption fails,
177                     ! solve cubic difference equation for downdraft theta                     ! solve cubic difference equation for downdraft theta
# Line 227  contains Line 210  contains
210                        endif                        endif
211                        mp(il, i) = amax1(0., mp(il, i))                        mp(il, i) = amax1(0., mp(il, i))
212    
213                        if (cvflag_grav) then                        ! Il y a vraisemblablement une erreur dans la
214                           ! Il y a vraisemblablement une erreur dans la                        ! ligne 2 suivante: il faut diviser par (mp(il,
215                           ! ligne 2 suivante: il faut diviser par (mp(il,                        ! i) * sigd * grav) et non par (mp(il, i) + sigd
216                           ! i) * sigd * grav) et non par (mp(il, i) + sigd                        ! * 0.1).  Et il faut bien revoir les facteurs
217                           ! * 0.1).  Et il faut bien revoir les facteurs                        ! 100.
218                           ! 100.                        b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
219                           b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &                             - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
220                                - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &                             - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
221                                - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &                             / (lvcp(il, i) * sigd * th(il, i))
                               / (lvcp(il, i) * sigd * th(il, i))  
                       else  
                          b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &  
                               - 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  
222                        b(il, i - 1) = amax1(b(il, i - 1), 0.)                        b(il, i - 1) = amax1(b(il, i - 1), 0.)
223                     endif                     endif
224    
# Line 268  contains Line 244  contains
244                     rp(il, i) = rr(il, i)                     rp(il, i) = rr(il, i)
245    
246                     if (mp(il, i) > mp(il, i + 1)) then                     if (mp(il, i) > mp(il, i + 1)) then
247                        if (cvflag_grav) then                        rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
248                           rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &                             * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
249                                * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &                             * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
250                                * 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  
251                        rp(il, i) = rp(il, i) / mp(il, i)                        rp(il, i) = rp(il, i) / mp(il, i)
252                        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) &
253                             * (mp(il, i) - mp(il, i + 1))                             * (mp(il, i) - mp(il, i + 1))
# Line 288  contains Line 257  contains
257                        vp(il, i) = vp(il, i) / mp(il, i)                        vp(il, i) = vp(il, i) / mp(il, i)
258                     else                     else
259                        if (mp(il, i + 1) > 1e-16) then                        if (mp(il, i + 1) > 1e-16) then
260                           if (cvflag_grav) then                           rp(il, i) = rp(il, i + 1) &
261                              rp(il, i) = rp(il, i + 1) &                                + 100. * ginv * 0.5 * sigd * (ph(il, i) &
262                                   + 100. * ginv * 0.5 * sigd * (ph(il, i) &                                - ph(il, i + 1)) &
263                                   - ph(il, i + 1)) &                                * (evap(il, i + 1) + evap(il, i)) &
264                                   * (evap(il, i + 1) + evap(il, i)) &                                / mp(il, i + 1)
                                  / 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  
265                           up(il, i) = up(il, i + 1)                           up(il, i) = up(il, i + 1)
266                           vp(il, i) = vp(il, i + 1)                           vp(il, i) = vp(il, i + 1)
267                        endif                        endif
Legend:
Removed from v.186  
changed lines
  Added in v.187
  ViewVC Help
Powered by ViewVC 1.1.21