anchor.diff:  			excess=qv-qvsat
   ! 			temperature change if all vapor is converted
   			deltat=excess*vapor2temp
  -! 			Now calculate new saturated mixing ratio at the new hypothetical tempera
  +! 			Now calculate new saturated mixing ratio at the new hypothetical temperature
   			qvsat=sat_mr(t+deltat*0.5,p)
   			excess=qv-qvsat
   			t=t+(excess*vapor2temp)
   			qv=qv-excess
   			qc=qc+excess
   			
  -! 		if unsaturated anc clouds exist, evaporate clouds
  +! 		if unsaturated and clouds exist, evaporate clouds
   		else if (qc>0) then
   			excess=qvsat-qv
   			if (excess<qc) then
  
positive.code: 							rain(i,j),snow(i,j),&
  							dt,dz(i,:,j),nz,((i==(nx/2+20)).and.(j==2)))
  				th(i,:,j)=t/pii(i,:,j)
  ! 				print*, maxval(t)
  ! 				if (minval(qc(i,:,j))<0) then
  ! 					print*,i,j,minval(qc(i,:,j)), "cloud"
  ! 				endif
  ! 				if (minval(qv(i,:,j))<0) then
  ! 					print*,i,j,minval(qv(i,:,j)), "vapor"
  ! 				endif
  ! 				if (minval(qr(i,:,j))<0) then
  ! 					print*,i,j,minval(qr(i,:,j)), "rain"
  ! 				endif
  ! 				if (minval(qs(i,:,j))<0) then
  ! 					print*,i,j,minval(qs(i,:,j)), "snow"
  ! 				endif
  				
  			enddo
  		enddo
  
negative.code: ! 			Murray, F. W., 1967: On the computation of saturation vapor pressure. 
  ! 				Journal of Applied Meteorology, Vol. 6, pp. 203-204.
  ! 		Also notes a 6th order polynomial and look up table as viable options. 
  		if (t<freezing_threshold) then
  			a=21.8745584
  			b=7.66
  		else
  