      subroutine distance(glat1,glon1,glat2,glon2,hkm,dkm)
              PI=ATAN(1.0)*4.
        UMR=PI/180.
	ss=sin(umr*glat1)*sin(umr*glat2)
	cc=cos(umr*glat1)*cos(umr*glat2)*cos(umr*(glon1-glon2))
	cres=ss+cc
	dkm=acos(cres)
	 dkm=dkm*(6370.0+hkm)
	return
	end
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.......................................................Sep. 2019
      subroutine subarea(glat1,glat2,dglon,skm2)
C     Area, km, for cell 
C
      Re=6371.0  ! Re, km
        PI=ATAN(1.0)*4.
        UMR=PI/180.
         rlat1=umr*glat1
	rlat2=umr*glat2
	deltaphi=abs(rlat2-rlat1)
	deltalon=abs(dglon*umr)
	skm2=0.5*Re**2*(deltaphi*(cos(rlat1)+cos(rlat2))*deltalon)
	RETURN
	ENd
