anchor.diff:              if (qv<qvsat) then
                   if (qr>SMALL_VALUE) then
                       ! evaporate rain
  -                    call phase_change(p,t,qr,qvsat,qv,L_evap,1.0) !cloud2rain)
  +                    call phase_change(p,t,qr,qvsat,qv,L_evap,cloud2rain)
                   endif
                   if (qs>SMALL_VALUE) then
                       ! sublimate snow
  -                    call phase_change(p,t,qs,qvsat,qv,L_subl,1.0) !cloud2snow)
  +                    call phase_change(p,t,qs,qvsat,qv,L_subl,cloud2snow)
                   endif
               endif
           endif
  
positive.code:     real,parameter :: freezing_threshold=273.15          ! [K]
      real,parameter :: snow_fall_rate=1.5                 ! [m/s]   for a water vapor scale height of 3750m corresponds to tau_f = 2500
      real,parameter :: rain_fall_rate=10.0                ! [m/s]   for a water vapor scale height of 3750m corresponds to tau_f = 375
      real,parameter :: snow_cloud_init=0.00001            ! [kg/kg] cloud ice content before snow will start to form
      real,parameter :: rain_cloud_init=0.00001            ! [kg/kg] cloud water content before rain will start to form
      
      
  !   these are recalculated every call because they are a function of dt
  
negative.code:                 th(i,:,j)=t/pii(i,:,j)
                  
              enddo
  !             where(qs(:,:,j)<0) qs(:,:,j)=0
  !             where(qr(:,:,j)<0) qr(:,:,j)=0
  !             where(qc(:,:,j)<0) qc(:,:,j)=0
  !             where(qv(:,:,j)<0) qv(:,:,j)=0
          enddo
          !$omp end do
          deallocate(t)
  