/[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 194 by guez, Thu May 12 14:35:35 2016 UTC revision 199 by guez, Tue May 31 16:22:42 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: cpd, ginv
14      USE dimphy, ONLY: klon, klev      use SUPHEC_M, only: rg
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(:, :), q(:, :), qs(:, :) ! (ncum, nl)
24      real, intent(in):: gz(:, :) ! (klon, klev)      real, intent(in):: gz(:, :) ! (klon, klev)
25      real, intent(in):: u(:, :), v(:, :) ! (klon, klev)      real, intent(in):: u(:, :), v(:, :) ! (ncum, nl)
26      real, intent(in):: p(klon, klev), ph(klon, klev + 1)      real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa
27      real, intent(in):: th(klon, klev)      real, intent(in):: ph(:, :) ! (ncum, klev + 1)
28      real, intent(in):: tv(klon, klev)      real, intent(in):: th(:, :) ! (ncum, nl - 1)
29      real, intent(in):: lv(klon, klev)      real, intent(in):: tv(:, :) ! (klon, klev)
30      real, intent(in):: cpn(klon, klev)      real, intent(in):: lv(:, :) ! (ncum, nl)
31      real, intent(in):: ep(:, :), sigp(:, :), clw(:, :) ! (ncum, klev)      real, intent(in):: cpn(:, :) ! (klon, klev)
32        real, intent(in):: ep(:, :) ! (ncum, klev)
33        real, intent(in):: clw(:, :) ! (ncum, klev)
34      real, intent(in):: m(:, :) ! (ncum, klev)      real, intent(in):: m(:, :) ! (ncum, klev)
35      real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)      real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)
36      real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)      real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)
37      real, intent(in):: delt      real, intent(in):: delt
38      real, intent(in):: plcl(klon)      real, intent(in):: plcl(:) ! (ncum)
39    
40        real, intent(out):: mp(:, :) ! (klon, klev)
41        ! mass flux of the unsaturated downdraft, defined positive downward
42        ! M_p in Emanuel (1991 928)
43    
     ! outputs:  
     real, intent(out):: mp(klon, klev)  
44      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)      real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
45      real wt(klon, klev), water(klon, klev), evap(klon, klev)      real, intent(out):: wt(:, :) ! (ncum, nl)
46    
47        real, intent(out):: water(:, :) ! (ncum, nl)
48        ! precipitation mixing ratio, l_p in Emanuel (1991 928)
49    
50        real, intent(out):: evap(:, :) ! (ncum, nl)
51        ! sigt * rate of evaporation of precipitation, in s-1
52        ! \sigma_s E in Emanuel (1991 928)
53    
54      real, intent(out):: b(:, :) ! (ncum, nl - 1)      real, intent(out):: b(:, :) ! (ncum, nl - 1)
55    
56      ! Local:      ! Local:
57    
58        real, parameter:: sigp = 0.15
59        ! fraction of precipitation falling outside of cloud, \sig_s in
60        ! Emanuel (1991 928)
61    
62      integer ncum      integer ncum
63      integer i, il, imax      integer i, il, imax
64      real tinv, delti      real tinv, delti
65      real afac, afac1, afac2, bfac      real afac, afac1, afac2, bfac
66      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real pr1, sigt, b6, c6, revap, tevap, delth
67      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
68      real ampmax      real ampmax
69      real lvcp(klon, klev)      real lvcp(size(icb), nl) ! (ncum, nl)
70      real wdtrain(size(icb)) ! (ncum)      real wdtrain(size(icb)) ! (ncum)
71      logical lwork(size(icb)) ! (ncum)      logical lwork(size(icb)) ! (ncum)
72    
# Line 86  contains Line 105  contains
105         ! and condensed water flux         ! and condensed water flux
106    
107         ! Calculate detrained precipitation         ! Calculate detrained precipitation
108         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 &
109              * (ep(il, i) * m(il, i) * clw(il, i) &              * (ep(il, i) * m(il, i) * clw(il, i) &
110              + 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.) &
111              * ment(il, :i - 1, i)))              * ment(il, :i - 1, i)))
# Line 94  contains Line 113  contains
113         ! Find rain water and evaporation using provisional         ! Find rain water and evaporation using provisional
114         ! estimates of qp(i) and qp(i - 1)         ! estimates of qp(i) and qp(i - 1)
115    
116         do il = 1, ncum         loop_horizontal: do il = 1, ncum
117            if (i <= inb(il) .and. lwork(il)) then            if (i <= inb(il) .and. lwork(il)) then
118               wt(il, i) = 45.               wt(il, i) = 45.
119    
# Line 128  contains Line 147  contains
147               afac = max(afac, 0.)               afac = max(afac, 0.)
148               bfac = 1. / (sigd * wt(il, i))               bfac = 1. / (sigd * wt(il, i))
149    
150               ! Prise en compte de la variation progressive de sigt dans               if (i <= icb(il)) then
151               ! les couches icb et icb - 1:                  ! Prise en compte de la variation progressive de sigt dans
152               ! pour plcl < ph(i + 1), pr1 = 0 et pr2 = 1                  ! les couches icb et icb - 1 :
153               ! pour plcl > ph(i), pr1 = 1 et pr2 = 0                  ! pour plcl <= ph(i + 1), pr1 = 0
154               ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion \`a cheval                  ! pour plcl >= ph(i), pr1 = 1
155               ! sur le nuage, et pr2 est la proportion sous la base du                  ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion
156               ! nuage.                  ! \`a cheval sur le nuage.
157               pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))                  pr1 = max(0., min(1., &
158               pr1 = max(0., min(1., pr1))                       (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))
159               pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))                  sigt = sigp * pr1 + 1. - pr1
160               pr2 = max(0., min(1., pr2))               else
161               sigt = sigp(il, i) * pr1 + pr2                  ! {i >= icb(il) + 1}
162                    sigt = sigp
163                 end if
164    
165               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
166               c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &               c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
167                    * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)                    * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
168    
169               if (c6 > 0.) then               if (c6 > 0.) then
170                  revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))                  revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
171                  evap(il, i) = sigt * afac * revap                  evap(il, i) = sigt * afac * revap
# Line 167  contains Line 189  contains
189                  ! difference equation for downdraft theta and mass                  ! difference equation for downdraft theta and mass
190                  ! flux from two simultaneous differential equations                  ! flux from two simultaneous differential equations
191    
192                  amfac = sigd * sigd * 70. * ph(il, i) &                  amfac = sigd**2 * 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) &
                      * (p(il, i - 1) - p(il, i)) &  
193                       * (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))
194                  amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &                  amp2 = abs(mp(il, i + 1)**2 - mp(il, i)**2)
                      * mp(il, i))  
195    
196                  if (amp2 > 0.1 * amfac) then                  if (amp2 > 0.1 * amfac) then
197                     xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))                     xf = 100. * sigd**3 * (ph(il, i) - ph(il, i + 1))
198                     tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &                     tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
199                          * t(il, i) / (lvcp(il, i) * sigd * th(il, i))                          * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
200                     af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv                     af = xf * tf + mp(il, i + 1)**2 * tinv
201                     bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &                     bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
202                          * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &                          * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
203                          - p(il, i)) * xf * tevap                          - p(il, i)) * xf * tevap
# Line 201  contains Line 221  contains
221    
222                     mp(il, i) = max(0., mp(il, i))                     mp(il, i) = max(0., mp(il, i))
223    
224                     ! Il y a vraisemblablement une erreur dans la                     ! Il y a vraisemblablement une erreur dans la ligne
225                     ! ligne suivante : il faut diviser par (mp(il,                     ! suivante : il faut diviser par (mp(il, i) * sigd
226                     ! i) * sigd * grav) et non par (mp(il, i) + sigd                     ! * rg) et non par (mp(il, i) + sigd * 0.1).  Et il
227                     ! * 0.1).  Et il faut bien revoir les facteurs                     ! faut bien revoir les facteurs 100.
                    ! 100.  
228                     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &                     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &
229                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
230                          - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &                          - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &
# Line 256  contains Line 275  contains
275                  qp(il, i) = max(qp(il, i), 0.)                  qp(il, i) = max(qp(il, i), 0.)
276               endif               endif
277            endif            endif
278         end do         end do loop_horizontal
279      end DO downdraft_loop      end DO downdraft_loop
280    
281    end SUBROUTINE cv30_unsat    end SUBROUTINE cv30_unsat
Legend:
Removed from v.194  
changed lines
  Added in v.199
  ViewVC Help
Powered by ViewVC 1.1.21