phylmd/CV30_routines/cv30_undilute1.f

module cv30_undilute1_m

  implicit none

contains

  SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, nk1, icb1, tp1, &
       tvp1, clw1, icbs1)

    ! UNDILUTE (ADIABATIC) UPDRAFT / 1st part
    ! (up through ICB1 + 1)
    ! Calculates the lifted parcel virtual temperature at nk1, the
    ! actual temperature, and the adiabatic liquid water content.

    ! Equivalent de TLIFT entre NK1 et ICB1+1 inclus

    ! Differences with convect4:
    ! - icbs1 is the first level above LCL (may differ from icb1)
    ! - in the iterations, used x(icbs1) instead x(icb1)
    ! - tvp1 is computed in only one time
    ! - icbs1: first level above Plcl1 (IMIN de TLIFT) in output
    ! - if icbs1=icb1, compute also tp1(icb1+1), tvp1(icb1+1) & clw1(icb1+1)

    use cv30_param_m, only: minorig, nl
    use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
    USE dimphy, ONLY: klev, klon

    ! inputs:
    integer, intent(in):: nk1(klon), icb1(klon)
    real, intent(in):: t1(klon, klev)
    real, intent(in):: q1(klon, klev), qs1(klon, klev), gz1(klon, klev)
    real, intent(in):: p1(klon, klev)
    real, intent(in):: plcl1(klon)

    ! outputs:
    real tp1(klon, klev), tvp1(klon, klev), clw1(klon, klev)

    ! local variables:
    integer i, k
    integer icbs1(klon), icbsmax2
    real tg, qg, alv, s, ahg, tc, denom, es
    real ah0(klon), cpp(klon)
    real tnk(klon), qnk(klon), gznk(klon), ticb(klon), gzicb(klon)
    real qsicb(klon)
    real cpinv(klon)

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

    !  Calculates the lifted parcel virtual temperature at nk1,
    !  the actual temperature, and the adiabatic
    !  liquid water content. The procedure is to solve the equation.
    ! cp*tp1+L*qp+phi=cp*tnk+L*qnk+gznk.

    do i=1, klon
       tnk(i)=t1(i, nk1(i))
       qnk(i)=q1(i, nk1(i))
       gznk(i)=gz1(i, nk1(i))
    end do

    ! *** Calculate certain parcel quantities, including static energy ***

    do i=1, klon
       ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) &
            +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
       cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
       cpinv(i)=1./cpp(i)
    end do

    ! *** Calculate lifted parcel quantities below cloud base ***

    do i=1, klon
       ! if icb1 is below LCL, start loop at ICB1+1:
       ! (icbs1 est le premier niveau au-dessus du LCL)
       icbs1(i)=MIN(max(icb1(i), 2), nl)
       if (plcl1(i) < p1(i, icbs1(i))) then
          icbs1(i)=MIN(icbs1(i)+1, nl)
       endif
    enddo

    do i=1, klon
       ticb(i)=t1(i, icbs1(i))
       gzicb(i)=gz1(i, icbs1(i))
       qsicb(i)=qs1(i, icbs1(i))
    enddo

    ! Re-compute icbsmax (icbsmax2):
    icbsmax2=2
    do i=1, klon
       icbsmax2=max(icbsmax2, icbs1(i))
    end do

    ! initialization outputs:

    do k=1, icbsmax2
       do i=1, klon
          tp1(i, k) = 0.0
          tvp1(i, k) = 0.0
          clw1(i, k) = 0.0
       enddo
    enddo

    ! tp1 and tvp1 below cloud base:

    do k=minorig, icbsmax2-1
       do i=1, klon
          tp1(i, k)=tnk(i)-(gz1(i, k)-gznk(i))*cpinv(i)
          tvp1(i, k)=tp1(i, k)*(1.+qnk(i)/eps-qnk(i))
       end do
    end do

    ! *** Find lifted parcel quantities above cloud base ***

    do i=1, klon
       tg=ticb(i)
       qg=qsicb(i)
       !debug alv=lv0-clmcpv*(ticb(i)-t0)
       alv=lv0-clmcpv*(ticb(i)-273.15)

       ! First iteration.

       s=cpd*(1.-qnk(i))+cl*qnk(i) &
            +alv*alv*qg/(rrv*ticb(i)*ticb(i))
       s=1./s

       ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
       tg=tg+s*(ah0(i)-ahg)

       !debug tc=tg-t0
       tc=tg-273.15
       denom=243.5+tc
       denom=MAX(denom, 1.0)

       es=6.112*exp(17.67*tc/denom)
       qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))

       ! Second iteration.

       ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
       tg=tg+s*(ah0(i)-ahg)

       !debug tc=tg-t0
       tc=tg-273.15
       denom=243.5+tc
       denom=MAX(denom, 1.0)

       es=6.112*exp(17.67*tc/denom)

       qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps))

       alv=lv0-clmcpv*(ticb(i)-273.15)

       ! no approximation:
       tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) &
            /(cpd+(cl-cpd)*qnk(i))

       clw1(i, icbs1(i))=qnk(i)-qg
       clw1(i, icbs1(i))=max(0.0, clw1(i, icbs1(i)))

       ! (qg utilise au lieu du vrai mixing ratio rg)
       tvp1(i, icbs1(i))=tp1(i, icbs1(i))*(1.+qg/eps-qnk(i))

    end do

    ! * icbs1 is the first level above the LCL:
    ! if plcl1<p1(icb1), then icbs1=icb1+1
    ! if plcl1>p1(icb1), then icbs1=icb1

    ! * the routine above computes tvp1 from minorig to icbs1 (included).

    ! * to compute buoybase (in cv30_trigger.F), both tvp1(icb1) and
    ! tvp1(icb1+1) must be known. This is the case if icbs1=icb1+1,
    ! but not if icbs1=icb1.

    ! * therefore, in the case icbs1=icb1, we compute tvp1 at level icb1+1
    ! (tvp1 at other levels will be computed in cv30_undilute2.F)

    do i=1, klon
       ticb(i)=t1(i, icb1(i)+1)
       gzicb(i)=gz1(i, icb1(i)+1)
       qsicb(i)=qs1(i, icb1(i)+1)
    enddo

    do i=1, klon
       tg=ticb(i)
       qg=qsicb(i)
       !debug alv=lv0-clmcpv*(ticb(i)-t0)
       alv=lv0-clmcpv*(ticb(i)-273.15)

       ! First iteration.

       s=cpd*(1.-qnk(i))+cl*qnk(i) &
            +alv*alv*qg/(rrv*ticb(i)*ticb(i))
       s=1./s

       ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
       tg=tg+s*(ah0(i)-ahg)

       !debug tc=tg-t0
       tc=tg-273.15
       denom=243.5+tc
       denom=MAX(denom, 1.0)

       es=6.112*exp(17.67*tc/denom)

       qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))

       ! Second iteration.

       ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i)
       tg=tg+s*(ah0(i)-ahg)

       !debug tc=tg-t0
       tc=tg-273.15
       denom=243.5+tc
       denom=MAX(denom, 1.0)

       es=6.112*exp(17.67*tc/denom)

       qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps))

       alv=lv0-clmcpv*(ticb(i)-273.15)

       ! no approximation:
       tp1(i, icb1(i)+1)=(ah0(i)-gz1(i, icb1(i)+1)-alv*qg) &
            /(cpd+(cl-cpd)*qnk(i))

       clw1(i, icb1(i)+1)=qnk(i)-qg
       clw1(i, icb1(i)+1)=max(0.0, clw1(i, icb1(i)+1))

       ! (qg utilise au lieu du vrai mixing ratio rg)
       tvp1(i, icb1(i)+1)=tp1(i, icb1(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
    end do

  end SUBROUTINE cv30_undilute1

end module cv30_undilute1_m
1	guez	195	module cv30_undilute1_m
2	guez	47
3	guez	195	implicit none
4	guez	47
5	guez	195	contains
6
7	guez	196	SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, nk1, icb1, tp1, &
8			tvp1, clw1, icbs1)
9	guez	195
10	guez	189	! UNDILUTE (ADIABATIC) UPDRAFT / 1st part
11	guez	196	! (up through ICB1 + 1)
12			! Calculates the lifted parcel virtual temperature at nk1, the
13	guez	189	! actual temperature, and the adiabatic liquid water content.
14
15	guez	196	! Equivalent de TLIFT entre NK1 et ICB1+1 inclus
16	guez	47
17	guez	195	! Differences with convect4:
18	guez	196	! - icbs1 is the first level above LCL (may differ from icb1)
19			! - in the iterations, used x(icbs1) instead x(icb1)
20			! - tvp1 is computed in only one time
21			! - icbs1: first level above Plcl1 (IMIN de TLIFT) in output
22			! - if icbs1=icb1, compute also tp1(icb1+1), tvp1(icb1+1) & clw1(icb1+1)
23	guez	47
24	guez	195	use cv30_param_m, only: minorig, nl
25			use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv
26			USE dimphy, ONLY: klev, klon
27	guez	47
28	guez	195	! inputs:
29	guez	196	integer, intent(in):: nk1(klon), icb1(klon)
30			real, intent(in):: t1(klon, klev)
31			real, intent(in):: q1(klon, klev), qs1(klon, klev), gz1(klon, klev)
32			real, intent(in):: p1(klon, klev)
33			real, intent(in):: plcl1(klon)
34	guez	47
35	guez	195	! outputs:
36	guez	196	real tp1(klon, klev), tvp1(klon, klev), clw1(klon, klev)
37	guez	47
38	guez	195	! local variables:
39			integer i, k
40	guez	196	integer icbs1(klon), icbsmax2
41	guez	195	real tg, qg, alv, s, ahg, tc, denom, es
42			real ah0(klon), cpp(klon)
43			real tnk(klon), qnk(klon), gznk(klon), ticb(klon), gzicb(klon)
44	guez	196	real qsicb(klon)
45			real cpinv(klon)
46	guez	47
47	guez	195	!-------------------------------------------------------------------
48	guez	47
49	guez	196	! Calculates the lifted parcel virtual temperature at nk1,
50	guez	195	! the actual temperature, and the adiabatic
51			! liquid water content. The procedure is to solve the equation.
52	guez	196	! cptp1+Lqp+phi=cptnk+Lqnk+gznk.
53	guez	47
54	guez	195	do i=1, klon
55	guez	196	tnk(i)=t1(i, nk1(i))
56			qnk(i)=q1(i, nk1(i))
57			gznk(i)=gz1(i, nk1(i))
58	guez	195	end do
59	guez	47
60	guez	195	! * Calculate certain parcel quantities, including static energy *
61	guez	47
62	guez	195	do i=1, klon
63			ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i) &
64			+qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
65			cpp(i)=cpd(1.-qnk(i))+qnk(i)cpv
66			cpinv(i)=1./cpp(i)
67			end do
68	guez	47
69	guez	195	! * Calculate lifted parcel quantities below cloud base *
70	guez	47
71	guez	196	do i=1, klon
72			! if icb1 is below LCL, start loop at ICB1+1:
73			! (icbs1 est le premier niveau au-dessus du LCL)
74			icbs1(i)=MIN(max(icb1(i), 2), nl)
75			if (plcl1(i) < p1(i, icbs1(i))) then
76			icbs1(i)=MIN(icbs1(i)+1, nl)
77	guez	195	endif
78	guez	196	enddo
79	guez	195
80	guez	196	do i=1, klon
81			ticb(i)=t1(i, icbs1(i))
82			gzicb(i)=gz1(i, icbs1(i))
83			qsicb(i)=qs1(i, icbs1(i))
84			enddo
85	guez	195
86	guez	196	! Re-compute icbsmax (icbsmax2):
87			icbsmax2=2
88			do i=1, klon
89			icbsmax2=max(icbsmax2, icbs1(i))
90	guez	195	end do
91
92			! initialization outputs:
93
94	guez	196	do k=1, icbsmax2
95			do i=1, klon
96			tp1(i, k) = 0.0
97			tvp1(i, k) = 0.0
98			clw1(i, k) = 0.0
99			enddo
100			enddo
101	guez	195
102	guez	196	! tp1 and tvp1 below cloud base:
103	guez	195
104			do k=minorig, icbsmax2-1
105			do i=1, klon
106	guez	196	tp1(i, k)=tnk(i)-(gz1(i, k)-gznk(i))*cpinv(i)
107			tvp1(i, k)=tp1(i, k)*(1.+qnk(i)/eps-qnk(i))
108	guez	195	end do
109			end do
110
111			! * Find lifted parcel quantities above cloud base *
112
113			do i=1, klon
114			tg=ticb(i)
115	guez	196	qg=qsicb(i)
116	guez	195	!debug alv=lv0-clmcpv*(ticb(i)-t0)
117			alv=lv0-clmcpv*(ticb(i)-273.15)
118
119			! First iteration.
120
121			s=cpd(1.-qnk(i))+clqnk(i) &
122	guez	196	+alvalvqg/(rrvticb(i)ticb(i))
123	guez	195	s=1./s
124	guez	47
125	guez	196	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i)
126	guez	195	tg=tg+s*(ah0(i)-ahg)
127	guez	47
128	guez	195	!debug tc=tg-t0
129			tc=tg-273.15
130			denom=243.5+tc
131	guez	196	denom=MAX(denom, 1.0)
132	guez	47
133	guez	195	es=6.112exp(17.67tc/denom)
134	guez	196	qg=epses/(p1(i, icbs1(i))-es(1.-eps))
135	guez	47
136	guez	195	! Second iteration.
137	guez	47
138	guez	196	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i)
139	guez	195	tg=tg+s*(ah0(i)-ahg)
140	guez	47
141	guez	195	!debug tc=tg-t0
142			tc=tg-273.15
143			denom=243.5+tc
144	guez	196	denom=MAX(denom, 1.0)
145	guez	47
146	guez	195	es=6.112exp(17.67tc/denom)
147	guez	47
148	guez	196	qg=epses/(p1(i, icbs1(i))-es(1.-eps))
149	guez	47
150	guez	195	alv=lv0-clmcpv*(ticb(i)-273.15)
151	guez	47
152	guez	196	! no approximation:
153			tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) &
154	guez	195	/(cpd+(cl-cpd)*qnk(i))
155
156	guez	196	clw1(i, icbs1(i))=qnk(i)-qg
157			clw1(i, icbs1(i))=max(0.0, clw1(i, icbs1(i)))
158	guez	195
159	guez	196	! (qg utilise au lieu du vrai mixing ratio rg)
160			tvp1(i, icbs1(i))=tp1(i, icbs1(i))*(1.+qg/eps-qnk(i))
161	guez	195
162			end do
163
164	guez	196	! * icbs1 is the first level above the LCL:
165			! if plcl1<p1(icb1), then icbs1=icb1+1
166			! if plcl1>p1(icb1), then icbs1=icb1
167	guez	195
168	guez	196	! * the routine above computes tvp1 from minorig to icbs1 (included).
169	guez	195
170	guez	196	! * to compute buoybase (in cv30_trigger.F), both tvp1(icb1) and
171			! tvp1(icb1+1) must be known. This is the case if icbs1=icb1+1,
172			! but not if icbs1=icb1.
173	guez	195
174	guez	196	! * therefore, in the case icbs1=icb1, we compute tvp1 at level icb1+1
175			! (tvp1 at other levels will be computed in cv30_undilute2.F)
176	guez	195
177			do i=1, klon
178	guez	196	ticb(i)=t1(i, icb1(i)+1)
179			gzicb(i)=gz1(i, icb1(i)+1)
180			qsicb(i)=qs1(i, icb1(i)+1)
181	guez	195	enddo
182
183			do i=1, klon
184			tg=ticb(i)
185	guez	196	qg=qsicb(i)
186	guez	195	!debug alv=lv0-clmcpv*(ticb(i)-t0)
187			alv=lv0-clmcpv*(ticb(i)-273.15)
188
189			! First iteration.
190
191			s=cpd(1.-qnk(i))+clqnk(i) &
192	guez	196	+alvalvqg/(rrvticb(i)ticb(i))
193	guez	195	s=1./s
194	guez	47
195	guez	196	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i)
196	guez	195	tg=tg+s*(ah0(i)-ahg)
197	guez	47
198	guez	195	!debug tc=tg-t0
199			tc=tg-273.15
200			denom=243.5+tc
201	guez	196	denom=MAX(denom, 1.0)
202	guez	47
203	guez	195	es=6.112exp(17.67tc/denom)
204	guez	47
205	guez	196	qg=epses/(p1(i, icb1(i)+1)-es(1.-eps))
206	guez	47
207	guez	195	! Second iteration.
208	guez	47
209	guez	196	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i)
210	guez	195	tg=tg+s*(ah0(i)-ahg)
211	guez	47
212	guez	195	!debug tc=tg-t0
213			tc=tg-273.15
214			denom=243.5+tc
215	guez	196	denom=MAX(denom, 1.0)
216	guez	47
217	guez	195	es=6.112exp(17.67tc/denom)
218
219	guez	196	qg=epses/(p1(i, icb1(i)+1)-es(1.-eps))
220	guez	195
221			alv=lv0-clmcpv*(ticb(i)-273.15)
222
223	guez	196	! no approximation:
224			tp1(i, icb1(i)+1)=(ah0(i)-gz1(i, icb1(i)+1)-alv*qg) &
225	guez	195	/(cpd+(cl-cpd)*qnk(i))
226
227	guez	196	clw1(i, icb1(i)+1)=qnk(i)-qg
228			clw1(i, icb1(i)+1)=max(0.0, clw1(i, icb1(i)+1))
229	guez	195
230	guez	196	! (qg utilise au lieu du vrai mixing ratio rg)
231			tvp1(i, icb1(i)+1)=tp1(i, icb1(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
232	guez	195	end do
233
234			end SUBROUTINE cv30_undilute1
235
236			end module cv30_undilute1_m