/[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 198 by guez, Tue May 31 16:17:35 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 201  contains Line 223  contains
223    
224                     mp(il, i) = max(0., mp(il, i))                     mp(il, i) = max(0., mp(il, i))
225    
226                     ! Il y a vraisemblablement une erreur dans la                     ! Il y a vraisemblablement une erreur dans la ligne
227                     ! ligne suivante : il faut diviser par (mp(il,                     ! suivante : il faut diviser par (mp(il, i) * sigd
228                     ! i) * sigd * grav) et non par (mp(il, i) + sigd                     ! * rg) et non par (mp(il, i) + sigd * 0.1).  Et il
229                     ! * 0.1).  Et il faut bien revoir les facteurs                     ! faut bien revoir les facteurs 100.
                    ! 100.  
230                     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)) &
231                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &                          * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
232                          - 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 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.194  
changed lines
  Added in v.198
  ViewVC Help
Powered by ViewVC 1.1.21