--- trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f 2016/05/18 17:56:44 195 +++ trunk/Sources/phylmd/CV30_routines/cv30_undilute1.f 2016/06/06 17:42:15 201 @@ -4,111 +4,108 @@ contains - SUBROUTINE cv30_undilute1(t, q, qs, gz, plcl, p, nk, icb, tp, tvp, clw, icbs) + SUBROUTINE cv30_undilute1(t1, q1, qs1, gz1, plcl1, p1, 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 minorig, the ! actual temperature, and the adiabatic liquid water content. - ! Equivalent de TLIFT entre NK et ICB+1 inclus + ! Equivalent de TLIFT entre MINORIG et ICB1+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) + ! - 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 cv_thermo_m, only: clmcpv, eps USE dimphy, ONLY: klev, klon + use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rv ! 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 + 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 tp(klon, klev), tvp(klon, klev), clw(klon, klev) + real tp1(klon, klev), tvp1(klon, klev), clw1(klon, klev) ! local variables: integer i, k - integer icb1(klon), icbs(klon), icbsmax2 ! convect3 + 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) ! convect3 - real cpinv(klon) ! convect3 + real qsicb(klon) + real cpinv(klon) !------------------------------------------------------------------- - ! Calculates the lifted parcel virtual temperature at nk, + ! 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*tp+L*qp+phi=cp*tnk+L*qnk+gznk. + ! cp*tp1+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)) + 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)=(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 + 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 !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 + 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 !convect3 + enddo + + do i=1, klon + ticb(i)=t1(i, icbs1(i)) + gzicb(i)=gz1(i, icbs1(i)) + qsicb(i)=qs1(i, icbs1(i)) + enddo - 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 + ! Re-compute icbsmax (icbsmax2): + icbsmax2=2 + do i=1, klon + icbsmax2=max(icbsmax2, icbs1(i)) 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 + 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 - ! tp and tvp below cloud base: + ! tp1 and tvp1 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) + 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 @@ -116,124 +113,123 @@ 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) + qg=qsicb(i) + !debug alv=rlvtt-clmcpv*(ticb(i)-t0) + alv=rlvtt-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=rcpd*(1.-qnk(i))+rcw*qnk(i) & + +alv*alv*qg/(rv*ticb(i)*ticb(i)) s=1./s - ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 + 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) ! convect3 + denom=MAX(denom, 1.0) es=6.112*exp(17.67*tc/denom) - qg=eps*es/(p(i, icbs(i))-es*(1.-eps)) + qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps)) ! Second iteration. - ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 + 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) ! convect3 + denom=MAX(denom, 1.0) es=6.112*exp(17.67*tc/denom) - qg=eps*es/(p(i, icbs(i))-es*(1.-eps)) + qg=eps*es/(p1(i, icbs1(i))-es*(1.-eps)) - alv=lv0-clmcpv*(ticb(i)-273.15) + alv=rlvtt-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)) + ! no approximation: + tp1(i, icbs1(i))=(ah0(i)-gz1(i, icbs1(i))-alv*qg) & + /(rcpd+(rcw-rcpd)*qnk(i)) - clw(i, icbs(i))=qnk(i)-qg - clw(i, icbs(i))=max(0.0, clw(i, icbs(i))) + clw1(i, icbs1(i))=qnk(i)-qg + clw1(i, icbs1(i))=max(0.0, clw1(i, icbs1(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 + ! (qg utilise au lieu du vrai mixing ratio rg) + tvp1(i, icbs1(i))=tp1(i, icbs1(i))*(1.+qg/eps-qnk(i)) end do - ! The following is only for convect3: - - ! * icbs is the first level above the LCL: - ! if plclp(icb), then icbs=icb + ! * icbs1 is the first level above the LCL: + ! if plcl1p1(icb1), then icbs1=icb1 - ! * the routine above computes tvp from minorig to icbs (included). + ! * the routine above computes tvp1 from minorig to icbs1 (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. + ! * 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 icbs=icb, we compute tvp at level icb+1 - ! (tvp at other levels will be computed in cv30_undilute2.F) + ! * 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)=t(i, icb(i)+1) - gzicb(i)=gz(i, icb(i)+1) - qsicb(i)=qs(i, icb(i)+1) + 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) ! convect3 - !debug alv=lv0-clmcpv*(ticb(i)-t0) - alv=lv0-clmcpv*(ticb(i)-273.15) + qg=qsicb(i) + !debug alv=rlvtt-clmcpv*(ticb(i)-t0) + alv=rlvtt-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=rcpd*(1.-qnk(i))+rcw*qnk(i) & + +alv*alv*qg/(rv*ticb(i)*ticb(i)) s=1./s - ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 + 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) ! convect3 + denom=MAX(denom, 1.0) es=6.112*exp(17.67*tc/denom) - qg=eps*es/(p(i, icb(i)+1)-es*(1.-eps)) + 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) ! convect3 + 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) ! convect3 + denom=MAX(denom, 1.0) es=6.112*exp(17.67*tc/denom) - qg=eps*es/(p(i, icb(i)+1)-es*(1.-eps)) + qg=eps*es/(p1(i, icb1(i)+1)-es*(1.-eps)) - alv=lv0-clmcpv*(ticb(i)-273.15) + alv=rlvtt-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)) + ! no approximation: + tp1(i, icb1(i)+1)=(ah0(i)-gz1(i, icb1(i)+1)-alv*qg) & + /(rcpd+(rcw-rcpd)*qnk(i)) - clw(i, icb(i)+1)=qnk(i)-qg - clw(i, icb(i)+1)=max(0.0, clw(i, icb(i)+1)) + clw1(i, icb1(i)+1)=qnk(i)-qg + clw1(i, icb1(i)+1)=max(0.0, clw1(i, icb1(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 + ! (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