/[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 200 by guez, Tue May 31 16:22:42 2016 UTC revision 201 by guez, Mon Jun 6 17:42:15 2016 UTC
# Line 10  contains Line 10  contains
10      ! Unsaturated (precipitating) downdrafts      ! 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      use cv_thermo_m, only: ginv
14      use SUPHEC_M, only: rg      use SUPHEC_M, only: rg, rcpd
15    
16      integer, intent(in):: icb(:) ! (ncum)      integer, intent(in):: icb(:) ! (ncum)
17      ! {2 <= icb <= nl - 3}      ! {2 <= icb <= nl - 3}
# Line 20  contains Line 20  contains
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(:, :) ! (ncum, nl)      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(:, :) ! (ncum, nl)      real, intent(in):: u(:, :), v(:, :) ! (ncum, nl)
27      real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa      real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa
28      real, intent(in):: ph(:, :) ! (ncum, klev + 1)      real, intent(in):: ph(:, :) ! (ncum, klev + 1)
29      real, intent(in):: th(:, :) ! (ncum, nl - 1)      real, intent(in):: th(:, :) ! (ncum, nl - 1) potential temperature, in K
30      real, intent(in):: tv(:, :) ! (klon, klev)      real, intent(in):: tv(:, :) ! (klon, klev)
31    
32      real, intent(in):: lv(:, :) ! (ncum, nl)      real, intent(in):: lv(:, :) ! (ncum, nl)
33      real, intent(in):: cpn(:, :) ! (klon, klev)      ! 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)      real, intent(in):: ep(:, :) ! (ncum, klev)
39      real, intent(in):: clw(:, :) ! (ncum, klev)      real, intent(in):: clw(:, :) ! (ncum, klev)
40      real, intent(in):: m(:, :) ! (ncum, klev)      real, intent(in):: m(:, :) ! (ncum, klev)
# Line 37  contains Line 43  contains
43      real, intent(in):: delt      real, intent(in):: delt
44      real, intent(in):: plcl(:) ! (ncum)      real, intent(in):: plcl(:) ! (ncum)
45    
46      real, intent(out):: mp(:, :) ! (klon, klev)      real, intent(out):: mp(:, :)
47      ! mass flux of the unsaturated downdraft, defined positive downward      ! (ncum, nl) Mass flux of the unsaturated downdraft, defined
48      ! M_p in Emanuel (1991 928)      ! positive downward, in kg m-2 s-1.  M_p in Emanuel (1991 928).
49    
50      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
51      real, intent(out):: wt(:, :) ! (ncum, nl)      real, intent(out):: wt(:, :) ! (ncum, nl)
# Line 61  contains Line 67  contains
67    
68      integer ncum      integer ncum
69      integer i, il, imax      integer i, il, imax
70      real tinv, delti      real, parameter:: tinv = 1. / 3.
71        real  delti
72      real afac, afac1, afac2, bfac      real afac, afac1, afac2, bfac
73      real pr1, sigt, b6, c6, revap, tevap, delth      real pr1, sigt, b6, c6, revap, tevap, delth
74      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real xf, tf, fac2, ur, sru, fac, d, af, bf
75      real ampmax      real ampmax
76      real lvcp(size(icb), nl) ! (ncum, nl)      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 74  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 118  contains Line 124  contains
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 131  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 188  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**2 * 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) &                       * 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) &
199                       * (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))) &
200                  amp2 = abs(mp(il, i + 1)**2 - mp(il, i)**2)                       then
   
                 if (amp2 > 0.1 * amfac) then  
201                     xf = 100. * sigd**3 * (ph(il, i) - ph(il, i + 1))                     xf = 100. * sigd**3 * (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))
# Line 214  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
# Line 225  contains Line 229  contains
229                     ! suivante : il faut diviser par (mp(il, i) * sigd                     ! suivante : il faut diviser par (mp(il, i) * sigd
230                     ! * rg) et non par (mp(il, i) + sigd * 0.1).  Et il                     ! * rg) et non par (mp(il, i) + sigd * 0.1).  Et il
231                     ! faut bien revoir les facteurs 100.                     ! faut bien revoir les facteurs 100.
232                     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &                     b(il, i - 1) = max(b(il, i) + 100. * (p(il, i - 1) &
233                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &                          - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) - 10. &
234                          - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &                          * (th(il, i) - th(il, i - 1)) * t(il, i) &
235                          * th(il, i))                          / (lvcp(il, i) * sigd * th(il, i)), 0.)
                    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
Legend:
Removed from v.200  
changed lines
  Added in v.201
  ViewVC Help
Powered by ViewVC 1.1.21