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	guez	185	module cv30_undilute2_m
2	guez	47
3	guez	183	implicit none
4	guez	47
5	guez	183	contains
6	guez	47
7	guez	198	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	guez	47
10	guez	187	! 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	guez	47
15	guez	183	! Vertical profile of buoyancy computed here (use of buoybase).
16	guez	47
17	guez	213	use conf_phys_m, only: epmax
18	guez	195	use cv30_param_m, only: minorig, nl
19	guez	311	use cv_thermo, only: clmcpv, eps
20	guez	187	USE dimphy, ONLY: klon, klev
21	guez	201	use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rv
22	guez	47
23	guez	332	integer, intent(in):: icb(:) ! (ncum) {2 <= icb <= nl - 3}
24
25			integer, intent(in):: icbs(:) ! (ncum)
26	guez	187	! icbs is the first level above LCL (may differ from icb)
27
28	guez	198	real, intent(in):: tnk(:), qnk(:), gznk(:) ! (klon)
29	guez	187	real, intent(in):: t(klon, klev), qs(klon, klev), gz(klon, klev)
30			real, intent(in):: p(klon, klev), h(klon, klev)
31	guez	201	real, intent(in):: tv(klon, klev)
32			real, intent(in):: lv(:, :) ! (ncum, nl)
33	guez	196	real, intent(in):: pbase(:), buoybase(:), plcl(:) ! (ncum)
34	guez	187
35	guez	183	! outputs:
36	guez	187	integer, intent(out):: inb(:) ! (ncum)
37			! first model level above the level of neutral buoyancy of the
38	guez	192	! parcel (1 <= inb <= nl - 1)
39	guez	187
40			real tp(klon, klev), tvp(klon, klev), clw(klon, klev)
41	guez	183	! condensed water not removed from tvp
42	guez	195	real hp(klon, klev), ep(klon, klev)
43	guez	187	real buoy(klon, klev)
44	guez	47
45	guez	183	! Local:
46	guez	195
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	guez	183	integer i, k
60			real tg, qg, ahg, alv, s, tc, es, denom
61			real pden
62	guez	187	real ah0(klon)
63	guez	47
64	guez	183	!---------------------------------------------------------------------
65	guez	47
66	guez	195	ncum = size(icb)
67
68	guez	183	! SOME INITIALIZATIONS
69	guez	47
70	guez	183	do k = 1, nl
71			do i = 1, ncum
72	guez	195	ep(i, k) = 0.
73	guez	183	end do
74			end do
75	guez	47
76	guez	183	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
77	guez	47
78	guez	183	! The procedure is to solve the equation.
79			! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk.
80	guez	47
81	guez	183	! Calculate certain parcel quantities, including static energy
82	guez	47
83	guez	183	do i = 1, ncum
84	guez	201	ah0(i) = (rcpd * (1. - qnk(i)) + rcw * qnk(i)) * tnk(i) &
85	guez	198	+ qnk(i) * (rlvtt - clmcpv * (tnk(i) - 273.15)) + gznk(i)
86	guez	183	end do
87	guez	47
88	guez	183	! Find lifted parcel quantities above cloud base
89	guez	47
90	guez	183	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	guez	198	alv = rlvtt - clmcpv * (t(i, k) - 273.15)
96	guez	47
97	guez	183	! First iteration.
98	guez	47
99	guez	201	s = rcpd * (1. - qnk(i)) + rcw * qnk(i) &
100			+ alv * alv * qg / (rv * t(i, k) * t(i, k))
101	guez	183	s = 1. / s
102	guez	47
103	guez	201	ahg = rcpd * tg + (rcw - rcpd) * qnk(i) * tg + alv * qg + gz(i, k)
104	guez	183	tg = tg + s * (ah0(i) - ahg)
105	guez	47
106	guez	183	tc = tg - 273.15
107			denom = 243.5 + tc
108	guez	195	denom = MAX(denom, 1.)
109	guez	47
110	guez	183	es = 6.112 * exp(17.67 * tc / denom)
111	guez	47
112	guez	183	qg = eps * es / (p(i, k) - es * (1. - eps))
113	guez	47
114	guez	183	! Second iteration.
115	guez	47
116	guez	201	ahg = rcpd * tg + (rcw - rcpd) * qnk(i) * tg + alv * qg + gz(i, k)
117	guez	183	tg = tg + s * (ah0(i) - ahg)
118	guez	47
119	guez	183	tc = tg - 273.15
120			denom = 243.5 + tc
121	guez	195	denom = MAX(denom, 1.)
122	guez	47
123	guez	183	es = 6.112 * exp(17.67 * tc / denom)
124
125			qg = eps * es / (p(i, k) - es * (1. - eps))
126
127	guez	198	alv = rlvtt - clmcpv * (t(i, k) - 273.15)
128	guez	183
129			! no approximation:
130			tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) &
131	guez	201	/ (rcpd + (rcw - rcpd) * qnk(i))
132	guez	183
133			clw(i, k) = qnk(i) - qg
134	guez	195	clw(i, k) = max(0., clw(i, k))
135	guez	183	! 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	guez	195	! SET THE PRECIPITATION EFFICIENCIES
142			! It MAY BE a FUNCTION OF TP(I), P(I) AND CLW(I)
143	guez	183	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	guez	195	ep(i, k) = max(ep(i, k), 0.)
148	guez	192	ep(i, k) = min(ep(i, k), epmax)
149	guez	183	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	guez	186	tp(i, nl + 1) = tp(i, nl)
159	guez	183	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	guez	187	! Compute inb:
182	guez	183
183	guez	187	inb = nl - 1
184	guez	183
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	guez	186	do k = 1, nl + 1
196	guez	183	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	guez	198	if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, minorig) &
204	guez	201	+ (lv(i, k) + (rcpd - rcpv) * t(i, k)) * ep(i, k) * clw(i, k)
205	guez	183	end do
206			end do
207
208	guez	185	end SUBROUTINE cv30_undilute2
209	guez	183
210	guez	185	end module cv30_undilute2_m