/[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 196 by guez, Mon May 23 13:50:39 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, grav      use cv_thermo_m, only: ginv
14        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}
     ! {ph(i, icb(i) + 1) < plcl(i) <= ph(i, icb(i))}  
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(:, :) ! (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(:, :) ! (klon, klev)      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      real, intent(in):: lv(:, :) ! (klon, klev)  
32      real, intent(in):: cpn(:, :) ! (klon, 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)      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        ! (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    
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)
52      real, intent(out):: water(:, :), evap(:, :) ! (ncum, nl)  
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:
# Line 51  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, pr2, 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 64  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 95  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) = grav &         forall (il = 1:ncum, inb(il) >= i .and. lwork(il)) wdtrain(il) = rg &
115              * (ep(il, i) * m(il, i) * clw(il, i) &              * (ep(il, i) * m(il, i) * clw(il, i) &
116              + sum(max(elij(il, :i - 1, i) - (1. - ep(il, i)) * clw(il, i), 0.) &              + sum(max(elij(il, :i - 1, i) - (1. - ep(il, i)) * clw(il, i), 0.) &
117              * ment(il, :i - 1, i)))              * ment(il, :i - 1, i)))
# Line 108  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 121  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 137  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 = max(0., min(1., &                  pr1 = max(0., min(1., &
164                    (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))                       (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))
165               pr2 = max(0., min(1., &                  sigt = sigp * pr1 + 1. - pr1
166                    (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))))               else
167               sigt = sigp * 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 &
# Line 176  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 204  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
Legend:
Removed from v.196  
changed lines
  Added in v.201
  ViewVC Help
Powered by ViewVC 1.1.21