--- trunk/Sources/phylmd/CV3_routines/cv3_undilute1.f 2016/03/11 18:47:26 178 +++ trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f 2016/05/23 13:50:39 196 @@ -1,282 +1,236 @@ +module cv30_undilute1_m - SUBROUTINE cv3_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb & - ,tp,tvp,clw,icbs) - use cv3_param_m - use cvthermo - implicit none - -!---------------------------------------------------------------- -! 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 plclp(icb), then icbs=icb -! -! * the routine above computes tvp from minorig to icbs (included). -! -! * to compute buoybase (in cv3_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 cv3_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)) + implicit none -! 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 +contains - 460 continue + SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, nk1, icb1, tp1, & + tvp1, clw1, icbs1) - return - end + ! 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 plcl1p1(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