phylmd/CV30_routines/cv30_undilute2.f

module cv30_undilute2_m

  implicit none

contains

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

    ! Undilute (adiabatic) updraft, second part. 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: minorig, nl
    use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
    USE dimphy, ONLY: klon, klev

    integer, intent(in):: icb(:), icbs(:) ! (ncum)
    ! icbs is the first level above LCL (may differ from icb)

    integer, intent(in):: nk(klon)
    real, intent(in):: tnk(klon), qnk(klon), gznk(klon)
    real, intent(in):: t(klon, klev), qs(klon, klev), gz(klon, klev)
    real, intent(in):: p(klon, klev), h(klon, klev)
    real, intent(in):: tv(klon, klev), lv(klon, klev)
    real, intent(in):: pbase(:), buoybase(:), plcl(:) ! (ncum)

    ! outputs:
    integer, intent(out):: inb(:) ! (ncum)
    ! first model level above the level of neutral buoyancy of the
    ! parcel (1 <= inb <= nl - 1)

    real tp(klon, klev), tvp(klon, klev), clw(klon, klev)
    ! condensed water not removed from tvp
    real hp(klon, klev), ep(klon, klev)
    real buoy(klon, klev)

    ! Local:

    integer ncum

    real, parameter:: pbcrit = 150.
    ! critical cloud depth (mbar) beneath which the precipitation
    ! efficiency is assumed to be zero

    real, parameter:: ptcrit = 500.
    ! cloud depth (mbar) above which the precipitation efficiency is
    ! assumed to be unity

    real, parameter:: dtovsh = - 0.2 ! dT for overshoot

    integer i, k
    real tg, qg, ahg, alv, s, tc, es, denom
    real pden
    real ah0(klon)

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

    ncum = size(icb)

    ! SOME INITIALIZATIONS

    do k = 1, nl
       do i = 1, ncum
          ep(i, k) = 0.
       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.)

             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.)

             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., 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
    ! It MAY BE a FUNCTION 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) = max(ep(i, k), 0.)
          ep(i, k) = min(ep(i, k), epmax)
       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

    ! Compute inb:

    inb = nl - 1

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