phylmd/CV30_routines/cv30_undilute2.f

module cv3_undilute2_m

  implicit none

contains

  SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
       qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
       ep, sigp, buoy)

    ! Purpose: find the rest of the lifted parcel temperatures;
    ! compute the precipitation efficiencies and the fraction of
    ! precipitation falling outside of cloud; find the level of
    ! neutral buoyancy.

    ! Vertical profile of buoyancy computed here (use of buoybase).

    use conema3_m, only: epmax
    use cv3_param_m, only: dtovsh, minorig, nl, nlp, pbcrit, ptcrit, spfac
    use cvthermo, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv

    ! inputs:
    integer, intent(in):: nloc, ncum, nd
    integer icb(nloc), icbs(nloc), nk(nloc)
    ! icbs (input) is the first level above LCL (may differ from icb)
    real tnk(nloc), qnk(nloc), gznk(nloc)
    real t(nloc, nd), qs(nloc, nd), gz(nloc, nd)
    real p(nloc, nd), h(nloc, nd)
    real tv(nloc, nd), lv(nloc, nd)
    real pbase(nloc), buoybase(nloc), plcl(nloc)

    ! outputs:
    integer inb(nloc)
    real tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
    ! condensed water not removed from tvp
    real hp(nloc, nd), ep(nloc, nd), sigp(nloc, nd)
    real buoy(nloc, nd)

    ! Local:
    integer i, k
    real tg, qg, ahg, alv, s, tc, es, denom
    real pden
    real ah0(nloc)

    !---------------------------------------------------------------------

    ! SOME INITIALIZATIONS

    do k = 1, nl
       do i = 1, ncum
          ep(i, k) = 0.0
          sigp(i, k) = spfac
       end do
    end do

    ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES

    ! The procedure is to solve the equation.
    ! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.

    ! Calculate certain parcel quantities, including static energy

    do i = 1, ncum
       ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) &
            + qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i)
    end do

    ! Find lifted parcel quantities above cloud base

    do k = minorig + 1, nl
       do i = 1, ncum
          if (k >= (icbs(i) + 1)) then
             tg = t(i, k)
             qg = qs(i, k)
             alv = lv0 - clmcpv * (t(i, k) - 273.15)

             ! First iteration.

             s = cpd * (1. - qnk(i)) + cl * qnk(i) &
                  + alv * alv * qg / (rrv * t(i, k) * t(i, k))
             s = 1. / s

             ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
             tg = tg + s * (ah0(i) - ahg)

             tc = tg - 273.15
             denom = 243.5 + tc
             denom = MAX(denom, 1.0)

             es = 6.112 * exp(17.67 * tc / denom)

             qg = eps * es / (p(i, k) - es * (1. - eps))

             ! Second iteration.

             ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
             tg = tg + s * (ah0(i) - ahg)

             tc = tg - 273.15
             denom = 243.5 + tc
             denom = MAX(denom, 1.0)

             es = 6.112 * exp(17.67 * tc / denom)

             qg = eps * es / (p(i, k) - es * (1. - eps))

             alv = lv0 - clmcpv * (t(i, k) - 273.15)

             ! no approximation:
             tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) &
                  / (cpd + (cl - cpd) * qnk(i))

             clw(i, k) = qnk(i) - qg
             clw(i, k) = max(0.0, clw(i, k))
             ! qg utilise au lieu du vrai mixing ratio rg:
             tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing
          endif
       end do
    end do

    ! SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
    ! PRECIPITATION FALLING OUTSIDE OF CLOUD
    ! THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
    do k = 1, nl
       do i = 1, ncum
          pden = ptcrit - pbcrit
          ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax
          ep(i, k) = amax1(ep(i, k), 0.0)
          ep(i, k) = amin1(ep(i, k), epmax)
          sigp(i, k) = spfac
       end do
    end do

    ! CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
    ! VIRTUAL TEMPERATURE

    ! tvp est calcule en une seule fois, et sans retirer
    ! l'eau condensee (~> reversible CAPE)
    do i = 1, ncum
       tp(i, nlp) = tp(i, nl)
    end do

    ! EFFECTIVE VERTICAL PROFILE OF BUOYANCY:

    ! first estimate of buoyancy:
    do i = 1, ncum
       do k = 1, nl
          buoy(i, k) = tvp(i, k) - tv(i, k)
       end do
    end do

    ! set buoyancy = buoybase for all levels below base
    ! for safety, set buoy(icb) = buoybase
    do i = 1, ncum
       do k = 1, nl
          if ((k >= icb(i)) .and. (k <= nl) .and. (p(i, k) >= pbase(i))) then
             buoy(i, k) = buoybase(i)
          endif
       end do
       buoy(icb(i), k) = buoybase(i)
    end do

    ! FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
    ! LEVEL OF NEUTRAL BUOYANCY

    do i = 1, ncum
       inb(i) = nl - 1
    end do

    do i = 1, ncum
       do k = 1, nl - 1
          if ((k >= icb(i)) .and. (buoy(i, k) < dtovsh)) then
             inb(i) = MIN(inb(i), k)
          endif
       end do
    end do

    ! CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL

    do k = 1, nlp
       do i = 1, ncum
          hp(i, k) = h(i, k)
       enddo
    enddo

    do k = minorig + 1, nl
       do i = 1, ncum
          if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, nk(i)) &
               + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
       end do
    end do

  end SUBROUTINE cv3_undilute2

end module cv3_undilute2_m
1	guez	183	module cv3_undilute2_m
2	guez	47
3	guez	183	implicit none
4	guez	47
5	guez	183	contains
6	guez	47
7	guez	183	SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
8			qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
9			ep, sigp, buoy)
10	guez	47
11	guez	183	! Purpose: find the rest of the lifted parcel temperatures;
12			! compute the precipitation efficiencies and the fraction of
13			! precipitation falling outside of cloud; find the level of
14			! neutral buoyancy.
15	guez	47
16	guez	183	! Vertical profile of buoyancy computed here (use of buoybase).
17	guez	47
18	guez	183	use conema3_m, only: epmax
19			use cv3_param_m, only: dtovsh, minorig, nl, nlp, pbcrit, ptcrit, spfac
20			use cvthermo, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
21	guez	47
22	guez	183	! inputs:
23			integer, intent(in):: nloc, ncum, nd
24			integer icb(nloc), icbs(nloc), nk(nloc)
25			! icbs (input) is the first level above LCL (may differ from icb)
26			real tnk(nloc), qnk(nloc), gznk(nloc)
27			real t(nloc, nd), qs(nloc, nd), gz(nloc, nd)
28			real p(nloc, nd), h(nloc, nd)
29			real tv(nloc, nd), lv(nloc, nd)
30			real pbase(nloc), buoybase(nloc), plcl(nloc)
31	guez	47
32	guez	183	! outputs:
33			integer inb(nloc)
34			real tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
35			! condensed water not removed from tvp
36			real hp(nloc, nd), ep(nloc, nd), sigp(nloc, nd)
37			real buoy(nloc, nd)
38	guez	47
39	guez	183	! Local:
40			integer i, k
41			real tg, qg, ahg, alv, s, tc, es, denom
42			real pden
43			real ah0(nloc)
44	guez	47
45	guez	183	!---------------------------------------------------------------------
46	guez	47
47	guez	183	! SOME INITIALIZATIONS
48	guez	47
49	guez	183	do k = 1, nl
50			do i = 1, ncum
51			ep(i, k) = 0.0
52			sigp(i, k) = spfac
53			end do
54			end do
55	guez	47
56	guez	183	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
57	guez	47
58	guez	183	! The procedure is to solve the equation.
59			! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.
60	guez	47
61	guez	183	! Calculate certain parcel quantities, including static energy
62	guez	47
63	guez	183	do i = 1, ncum
64			ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) &
65			+ qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i)
66			end do
67	guez	47
68	guez	183	! Find lifted parcel quantities above cloud base
69	guez	47
70	guez	183	do k = minorig + 1, nl
71			do i = 1, ncum
72			if (k >= (icbs(i) + 1)) then
73			tg = t(i, k)
74			qg = qs(i, k)
75			alv = lv0 - clmcpv * (t(i, k) - 273.15)
76	guez	47
77	guez	183	! First iteration.
78	guez	47
79	guez	183	s = cpd * (1. - qnk(i)) + cl * qnk(i) &
80			+ alv * alv * qg / (rrv * t(i, k) * t(i, k))
81			s = 1. / s
82	guez	47
83	guez	183	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
84			tg = tg + s * (ah0(i) - ahg)
85	guez	47
86	guez	183	tc = tg - 273.15
87			denom = 243.5 + tc
88			denom = MAX(denom, 1.0)
89	guez	47
90	guez	183	es = 6.112 * exp(17.67 * tc / denom)
91	guez	47
92	guez	183	qg = eps * es / (p(i, k) - es * (1. - eps))
93	guez	47
94	guez	183	! Second iteration.
95	guez	47
96	guez	183	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
97			tg = tg + s * (ah0(i) - ahg)
98	guez	47
99	guez	183	tc = tg - 273.15
100			denom = 243.5 + tc
101			denom = MAX(denom, 1.0)
102	guez	47
103	guez	183	es = 6.112 * exp(17.67 * tc / denom)
104
105			qg = eps * es / (p(i, k) - es * (1. - eps))
106
107			alv = lv0 - clmcpv * (t(i, k) - 273.15)
108
109			! no approximation:
110			tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) &
111			/ (cpd + (cl - cpd) * qnk(i))
112
113			clw(i, k) = qnk(i) - qg
114			clw(i, k) = max(0.0, clw(i, k))
115			! qg utilise au lieu du vrai mixing ratio rg:
116			tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing
117			endif
118			end do
119			end do
120
121			! SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
122			! PRECIPITATION FALLING OUTSIDE OF CLOUD
123			! THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
124			do k = 1, nl
125			do i = 1, ncum
126			pden = ptcrit - pbcrit
127			ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax
128			ep(i, k) = amax1(ep(i, k), 0.0)
129			ep(i, k) = amin1(ep(i, k), epmax)
130			sigp(i, k) = spfac
131			end do
132			end do
133
134			! CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
135			! VIRTUAL TEMPERATURE
136
137			! tvp est calcule en une seule fois, et sans retirer
138			! l'eau condensee (~> reversible CAPE)
139			do i = 1, ncum
140			tp(i, nlp) = tp(i, nl)
141			end do
142
143			! EFFECTIVE VERTICAL PROFILE OF BUOYANCY:
144
145			! first estimate of buoyancy:
146			do i = 1, ncum
147			do k = 1, nl
148			buoy(i, k) = tvp(i, k) - tv(i, k)
149			end do
150			end do
151
152			! set buoyancy = buoybase for all levels below base
153			! for safety, set buoy(icb) = buoybase
154			do i = 1, ncum
155			do k = 1, nl
156			if ((k >= icb(i)) .and. (k <= nl) .and. (p(i, k) >= pbase(i))) then
157			buoy(i, k) = buoybase(i)
158			endif
159			end do
160			buoy(icb(i), k) = buoybase(i)
161			end do
162
163			! FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
164			! LEVEL OF NEUTRAL BUOYANCY
165
166			do i = 1, ncum
167			inb(i) = nl - 1
168			end do
169
170			do i = 1, ncum
171			do k = 1, nl - 1
172			if ((k >= icb(i)) .and. (buoy(i, k) < dtovsh)) then
173			inb(i) = MIN(inb(i), k)
174			endif
175			end do
176			end do
177
178			! CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
179
180			do k = 1, nlp
181			do i = 1, ncum
182			hp(i, k) = h(i, k)
183			enddo
184			enddo
185
186			do k = minorig + 1, nl
187			do i = 1, ncum
188			if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, nk(i)) &
189			+ (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
190			end do
191			end do
192
193			end SUBROUTINE cv3_undilute2
194
195			end module cv3_undilute2_m