phylmd/CV30_routines/cv30_undilute1.f

module cv30_undilute1_m

  implicit none

contains

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

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

    ! Equivalent de TLIFT entre MINORIG 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: clmcpv, eps
    USE dimphy, ONLY: klev, klon
    use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rv

    ! inputs:
    integer, intent(in):: 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 minorig,
    !  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, minorig)
       qnk(i)=q1(i, minorig)
       gznk(i)=gz1(i, minorig)
    end do

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

    do i=1, klon
       ah0(i)=(rcpd*(1.-qnk(i))+rcw*qnk(i))*tnk(i) &
            +qnk(i)*(rlvtt-clmcpv*(tnk(i)-273.15))+gznk(i)
       cpp(i)=rcpd*(1.-qnk(i))+qnk(i)*rcpv
       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=rlvtt-clmcpv*(ticb(i)-t0)
       alv=rlvtt-clmcpv*(ticb(i)-273.15)

       ! First iteration.

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

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

       ! no approximation:
       tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) &
            /(rcpd+(rcw-rcpd)*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=rlvtt-clmcpv*(ticb(i)-t0)
       alv=rlvtt-clmcpv*(ticb(i)-273.15)

       ! First iteration.

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

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

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