fsaut(:ncol,:) = 0._r8   fracw(:ncol,:) = 0._r8   fsacw(:ncol,:) = 0._r8   fsaci(:ncol,:) = 0._r8!! find the wet bulb temp and saturation value! for the provisional t and q without condensation!   call findsp (lchnk, ncol, qn, tn, p, tsp, qsp)   do 800 k = k1mb,pver      call vqsatd (t(1,k), p(1,k), es, qs, gamma, ncol)      do i = 1,ncol         relhum(i) = q(i,k)/qs(i)!         cldm(i) = max(cldn(i,k),mincld)!! the max cloud fraction above this level!         cldmax(i) = max(cldmax(i), cldm(i))! define the coefficients for C - E calculation         calpha(i) = 1.0/qs(i)         cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl         cbetah(i) = 1.0/qs(i)*gamma(i)*cpohl         cgamma(i) = calpha(i)+hlatv*cbeta(i)/cp         cgamah(i) = calpha(i)+hlatv*cbetah(i)/cp         rcgama(i) = cgamma(i)/cgamah(i)         if(cldm(i) > mincld) then            icwc(i) = max(0._r8,cwat(i,k)/cldm(i))         else            icwc(i) = 0.0         endif!! initial guess of evaporation, will be updated within iteration!         evapr(i,k) = conke*(1. - cldm(i))*sqrt(precab(i)) &                        *(1. - min(relhum(i),1._r8))!! zero cmeres before iteration for each level!         cmeres(i)=0.0      end do      do i = 1,ncol!! fractions of ice at this level!         tc = t(i,k) - t0         fice(i) = max(0._r8,min(-tc*0.05,1.0_r8))!! calculate the cooling due to a phase change of the rainwater! from above!         if (t(i,k) >= t0) then!              rmelt(i,k) =  -hlatf/cp*iceab(i)*gravit/pdel(i,k)            rmelt(i,k) = 0.            iceab(i) = 0.         else            rmelt(i,k) = 0.         endif      end do!! calculate cme and formation of precip. !! The cloud microphysics is highly nonlinear and coupled with cme! Both rain processes and cme are calculated iteratively.!       do 100 l = 1,iter        do i = 1,ncol!! calculation of cme has 4 scenarios! ==================================!           if(relhum(i) > rhu00(i,k)) then               ! 1. whole grid saturation           ! ========================              if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then                 cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)           ! 2. fractional saturation           ! ========================              else                  csigma(i) = 1.0/(rhdfda(i,k)+cgamma(i)*icwc(i))                  cmec1(i) = (1.0-cldm(i))*csigma(i)*rhdfda(i,k)                  cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0-rcgama(i)*cldm(i))*   &                             csigma(i)*calpha(i)*icwc(i)                  cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) +  &                           (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)                  cmec4(i) = csigma(i)*cgamma(i)*icwc(i)                  ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er                    cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k)  &                             -cmec3(i)*ttend(i,k) + cmec4(i)*evapr(i,k)               endif           ! 3. when rh < rhu00, evaporate existing cloud water           ! ==================================================            else if(cwat(i,k) > 0.0)then              ! liquid water should be evaporated but not to exceed               ! saturation point. if qn > qsp, not to evaporate cwat              cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat            ! 4. no condensation nor evaporation           ! ==================================           else              cme(i,k)=0.0           endif          end do    !end loop for cme update! Because of the finite time step, ! place a bound here not to exceed wet bulb point! and not to evaporate more than available water!         do i = 1, ncol            qtmp = qn(i,k) - cme(i,k)*deltat! possibilities to have qtmp > qsp!!   1. if qn > qs(tn), it condenses; !      if after applying cme,  qtmp > qsp,  more condensation is applied. !      !   2. if qn < qs, evaporation should not exceed qsp,                if(qtmp > qsp(i,k)) then              cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat            endif!! if net evaporation, it should not exceed available cwat!            if(cme(i,k) < -cwat(i,k)/deltat)  &               cme(i,k) = -cwat(i,k)/deltat!! addition of residual condensation from previous step of iteration!            cme(i,k) = cme(i,k) + cmeres(i)         end do         do i = 1,ncol!! as a safe limit, condensation should not reduce grid mean rh below rhu00!            if(cme(i,k) > 0.0 .and. relhum(i) > rhu00(i,k) )  &              cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu00(i,k))/deltat)!! initial guess for cwm (mean cloud water over time step) if 1st iteration!           if(l < 2) then             cwm(i) = max(cwat(i,k)+cme(i,k)*dto2,  0._r8)           endif         enddo! provisional precipitation falling through model layer         do i = 1,ncol            prprov(i) = prect(i) + prain(i,k)*pdel(i,k)/gravit         end do! calculate conversion of condensate to precipitation by cloud microphysics          call findmcnew (lchnk   ,ncol    , &                         k       ,prprov  ,t       ,p       , &                         cwm     ,cldm    ,cldmax  ,fice    ,coef    , &                         fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &                         icefrac)!! calculate the precip rate!         do i = 1,ncol            if (cldm(i) > 0) then  !! first predict the cloud water!               cdt = coef(i)*deltat               if (cdt > 0.01) then                  pol = cme(i,k)/coef(i) ! production over loss                  cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)               else                  cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))               endif!! now back out the tendency of net rain production!               prain(i,k) =  max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)            else               prain(i,k) = 0.0               cwn(i) = 0.            endif!! update any remaining  provisional values!            cwm(i) = (cwn(i) + cwat(i,k))*0.5!! update in cloud water!            if(cldm(i) > mincld) then               icwc(i) = cwm(i)/cldm(i)            else               icwc(i) = 0.0            endif         end do              ! end of do i = 1,ncol!! calculate provisional value of cloud water for! evaporation of rain (evapr) calculation!      do i = 1,ncol         qtmp = qn(i,k) - cme(i,k)*deltat         ttmp = tn(i,k) + rmelt(i,k)*deltat + hlocp*deltat*cme(i,k)         esn = estblf(ttmp)         qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)         qtl(i) = max((qsn - qtmp)/deltat,0._r8)         relhum1(i) = qtmp/qsn      end do!      do i = 1,ncol#ifdef PERGRO         evapr(i,k) = conke*(1. - max(cldm(i),mincld))* &                      sqrt(precab(i))*(1. - min(relhum1(i),1._r8))#else         evapr(i,k) = conke*(1. - cldm(i))*sqrt(precab(i)) &                      *(1. - min(relhum1(i),1._r8))#endif!! limit the evaporation to the amount which is entering the box! or saturates the box!         prtmp = precab(i)*gravit/pdel(i,k)         evapr(i,k) = min(evapr(i,k), prtmp, qtl(i))*omsm#ifdef PERGRO!           zeroing needed for pert growth         evapr(i,k) = 0.#endif      end do! now remove the residual of any over-saturation. Normally,! the oversaturated water vapor should have been removed by ! cme formulation plus constraints by wet bulb tsp/qsp! as computed above. However, because of non-linearity,! addition of (cme-evapr) to update t and q may still cause! a very small amount of over saturation. It is called a! residual of over-saturation because theoretically, cme! should have taken care of all of large scale condensation.!        do i = 1,ncol          qtmp = qn(i,k)-(cme(i,k)-evapr(i,k))*deltat          ttmp = tn(i,k)+(rmelt(i,k)+hlocp*(cme(i,k)-evapr(i,k)) )  &                      *deltat          esn = estblf(ttmp)          qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)          !          !Upper stratosphere and mesosphere, qsn calculated          !above may be negative. Here just to skip it instead          !of resetting it to 1 as in aqsat          !          if(qtmp > qsn .and. qsn > 0) then             !calculate dqsdt, a more precise calculation             !which taking into account different range of T              !can be found in aqsatd.F. Here follows             !cond.F to calculate it.             !             denom = (p(i,k)-omeps*esn)*ttmp*ttmp             dqsdt = clrh2o*qsn*p(i,k)/denom             !             !now extra condensation to bring air to just saturation             !             ctmp = (qtmp-qsn)/(1.+hlocp*dqsdt)/deltat             cme(i,k) = cme(i,k)+ctmp!! save residual on cmeres to addtion to cme on entering next iteration! cme exit here contain the residual but overrided if back to iteration!             cmeres(i) = ctmp          else             cmeres(i) = 0.0          endif       end do               100 continue              ! end of do l = 1,iter!! precipitation!      do i = 1,ncol         prtmp = pdel(i,k) / gravit *(prain(i,k) - evapr(i,k))         iceab(i) = iceab(i) + fice(i)*prtmp         precab(i) = precab(i) + prtmp         prect(i) = prect(i) + prtmp + pcflx(i,k+1)         if ((precab(i)) < 1.e-10) then                  precab(i) = 0.         endif         if ((prect(i)) < 1.e-10) then                  prect(i) = 0.         endif      end do 800 continue                ! level loop (k=1,pver)   returnend subroutine pcond!##############################################################################subroutine findmcnew (lchnk   ,ncol    , &                      k       ,precab  ,t       ,p       , &                      cwm     ,cldm    ,cldmax  ,fice    ,coef    , &                      fwaut   ,fsaut   ,fracw   ,fsacw   ,fsaci   , &                      icefrac )!----------------------------------------------------------------------- ! ! Purpose: ! calculate the conversion of condensate to precipitate! ! Method: ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3!  model climate using diagnosed and !  predicted condensate parameterizations, 1998, J. Clim., 11,!  pp1587---1614.! ! Author: P. Rasch! !-----------------------------------------------------------------------   use phys_grid, only: get_rlat_all_p