phylmd/CV30_routines/cv30_undilute2.f

module cv30_undilute2_m

  implicit none

contains

  SUBROUTINE cv30_undilute2(ncum, icb, icbs, nk, tnk, qnk, gznk, t, qs, gz, &
       p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, ep, sigp, &
       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: dtovsh, minorig, nl, pbcrit, ptcrit, spfac
    use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
    USE dimphy, ONLY: klon, klev

    integer, intent(in):: ncum

    integer, intent(in):: icb(klon), icbs(klon)
    ! 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(klon), buoybase(klon), plcl(klon)

    ! 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), sigp(klon, klev)
    real buoy(klon, klev)

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

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

    ! 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) = max(ep(i, k), 0.0)
          ep(i, k) = min(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

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