--- trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f 2016/03/29 15:20:23 189 +++ trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f 2016/05/23 13:50:39 196 @@ -1,287 +1,236 @@ +module cv30_undilute1_m - SUBROUTINE cv30_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb & - ,tp,tvp,clw,icbs) - use cv30_param_m - use cvthermo - implicit none + 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 ICB for convect4, up through ICB + 1 for convect3) - ! Calculates the lifted parcel virtual temperature at nk, the + ! (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 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) -!---------------------------------------------------------------- - - -! inputs: - integer, intent(in):: len, nd - integer nk(len), icb(len) - real, intent(in):: t(len,nd) - real, intent(in):: q(len,nd), qs(len,nd), gz(len,nd) - real, intent(in):: p(len,nd) - real plcl(len) ! convect3 - -! outputs: - real tp(len,nd), tvp(len,nd), clw(len,nd) - -! local variables: - integer i, k - integer icb1(len), icbs(len), icbsmax2 ! convect3 - real tg, qg, alv, s, ahg, tc, denom, es - real ah0(len), cpp(len) - real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len) - real qsicb(len) ! convect3 - real cpinv(len) ! 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 320 i=1,len - tnk(i)=t(i,nk(i)) - qnk(i)=q(i,nk(i)) - gznk(i)=gz(i,nk(i)) -! ori ticb(i)=t(i,icb(i)) -! ori gzicb(i)=gz(i,icb(i)) - 320 continue -! -! *** Calculate certain parcel quantities, including static energy *** -! - do 330 i=1,len - 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) - 330 continue -! -! *** Calculate lifted parcel quantities below cloud base *** -! - do i=1,len !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).lt.p(i,icb1(i))) then - icbs(i)=MIN(icbs(i)+1,nl) !convect3 - endif - enddo !convect3 - - do i=1,len !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 310 i=1,len !convect3 - icbsmax2=max(icbsmax2,icbs(i)) !convect3 - 310 continue !convect3 - -! initialization outputs: - - do k=1,icbsmax2 ! convect3 - do i=1,len ! 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 350 k=minorig,icbsmax2-1 - do 340 i=1,len - 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) - 340 continue - 350 continue -! -! *** Find lifted parcel quantities above cloud base *** -! - do 360 i=1,len - tg=ticb(i) -! ori qg=qs(i,icb(i)) - qg=qsicb(i) ! convect3 -!debug alv=lv0-clmcpv*(ticb(i)-t0) - alv=lv0-clmcpv*(ticb(i)-273.15) -! -! First iteration. -! -! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) - s=cpd*(1.-qnk(i))+cl*qnk(i) & - +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 - s=1./s -! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) - ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 - tg=tg+s*(ah0(i)-ahg) -! ori tg=max(tg,35.0) -!debug tc=tg-t0 - tc=tg-273.15 - denom=243.5+tc - denom=MAX(denom,1.0) ! convect3 -! ori if(tc.ge.0.0)then - es=6.112*exp(17.67*tc/denom) -! ori else -! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) -! ori endif -! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) - qg=eps*es/(p(i,icbs(i))-es*(1.-eps)) -! -! Second iteration. -! - -! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) -! ori s=1./s -! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) - ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 - tg=tg+s*(ah0(i)-ahg) -! ori tg=max(tg,35.0) -!debug tc=tg-t0 - tc=tg-273.15 - denom=243.5+tc - denom=MAX(denom,1.0) ! convect3 -! ori if(tc.ge.0.0)then - es=6.112*exp(17.67*tc/denom) -! ori else -! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) -! ori end if -! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) - qg=eps*es/(p(i,icbs(i))-es*(1.-eps)) - - alv=lv0-clmcpv*(ticb(i)-273.15) - -! ori c approximation here: -! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) -! ori & -gz(i,icb(i))-alv*qg)/cpd - -! convect3: no approximation: - tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg) & - /(cpd+(cl-cpd)*qnk(i)) - -! ori clw(i,icb(i))=qnk(i)-qg -! ori clw(i,icb(i))=max(0.0,clw(i,icb(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 - - 360 continue -! -! ori do 380 k=minorig,icbsmax2 -! ori do 370 i=1,len -! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i) -! ori 370 continue -! ori 380 continue -! - -! -- The following is only for convect3: -! -! * icbs is the first level above the LCL: -! 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,len
- ticb(i)=t(i,icb(i)+1)
- gzicb(i)=gz(i,icb(i)+1)
- qsicb(i)=qs(i,icb(i)+1)
- enddo
-
- do 460 i=1,len
- tg=ticb(i)
- qg=qsicb(i) ! convect3
-!debug alv=lv0-clmcpv*(ticb(i)-t0)
- alv=lv0-clmcpv*(ticb(i)-273.15)
-!
-! First iteration.
-!
-! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
- s=cpd*(1.-qnk(i))+cl*qnk(i) &
- +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
- s=1./s
-! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
- ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
- tg=tg+s*(ah0(i)-ahg)
-! ori tg=max(tg,35.0)
-!debug tc=tg-t0
- tc=tg-273.15
- denom=243.5+tc
- denom=MAX(denom,1.0) ! convect3
-! ori if(tc.ge.0.0)then
- es=6.112*exp(17.67*tc/denom)
-! ori else
-! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-! ori endif
-! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
- qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
-!
-! Second iteration.
-!
-
-! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
-! ori s=1./s
-! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
- ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
- tg=tg+s*(ah0(i)-ahg)
-! ori tg=max(tg,35.0)
-!debug tc=tg-t0
- tc=tg-273.15
- denom=243.5+tc
- denom=MAX(denom,1.0) ! convect3
-! ori if(tc.ge.0.0)then
- es=6.112*exp(17.67*tc/denom)
-! ori else
-! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-! ori end if
-! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
- qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
-
- alv=lv0-clmcpv*(ticb(i)-273.15)
-
-! ori c approximation here:
-! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
-! ori & -gz(i,icb(i))-alv*qg)/cpd
-
-! convect3: no approximation:
- tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg) &
- /(cpd+(cl-cpd)*qnk(i))
-
-! ori clw(i,icb(i))=qnk(i)-qg
-! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
- clw(i,icb(i)+1)=qnk(i)-qg
- clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1))
+ ! 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