phylmd/CV30_routines/cv30_undilute2.f90

module cv30_undilute2_m

  implicit none

contains

  SUBROUTINE cv30_undilute2(icb, icbs, 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 conf_phys_m, only: epmax
    use cv30_param_m, only: minorig, nl
    use cv_thermo, only: clmcpv, eps
    USE dimphy, ONLY: klon, klev
    use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rv

    integer, intent(in):: icb(:) ! (ncum) {2 <= icb <= nl - 3}

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

    real, intent(in):: tnk(:), qnk(:), 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)
    real, intent(in):: lv(:, :) ! (ncum, nl)
    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) = (rcpd * (1. - qnk(i)) + rcw * qnk(i)) * tnk(i) &
            + qnk(i) * (rlvtt - 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 = rlvtt - clmcpv * (t(i, k) - 273.15)

             ! First iteration.

             s = rcpd * (1. - qnk(i)) + rcw * qnk(i) &
                  + alv * alv * qg / (rv * t(i, k) * t(i, k))
             s = 1. / s

             ahg = rcpd * tg + (rcw - rcpd) * 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 = rcpd * tg + (rcw - rcpd) * 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 = rlvtt - clmcpv * (t(i, k) - 273.15)

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