#include "dims.h"
      SUBROUTINE COLATH
      implicit none
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"
!     integer,PARAMETER :: IMAXG=ZIMX,JMAXG=ZJMX,IMAXGP=IMAXG+1,
!    |                     JMAXGP=JMAXG+1
#include "consdyn.h"
      real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc,
     |  rlatmp,rlonmp,dlons
      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)
#include "ioncr.h"
!
! Local:
      integer :: i,j
      real :: ofdc,cosofc,sinofc,aslonc,sinlat,coslat
C     ****
C     ****     PICK RELEVANT NORTHERN PARAMETERS FROM /IONCR/
C     ****
      OFDC = SQRT(OFFC(2)**2+DSKOFC(2)**2)
      COSOFC = COS(OFDC)
      SINOFC = SIN(OFDC)
      ASLONC = ASIN(DSKOFC(2)/OFDC)
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.)
          if (p(i,j) < 0.) p(i,j) = 0.
!         P(I,J) = merge(1.,P(I,J),P(I,J)-1.>=0.)
          if (p(i,j) >= 1.) p(i,j) = 1.          
    2 CONTINUE
      RETURN
      END
C