      SUBROUTINE COLATH
C     ****
C     ****     FOR EACH NORTHERN HEMISPHERE GEOMAGNETIC GRID POINT,
C     ****       DETERMINE COLATITUDE IN NORTHERN AURORAL CIRCLE
C     ****       COORDINATES.  PLACE VALUES IN COLATC(IMAXMP,JMX0).
C     ****       ALSO CALCULATE FRACTIONAL PRESENCE OF DYNAMO
C     ****       EQUATION GIVEN CRITICAL COLATITUDES CRIT(2)
C     ****
      include "params.h"
      PARAMETER(IMAXG=ZIMX,JMAXG=ZJMX,IMAXGP=IMAXG+1,JMAXGP=JMAXG+1)
      include "consts.h"
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3  RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),
     1  DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),
     2  PHINM0(2),PSIM(2),PSIE(2),PCEN(2),RR1(2),COLATC(IMAXMP,JMX0),
     3  P(IMX0,JMX0),CRIT(2)
C     ****
C     ****     PICK RELEVANT NORTHERN PARAMETERS FROM /IONCR/
C     ****
c       write (6,"(1x,'colath: crit1,2=',2f7.2)") crit(1)*180./3.1416,
c    |   crit(2)*180./3.1416
      IF (ABS(OFFC(2)+DSKOFC(2)) .GT. 0.01) THEN
      OFDC = SQRT(OFFC(2)**2+DSKOFC(2)**2)
      COSOFC = COS(OFDC)
      SINOFC = SIN(OFDC)
      ASLONC = ASIN(DSKOFC(2)/OFDC)
      ELSE
	OFDC = 0.
	COSOFC = 1.
	SINOFC = 0.
	ASLONC = 0.
      ENDIF
C     ****
C     ****     FILL ARRAY COLATC WITH REQUIRED COLATITUDE VALUES
C     ****
      DO 1 J = 1,JMX0
        SINLAT = SIN(ABS(YLATM(J+JMX0-1)))
        COSLAT = COS(YLATM(J+JMX0-1))
        DO 1 I = 1,IMAXMP
          COLATC(I,J) = COS(YLONM(I)-DLONS(1)+ASLONC)
          COLATC(I,J) = COSOFC*SINLAT-SINOFC*COSLAT*COLATC(I,J)
          COLATC(I,J) = ACOS(COLATC(I,J))
    1 CONTINUE
C     ****
C     ****     CALCULATE FRACTIONAL PRESENCE OF DYNAMO EQUATION
C     ****       AT EACH NORTHERN HEMISPHERE GEOMAGNETIC GRID
C     ****       POINT.  PLACE VALUES IN P(IMX0,JMX0)
C     ****
      DO 2 J = 1,JMX0
        DO 2 I = 1,IMX0
          P(I,J) = (COLATC(I,J)-CRIT(1))/(CRIT(2)-CRIT(1))
          P(I,J) = merge(P(I,J),0.,P(I,J)>=0.)
          P(I,J) = merge(1.,P(I,J),P(I,J)-1.>=0.)
    2 CONTINUE
      RETURN
      END
C
