phylmd/CV30_routines/cv30_undilute2.f

module cv30_undilute2_m

  implicit none

contains

  SUBROUTINE cv30_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 cv30_param_m, only: dtovsh, minorig, nl, 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, nl + 1) = 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, nl + 1
       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 cv30_undilute2

end module cv30_undilute2_m
1	module cv30_undilute2_m
2
3	implicit none
4
5	contains
6
7	SUBROUTINE cv30_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
11	! 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
16	! Vertical profile of buoyancy computed here (use of buoybase).
17
18	use conema3_m, only: epmax
19	use cv30_param_m, only: dtovsh, minorig, nl, pbcrit, ptcrit, spfac
20	use cvthermo, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
21
22	! 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
32	! 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
39	! Local:
40	integer i, k
41	real tg, qg, ahg, alv, s, tc, es, denom
42	real pden
43	real ah0(nloc)
44
45	!---------------------------------------------------------------------
46
47	! SOME INITIALIZATIONS
48
49	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
56	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
57
58	! The procedure is to solve the equation.
59	! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.
60
61	! Calculate certain parcel quantities, including static energy
62
63	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
68	! Find lifted parcel quantities above cloud base
69
70	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
77	! First iteration.
78
79	s = cpd * (1. - qnk(i)) + cl * qnk(i) &
80	+ alv * alv * qg / (rrv * t(i, k) * t(i, k))
81	s = 1. / s
82
83	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
84	tg = tg + s * (ah0(i) - ahg)
85
86	tc = tg - 273.15
87	denom = 243.5 + tc
88	denom = MAX(denom, 1.0)
89
90	es = 6.112 * exp(17.67 * tc / denom)
91
92	qg = eps * es / (p(i, k) - es * (1. - eps))
93
94	! Second iteration.
95
96	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k)
97	tg = tg + s * (ah0(i) - ahg)
98
99	tc = tg - 273.15
100	denom = 243.5 + tc
101	denom = MAX(denom, 1.0)
102
103	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, nl + 1) = 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, nl + 1
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 cv30_undilute2
194
195	end module cv30_undilute2_m