phylmd/CV30_routines/cv30_undilute1.f

module cv30_undilute1_m

  implicit none

contains

  SUBROUTINE cv30_undilute1(t, q, qs, gz, plcl, p, nk, icb, tp, tvp, clw, icbs)

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

    ! Equivalent de TLIFT entre NK et ICB+1 inclus

    ! Differences with convect4:
    ! - specify plcl in input
    ! - icbs is the first level above LCL (may differ from icb)
    ! - in the iterations, used x(icbs) instead x(icb)
    ! - many minor differences in the iterations
    ! - tvp is computed in only one time
    ! - icbs: first level above Plcl (IMIN de TLIFT) in output
    ! - if icbs=icb, compute also tp(icb+1), tvp(icb+1) & clw(icb+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):: nk(klon), icb(klon)
    real, intent(in):: t(klon, klev)
    real, intent(in):: q(klon, klev), qs(klon, klev), gz(klon, klev)
    real, intent(in):: p(klon, klev)
    real, intent(in):: plcl(klon) ! convect3

    ! outputs:
    real tp(klon, klev), tvp(klon, klev), clw(klon, klev)

    ! local variables:
    integer i, k
    integer icb1(klon), icbs(klon), icbsmax2 ! convect3
    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) ! convect3
    real cpinv(klon) ! convect3

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

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

    do i=1, klon
       tnk(i)=t(i, nk(i))
       qnk(i)=q(i, nk(i))
       gznk(i)=gz(i, nk(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 !convect3
       icb1(i)=MAX(icb(i), 2) !convect3
       icb1(i)=MIN(icb(i), nl) !convect3
       ! if icb is below LCL, start loop at ICB+1:
       ! (icbs est le premier niveau au-dessus du LCL)
       icbs(i)=icb1(i) !convect3
       if (plcl(i) < p(i, icb1(i))) then
          icbs(i)=MIN(icbs(i)+1, nl) !convect3
       endif
    enddo !convect3

    do i=1, klon !convect3
       ticb(i)=t(i, icbs(i)) !convect3
       gzicb(i)=gz(i, icbs(i)) !convect3
       qsicb(i)=qs(i, icbs(i)) !convect3
    enddo !convect3

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

    ! initialization outputs:

    do k=1, icbsmax2 ! convect3
       do i=1, klon ! convect3
          tp(i, k) = 0.0 ! convect3
          tvp(i, k) = 0.0 ! convect3
          clw(i, k) = 0.0 ! convect3
       enddo ! convect3
    enddo ! convect3

    ! tp and tvp below cloud base:

    do k=minorig, icbsmax2-1
       do i=1, klon
          tp(i, k)=tnk(i)-(gz(i, k)-gznk(i))*cpinv(i)
          tvp(i, k)=tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
       end do
    end do

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

    do i=1, klon
       tg=ticb(i)
       qg=qsicb(i) ! convect3
       !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)) ! convect3
       s=1./s

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

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

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

       ! Second iteration.

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

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

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

       qg=eps*es/(p(i, icbs(i))-es*(1.-eps))

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

       ! convect3: no approximation:
       tp(i, icbs(i))=(ah0(i)-gz(i, icbs(i))-alv*qg) &
            /(cpd+(cl-cpd)*qnk(i))

       clw(i, icbs(i))=qnk(i)-qg
       clw(i, icbs(i))=max(0.0, clw(i, icbs(i)))

       ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
       tvp(i, icbs(i))=tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing

    end do

    ! The following is only for convect3:

    ! * icbs is the first level above the LCL:
    ! if plcl<p(icb), then icbs=icb+1
    ! if plcl>p(icb), then icbs=icb

    ! * the routine above computes tvp from minorig to icbs (included).

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

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

    do i=1, klon
       ticb(i)=t(i, icb(i)+1)
       gzicb(i)=gz(i, icb(i)+1)
       qsicb(i)=qs(i, icb(i)+1)
    enddo

    do i=1, klon
       tg=ticb(i)
       qg=qsicb(i) ! convect3
       !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)) ! convect3
       s=1./s

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

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

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

       qg=eps*es/(p(i, icb(i)+1)-es*(1.-eps))

       ! Second iteration.

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

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

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

       qg=eps*es/(p(i, icb(i)+1)-es*(1.-eps))

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

       ! convect3: no approximation:
       tp(i, icb(i)+1)=(ah0(i)-gz(i, icb(i)+1)-alv*qg) &
            /(cpd+(cl-cpd)*qnk(i))

       clw(i, icb(i)+1)=qnk(i)-qg
       clw(i, icb(i)+1)=max(0.0, clw(i, icb(i)+1))

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