subroutine gg_gm
     a  (flag,pole_lat,pole_lon,gglat,gglon,gmlat,gmlon)
c
c	s/r to convert between geographic and geomagnetic coords (both ways)
c	can also be used to convert between any other coord systems
c	written by leo mcnamara (from the IRI s/r GGM) 30-jul-90
c-----	last modified 31-jul-90
c
	real lat
	integer flag
c	flag =  1 for geographic to geomagnetic conversion (gg -> gm)
c	flag = -1 for geomagnetic to geographic conversion
c	pole_lat is the geogc latitude of the north geom pole
c	pole_lon is the geogc longitude of the north geom pole
c
	data pi,radeg,degrad / 3.141592654,57.29577951,0.017453292 /
c
	cp = cos((90.0 - pole_lat) * degrad)
	sp = sin((90.0 - pole_lat) * degrad)
c
	if(flag.eq.1) slat = gglat
	if(flag.eq.-1) slat = gmlat
	if(flag.eq.1) slon = gglon +(360.0 - pole_lon)
	if(flag.eq.-1) slon = gmlon
c
	cb = cos(slat*degrad)
	sb = sin(slat*degrad)
	cl = cos(slon*degrad)
	sl = sin(slon*degrad)
	sbb = sb * cp + flag * cb * cl * sp
	if(abs(sbb).gt.1.0) sbb = sign(1.0,sbb)
	lat = asin(sbb)
c
	if(flag.eq.1) gmlat = lat * radeg
	if(flag.eq.-1) gglat = lat * radeg
c
	cg = cos(lat)
	cll = (-flag * sb * sp + cb * cl * cp) / cg
ccc	if(abs(cll).gt.1.) print*,' GG_GM --  cll ', cll
	if(abs(cll).gt.1.) cll = sign(1.0,cll)
	flon = acos(cll)
	sll = cb * sl / cg
	if(sll.lt.0.) flon = 2.0 * pi - flon
c
	if(flag.eq.1) gmlon = flon * radeg
	if(flag.eq.-1) gglon = flon * radeg - (360.0 - pole_lon)
c
	if(gmlon.gt.360.) gmlon = gmlon - 360.
	if(gmlon.le.0.  ) gmlon = gmlon + 360.
	if(gglon.gt.360.) gglon = gglon - 360.
	if(gglon.le.0.  ) gglon = gglon + 360.
c
	return
	end