/[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 193 by guez, Thu May 12 13:22:19 2016 UTC revision 201 by guez, Mon Jun 6 17:42:15 2016 UTC
# Line 5  module cv30_unsat_m Line 5  module cv30_unsat_m
5  contains  contains
6    
7    SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &    SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &
8         ep, sigp, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, &         ep, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, evap, b)
9         evap, b)  
10        ! Unsaturated (precipitating) downdrafts
11    
12      use cv30_param_m, only: nl, sigd      use cv30_param_m, only: nl, sigd
13      use cv_thermo_m, only: cpd, ginv, grav      use cv_thermo_m, only: ginv
14      USE dimphy, ONLY: klon, klev      use SUPHEC_M, only: rg, rcpd
15    
16      integer, intent(in):: icb(:) ! (ncum)      integer, intent(in):: icb(:) ! (ncum)
17        ! {2 <= icb <= nl - 3}
18    
19      integer, intent(in):: inb(:) ! (ncum)      integer, intent(in):: inb(:) ! (ncum)
20      ! first model level above the level of neutral buoyancy of the      ! first model level above the level of neutral buoyancy of the
21      ! parcel (1 <= inb <= nl - 1)      ! parcel (1 <= inb <= nl - 1)
22    
23      real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (klon, klev)      real, intent(in):: t(:, :) ! (ncum, nl) temperature (K)
24        real, intent(in):: q(:, :), qs(:, :) ! (ncum, nl)
25      real, intent(in):: gz(:, :) ! (klon, klev)      real, intent(in):: gz(:, :) ! (klon, klev)
26      real, intent(in):: u(:, :), v(:, :) ! (klon, klev)      real, intent(in):: u(:, :), v(:, :) ! (ncum, nl)
27      real, intent(in):: p(klon, klev), ph(klon, klev + 1)      real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa
28      real, intent(in):: th(klon, klev)      real, intent(in):: ph(:, :) ! (ncum, klev + 1)
29      real, intent(in):: tv(klon, klev)      real, intent(in):: th(:, :) ! (ncum, nl - 1) potential temperature, in K
30      real, intent(in):: lv(klon, klev)      real, intent(in):: tv(:, :) ! (klon, klev)
31      real, intent(in):: cpn(klon, klev)  
32      real, intent(in):: ep(:, :), sigp(:, :), clw(:, :) ! (ncum, klev)      real, intent(in):: lv(:, :) ! (ncum, nl)
33        ! specific latent heat of vaporization of water, in J kg-1
34    
35        real, intent(in):: cpn(:, :) ! (ncum, nl)
36        ! specific heat capacity at constant pressure of humid air, in J K-1 kg-1
37    
38        real, intent(in):: ep(:, :) ! (ncum, klev)
39        real, intent(in):: clw(:, :) ! (ncum, klev)
40      real, intent(in):: m(:, :) ! (ncum, klev)      real, intent(in):: m(:, :) ! (ncum, klev)
41      real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)      real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)
42      real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)      real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)
43      real, intent(in):: delt      real, intent(in):: delt
44      real, intent(in):: plcl(klon)      real, intent(in):: plcl(:) ! (ncum)
45    
46        real, intent(out):: mp(:, :)
47        ! (ncum, nl) Mass flux of the unsaturated downdraft, defined
48        ! positive downward, in kg m-2 s-1.  M_p in Emanuel (1991 928).
49    
     ! outputs:  
     real, intent(out):: mp(klon, klev)  
50      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
51      real wt(klon, klev), water(klon, klev), evap(klon, klev)      real, intent(out):: wt(:, :) ! (ncum, nl)
52    
53        real, intent(out):: water(:, :) ! (ncum, nl)
54        ! precipitation mixing ratio, l_p in Emanuel (1991 928)
55    
56        real, intent(out):: evap(:, :) ! (ncum, nl)
57        ! sigt * rate of evaporation of precipitation, in s-1
58        ! \sigma_s E in Emanuel (1991 928)
59    
60      real, intent(out):: b(:, :) ! (ncum, nl - 1)      real, intent(out):: b(:, :) ! (ncum, nl - 1)
61    
62      ! Local:      ! Local:
63    
64        real, parameter:: sigp = 0.15
65        ! fraction of precipitation falling outside of cloud, \sig_s in
66        ! Emanuel (1991 928)
67    
68      integer ncum      integer ncum
69      integer i, j, il, imax      integer i, il, imax
70      real tinv, delti      real, parameter:: tinv = 1. / 3.
71      real awat, afac, afac1, afac2, bfac      real  delti
72      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real afac, afac1, afac2, bfac
73      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real pr1, sigt, b6, c6, revap, tevap, delth
74        real xf, tf, fac2, ur, sru, fac, d, af, bf
75      real ampmax      real ampmax
76      real lvcp(klon, klev)      real lvcp(size(icb), nl) ! (ncum, nl) L_v / C_p, in K
77      real wdtrain(size(icb)) ! (ncum)      real wdtrain(size(icb)) ! (ncum)
78      logical lwork(size(icb)) ! (ncum)      logical lwork(size(icb)) ! (ncum)
79    
# Line 55  contains Line 81  contains
81    
82      ncum = size(icb)      ncum = size(icb)
83      delti = 1. / delt      delti = 1. / delt
     tinv = 1. / 3.  
84      mp = 0.      mp = 0.
85      b = 0.      b = 0.
86    
# Line 86  contains Line 111  contains
111         ! and condensed water flux         ! and condensed water flux
112    
113         ! Calculate detrained precipitation         ! Calculate detrained precipitation
114           forall (il = 1:ncum, inb(il) >= i .and. lwork(il)) wdtrain(il) = rg &
115         do il = 1, ncum              * (ep(il, i) * m(il, i) * clw(il, i) &
116            if (i <= inb(il) .and. lwork(il)) then              + sum(max(elij(il, :i - 1, i) - (1. - ep(il, i)) * clw(il, i), 0.) &
117               wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)              * ment(il, :i - 1, i)))
           endif  
        enddo  
   
        if (i > 1) then  
           do j = 1, i - 1  
              do il = 1, ncum  
                 if (i <= inb(il) .and. lwork(il)) then  
                    awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)  
                    awat = max(awat, 0.)  
                    wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)  
                 endif  
              enddo  
           end do  
        endif  
118    
119         ! Find rain water and evaporation using provisional         ! Find rain water and evaporation using provisional
120         ! estimates of qp(i) and qp(i - 1)         ! estimates of qp(i) and qp(i - 1)
121    
122         do il = 1, ncum         loop_horizontal: do il = 1, ncum
123            if (i <= inb(il) .and. lwork(il)) then            if (i <= inb(il) .and. lwork(il)) then
124               wt(il, i) = 45.               wt(il, i) = 45.
125    
126               if (i < inb(il)) then               if (i < inb(il)) then
127                  qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, i)) &                  qp(il, i) = qp(il, i + 1) + (rcpd * (t(il, i + 1) - t(il, i)) &
128                       + gz(il, i + 1) - gz(il, i)) / lv(il, i)                       + gz(il, i + 1) - gz(il, i)) / lv(il, i)
129                  qp(il, i) = 0.5 * (qp(il, i) + q(il, i))                  qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
130               endif               endif
# Line 126  contains Line 137  contains
137                  afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &                  afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
138                       / (1e4 + 2000. * p(il, 1) * qs(il, 1))                       / (1e4 + 2000. * p(il, 1) * qs(il, 1))
139               else               else
140                  qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) - t(il, i - 1)) &                  qp(il, i - 1) = qp(il, i) + (rcpd * (t(il, i) - t(il, i - 1)) &
141                       + gz(il, i) - gz(il, i - 1)) / lv(il, i)                       + gz(il, i) - gz(il, i - 1)) / lv(il, i)
142                  qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))                  qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
143                  qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))                  qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
# Line 142  contains Line 153  contains
153               afac = max(afac, 0.)               afac = max(afac, 0.)
154               bfac = 1. / (sigd * wt(il, i))               bfac = 1. / (sigd * wt(il, i))
155    
156               ! Prise en compte de la variation progressive de sigt dans               if (i <= icb(il)) then
157               ! les couches icb et icb - 1:                  ! Prise en compte de la variation progressive de sigt dans
158               ! pour plcl < ph(i + 1), pr1 = 0 et pr2 = 1                  ! les couches icb et icb - 1 :
159               ! pour plcl > ph(i), pr1 = 1 et pr2 = 0                  ! pour plcl <= ph(i + 1), pr1 = 0
160               ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion \`a cheval                  ! pour plcl >= ph(i), pr1 = 1
161               ! sur le nuage, et pr2 est la proportion sous la base du                  ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion
162               ! nuage.                  ! \`a cheval sur le nuage.
163               pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))                  pr1 = max(0., min(1., &
164               pr1 = max(0., min(1., pr1))                       (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))
165               pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))                  sigt = sigp * pr1 + 1. - pr1
166               pr2 = max(0., min(1., pr2))               else
167               sigt = sigp(il, i) * pr1 + pr2                  ! {i >= icb(il) + 1}
168                    sigt = sigp
169                 end if
170    
171               b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac               b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac
172               c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &               c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
173                    * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)                    * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
174    
175               if (c6 > 0.) then               if (c6 > 0.) then
176                  revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))                  revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
177                  evap(il, i) = sigt * afac * revap                  evap(il, i) = sigt * afac * revap
# Line 180  contains Line 194  contains
194                  ! If hydrostatic assumption fails, solve cubic                  ! If hydrostatic assumption fails, solve cubic
195                  ! difference equation for downdraft theta and mass                  ! difference equation for downdraft theta and mass
196                  ! flux from two simultaneous differential equations                  ! flux from two simultaneous differential equations
197                    if (abs(mp(il, i + 1)**2 - mp(il, i)**2) > 0.1 * sigd**2 &
198                  amfac = sigd * sigd * 70. * ph(il, i) &                       * 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) &
199                       * (p(il, i - 1) - p(il, i)) &                       * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))) &
200                       * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))                       then
201                  amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &                     xf = 100. * sigd**3 * (ph(il, i) - ph(il, i + 1))
                      * mp(il, i))  
   
                 if (amp2 > 0.1 * amfac) then  
                    xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))  
202                     tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &                     tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
203                          * t(il, i) / (lvcp(il, i) * sigd * th(il, i))                          * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
204                     af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv                     af = xf * tf + mp(il, i + 1)**2 * tinv
205                     bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &                     bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
206                          * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &                          * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
207                          - p(il, i)) * xf * tevap                          - p(il, i)) * xf * tevap
# Line 208  contains Line 218  contains
218                             + sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv                             + sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv
219                     else                     else
220                        d = atan(2. * sqrt(- ur) / (bf + 1e-28))                        d = atan(2. * sqrt(- ur) / (bf + 1e-28))
221                        if (fac2 < 0.)d = 3.14159 - d                        if (fac2 < 0.) d = 3.14159 - d
222                        mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &                        mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &
223                             * cos(d * tinv)                             * cos(d * tinv)
224                     endif                     endif
225    
226                     mp(il, i) = max(0., mp(il, i))                     mp(il, i) = max(0., mp(il, i))
227    
228                     ! Il y a vraisemblablement une erreur dans la                     ! Il y a vraisemblablement une erreur dans la ligne
229                     ! ligne suivante : il faut diviser par (mp(il,                     ! suivante : il faut diviser par (mp(il, i) * sigd
230                     ! i) * sigd * grav) et non par (mp(il, i) + sigd                     ! * rg) et non par (mp(il, i) + sigd * 0.1).  Et il
231                     ! * 0.1).  Et il faut bien revoir les facteurs                     ! faut bien revoir les facteurs 100.
232                     ! 100.                     b(il, i - 1) = max(b(il, i) + 100. * (p(il, i - 1) &
233                     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &                          - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) - 10. &
234                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &                          * (th(il, i) - th(il, i - 1)) * t(il, i) &
235                          - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &                          / (lvcp(il, i) * sigd * th(il, i)), 0.)
                         * th(il, i))  
                    b(il, i - 1) = max(b(il, i - 1), 0.)  
236                  endif                  endif
237    
238                  ! Limit magnitude of mp(i) to meet CFL condition:                  ! Limit magnitude of mp to meet CFL condition:
239                  ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti                  ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
240                  amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti                  ampmax = min(ampmax, 2. * (ph(il, i - 1) - ph(il, i)) * delti)
                 ampmax = min(ampmax, amp2)  
241                  mp(il, i) = min(mp(il, i), ampmax)                  mp(il, i) = min(mp(il, i), ampmax)
242    
243                  ! Force mp to decrease linearly to zero between cloud                  ! Force mp to decrease linearly to zero between cloud
# Line 270  contains Line 277  contains
277                  qp(il, i) = max(qp(il, i), 0.)                  qp(il, i) = max(qp(il, i), 0.)
278               endif               endif
279            endif            endif
280         end do         end do loop_horizontal
281      end DO downdraft_loop      end DO downdraft_loop
282    
283    end SUBROUTINE cv30_unsat    end SUBROUTINE cv30_unsat
Legend:
Removed from v.193  
changed lines
  Added in v.201
  ViewVC Help
Powered by ViewVC 1.1.21