/[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 192 by guez, Thu May 12 13:00:07 2016 UTC revision 193 by guez, Thu May 12 13:22:19 2016 UTC
# Line 41  contains Line 41  contains
41    
42      ! Local:      ! Local:
43      integer ncum      integer ncum
44      integer i, j, il, num1      integer i, j, il, imax
45      real tinv, delti      real tinv, delti
46      real awat, afac, afac1, afac2, bfac      real awat, afac, afac1, afac2, bfac
47      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real pr1, pr2, sigt, b6, c6, revap, tevap, delth
48      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
49      real ampmax      real ampmax
50      real lvcp(klon, klev)      real lvcp(klon, klev)
51      real wdtrain(size(icb))      real wdtrain(size(icb)) ! (ncum)
52      logical lwork(size(icb))      logical lwork(size(icb)) ! (ncum)
53    
54      !------------------------------------------------------      !------------------------------------------------------
55    
# Line 71  contains Line 71  contains
71         enddo         enddo
72      enddo      enddo
73    
74      ! check whether ep(inb) = 0, if so, skip precipitating      ! Check whether ep(inb) = 0. If so, skip precipitating downdraft
75      ! downdraft calculation      ! calculation.
76    
77      forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4      forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
78    
79      wdtrain = 0.      imax = nl - 1
80        do while (.not. any(inb >= imax .and. lwork) .and. imax >= 1)
81           imax = imax - 1
82        end do
83    
84        downdraft_loop: DO i = imax, 1, - 1
85           ! Integrate liquid water equation to find condensed water
86           ! and condensed water flux
87    
88      downdraft_loop: DO i = nl - 1, 1, - 1         ! Calculate detrained precipitation
        num1 = 0  
89    
90         do il = 1, ncum         do il = 1, ncum
91            if (i <= inb(il) .and. lwork(il)) num1 = num1 + 1            if (i <= inb(il) .and. lwork(il)) then
92                 wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
93              endif
94         enddo         enddo
95    
96         if (num1 > 0) then         if (i > 1) then
97            ! integrate liquid water equation to find condensed water            do j = 1, i - 1
98            ! and condensed water flux               do il = 1, ncum
99                    if (i <= inb(il) .and. lwork(il)) then
100            ! calculate detrained precipitation                     awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
101                       awat = max(awat, 0.)
102            do il = 1, ncum                     wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
103               if (i <= inb(il) .and. lwork(il)) then                  endif
104                  wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)               enddo
105               endif            end do
106            enddo         endif
107    
108            if (i > 1) then         ! Find rain water and evaporation using provisional
109               do j = 1, i - 1         ! estimates of qp(i) and qp(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  
110    
111            ! find rain water and evaporation using provisional         do il = 1, ncum
112            ! estimates of qp(i) and qp(i - 1)            if (i <= inb(il) .and. lwork(il)) then
113                 wt(il, i) = 45.
114    
115            do il = 1, ncum               if (i < inb(il)) then
116               if (i <= inb(il) .and. lwork(il)) then                  qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, i)) &
117                  wt(il, i) = 45.                       + gz(il, i + 1) - gz(il, i)) / lv(il, i)
118                    qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
119                  if (i < inb(il)) then               endif
                    qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) &  
                         - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)  
                    qp(il, i) = 0.5 * (qp(il, i) + q(il, i))  
                 endif  
120    
121                  qp(il, i) = max(qp(il, i), 0.)               qp(il, i) = max(qp(il, i), 0.)
122                  qp(il, i) = min(qp(il, i), qs(il, i))               qp(il, i) = min(qp(il, i), qs(il, i))
123                  qp(il, inb(il)) = q(il, inb(il))               qp(il, inb(il)) = q(il, inb(il))
124    
125                 if (i == 1) then
126                    afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
127                         / (1e4 + 2000. * p(il, 1) * qs(il, 1))
128                 else
129                    qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) - t(il, i - 1)) &
130                         + gz(il, i) - gz(il, i - 1)) / lv(il, i)
131                    qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
132                    qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
133                    qp(il, i - 1) = max(qp(il, i - 1), 0.)
134                    afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &
135                         / (1e4 + 2000. * p(il, i) * qs(il, i))
136                    afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) &
137                         / (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1))
138                    afac = 0.5 * (afac1 + afac2)
139                 endif
140    
141                  if (i == 1) then               if (i == inb(il)) afac = 0.
142                     afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &               afac = max(afac, 0.)
143                          / (1e4 + 2000. * p(il, 1) * qs(il, 1))               bfac = 1. / (sigd * wt(il, i))
144                  else  
145                     qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) &               ! Prise en compte de la variation progressive de sigt dans
146                          - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)               ! les couches icb et icb - 1:
147                     qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))               ! pour plcl < ph(i + 1), pr1 = 0 et pr2 = 1
148                     qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))               ! pour plcl > ph(i), pr1 = 1 et pr2 = 0
149                     qp(il, i - 1) = max(qp(il, i - 1), 0.)               ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion \`a cheval
150                     afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &               ! sur le nuage, et pr2 est la proportion sous la base du
151                          / (1e4 + 2000. * p(il, i) * qs(il, i))               ! nuage.
152                     afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) &               pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))
153                          / (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1))               pr1 = max(0., min(1., pr1))
154                     afac = 0.5 * (afac1 + afac2)               pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))
155                 pr2 = max(0., min(1., pr2))
156                 sigt = sigp(il, i) * pr1 + pr2
157    
158                 b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac
159                 c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
160                      * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
161                 if (c6 > 0.) then
162                    revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
163                    evap(il, i) = sigt * afac * revap
164                    water(il, i) = revap * revap
165                 else
166                    evap(il, i) = - evap(il, i + 1) + 0.02 * (wdtrain(il) + sigd &
167                         * wt(il, i) * water(il, i + 1)) / (sigd * (ph(il, i) &
168                         - ph(il, i + 1)))
169                 end if
170    
171                 ! Calculate precipitating downdraft mass flux under
172                 ! hydrostatic approximation
173    
174                 test_above_surface: if (i /= 1) then
175                    tevap = max(0., evap(il, i))
176                    delth = max(0.001, (th(il, i) - th(il, i - 1)))
177                    mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
178                         * (p(il, i - 1) - p(il, i)) / delth
179    
180                    ! If hydrostatic assumption fails, solve cubic
181                    ! difference equation for downdraft theta and mass
182                    ! flux from two simultaneous differential equations
183    
184                    amfac = sigd * sigd * 70. * ph(il, i) &
185                         * (p(il, i - 1) - p(il, i)) &
186                         * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
187                    amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
188                         * mp(il, i))
189    
190                    if (amp2 > 0.1 * amfac) then
191                       xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))
192                       tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
193                            * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
194                       af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
195                       bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
196                            * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
197                            - p(il, i)) * xf * tevap
198                       fac2 = 1.
199                       if (bf < 0.) fac2 = - 1.
200                       bf = abs(bf)
201                       ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
202    
203                       if (ur >= 0.) then
204                          sru = sqrt(ur)
205                          fac = 1.
206                          if ((0.5 * bf - sru) < 0.) fac = - 1.
207                          mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf &
208                               + sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv
209                       else
210                          d = atan(2. * sqrt(- ur) / (bf + 1e-28))
211                          if (fac2 < 0.)d = 3.14159 - d
212                          mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &
213                               * cos(d * tinv)
214                       endif
215    
216                       mp(il, i) = max(0., mp(il, i))
217    
218                       ! Il y a vraisemblablement une erreur dans la
219                       ! ligne suivante : il faut diviser par (mp(il,
220                       ! i) * sigd * grav) et non par (mp(il, i) + sigd
221                       ! * 0.1).  Et il faut bien revoir les facteurs
222                       ! 100.
223                       b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &
224                            * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
225                            - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &
226                            * th(il, i))
227                       b(il, i - 1) = max(b(il, i - 1), 0.)
228                  endif                  endif
229    
230                  if (i == inb(il)) afac = 0.                  ! Limit magnitude of mp(i) to meet CFL condition:
231                  afac = max(afac, 0.)                  ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
232                  bfac = 1. / (sigd * wt(il, i))                  amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
233                    ampmax = min(ampmax, amp2)
234                  ! prise en compte de la variation progressive de sigt dans                  mp(il, i) = min(mp(il, i), ampmax)
235                  ! les couches icb et icb - 1:  
236                  ! pour plcl < ph(i + 1), pr1 = 0 & pr2 = 1                  ! Force mp to decrease linearly to zero between cloud
237                  ! pour plcl > ph(i), pr1 = 1 & pr2 = 0                  ! base and the surface:
238                  ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion a cheval                  if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
239                  ! sur le nuage, et pr2 est la proportion sous la base du                       * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
240                  ! nuage.               endif test_above_surface
241                  pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))  
242                  pr1 = max(0., min(1., pr1))               ! Find mixing ratio of precipitating downdraft
243                  pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))  
244                  pr2 = max(0., min(1., pr2))               if (i /= inb(il)) then
245                  sigt = sigp(il, i) * pr1 + pr2                  qp(il, i) = q(il, i)
246    
247                  b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt &                  if (mp(il, i) > mp(il, i + 1)) then
248                       * afac                     qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
249                  c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &                          * (mp(il, i) - mp(il, i + 1)) + 100. * ginv * 0.5 &
250                       * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)                          * sigd * (ph(il, i) - ph(il, i + 1)) &
251                  if (c6 > 0.) then                          * (evap(il, i + 1) + evap(il, i))
252                     revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))                     qp(il, i) = qp(il, i) / mp(il, i)
253                     evap(il, i) = sigt * afac * revap                     up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
254                     water(il, i) = revap * revap                          * (mp(il, i) - mp(il, i + 1))
255                       up(il, i) = up(il, i) / mp(il, i)
256                       vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
257                            * (mp(il, i) - mp(il, i + 1))
258                       vp(il, i) = vp(il, i) / mp(il, i)
259                  else                  else
260                     evap(il, i) = - evap(il, i + 1) + 0.02 * (wdtrain(il) &                     if (mp(il, i + 1) > 1e-16) then
261                          + sigd * wt(il, i) * water(il, i + 1)) &                        qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
262                          / (sigd * (ph(il, i) - ph(il, i + 1)))                             * (ph(il, i) - ph(il, i + 1)) * (evap(il, i + 1) &
263                  end if                             + evap(il, i)) / mp(il, i + 1)
264                          up(il, i) = up(il, i + 1)
265                  ! calculate precipitating downdraft mass flux under                        vp(il, i) = vp(il, i + 1)
                 ! hydrostatic approximation  
   
                 if (i /= 1) then  
                    tevap = max(0., evap(il, i))  
                    delth = max(0.001, (th(il, i) - th(il, i - 1)))  
                    mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &  
                         * (p(il, i - 1) - p(il, i)) / delth  
   
                    ! If hydrostatic assumption fails, solve cubic  
                    ! difference equation for downdraft theta and mass  
                    ! flux from two simultaneous differential equations  
   
                    amfac = sigd * sigd * 70. * ph(il, i) &  
                         * (p(il, i - 1) - p(il, i)) &  
                         * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))  
                    amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &  
                         * mp(il, i))  
   
                    if (amp2 > 0.1 * amfac) then  
                       xf = 100. * sigd * sigd * sigd * (ph(il, i) &  
                            - ph(il, i + 1))  
                       tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &  
                            * t(il, i) / (lvcp(il, i) * sigd * th(il, i))  
                       af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv  
                       bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &  
                            * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &  
                            - p(il, i)) * xf * tevap  
                       fac2 = 1.  
                       if (bf < 0.) fac2 = - 1.  
                       bf = abs(bf)  
                       ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv  
   
                       if (ur >= 0.) then  
                          sru = sqrt(ur)  
                          fac = 1.  
                          if ((0.5 * bf - sru) < 0.) fac = - 1.  
                          mp(il, i) = mp(il, i + 1) * tinv &  
                               + (0.5 * bf + sru)**tinv &  
                               + fac * (abs(0.5 * bf - sru))**tinv  
                       else  
                          d = atan(2. * sqrt(- ur) / (bf + 1e-28))  
                          if (fac2 < 0.)d = 3.14159 - d  
                          mp(il, i) = mp(il, i + 1) * tinv + 2. &  
                               * sqrt(af * tinv) * cos(d * tinv)  
                       endif  
   
                       mp(il, i) = max(0., mp(il, i))  
   
                       ! Il y a vraisemblablement une erreur dans la  
                       ! ligne suivante : il faut diviser par (mp(il,  
                       ! i) * sigd * grav) et non par (mp(il, i) + sigd  
                       ! * 0.1).  Et il faut bien revoir les facteurs  
                       ! 100.  
                       b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &  
                            - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &  
                            - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &  
                            / (lvcp(il, i) * sigd * th(il, i))  
                       b(il, i - 1) = max(b(il, i - 1), 0.)  
                    endif  
   
                    ! limit magnitude of mp(i) to meet cfl condition  
                    ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti  
                    amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti  
                    ampmax = min(ampmax, amp2)  
                    mp(il, i) = min(mp(il, i), ampmax)  
   
                    ! force mp to decrease linearly to zero  
                    ! between cloud base and the surface  
                    if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &  
                         * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))  
                 endif ! i == 1  
   
                 ! find mixing ratio of precipitating downdraft  
   
                 if (i /= inb(il)) then  
                    qp(il, i) = q(il, i)  
   
                    if (mp(il, i) > mp(il, i + 1)) then  
                       qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &  
                            * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &  
                            * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &  
                            * (evap(il, i + 1) + evap(il, i))  
                       qp(il, i) = qp(il, i) / mp(il, i)  
                       up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &  
                            * (mp(il, i) - mp(il, i + 1))  
                       up(il, i) = up(il, i) / mp(il, i)  
                       vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &  
                            * (mp(il, i) - mp(il, i + 1))  
                       vp(il, i) = vp(il, i) / mp(il, i)  
                    else  
                       if (mp(il, i + 1) > 1e-16) then  
                          qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &  
                               * (ph(il, i) - ph(il, i + 1)) &  
                               * (evap(il, i + 1) + evap(il, i)) / mp(il, i + 1)  
                          up(il, i) = up(il, i + 1)  
                          vp(il, i) = vp(il, i + 1)  
                       endif  
266                     endif                     endif
   
                    qp(il, i) = min(qp(il, i), qs(il, i))  
                    qp(il, i) = max(qp(il, i), 0.)  
267                  endif                  endif
268    
269                    qp(il, i) = min(qp(il, i), qs(il, i))
270                    qp(il, i) = max(qp(il, i), 0.)
271               endif               endif
272            end do            endif
273         end if         end do
274      end DO downdraft_loop      end DO downdraft_loop
275    
276    end SUBROUTINE cv30_unsat    end SUBROUTINE cv30_unsat

Legend:
Removed from v.192  
changed lines
  Added in v.193

  ViewVC Help
Powered by ViewVC 1.1.21