CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE GRDSET(ALATM,ALONM,IM,JM,DIM,DJM,L,IMAXG,JMAXG,IMAXM,          
     1  JMAXM,YLATM,FLONM)                                                      
C     ****                                                                      
C     ****     PROCEDURE TO SET UP BILINEAR TRANSFORMATION FOR                  
C     ****     TRANSFORMING A SCALAR FIELD FROM GEOMAGNETIC TO                  
C     ****     GEOGRAPHIC COORDINATES                                           
C     ****                                                                      
C     ****     GEOGRAPHIC GRID:                                                 
C     ****       IMAXG EQUALLY SPACED LONGITUDE POINTS.  RANGE                  
C     ****         -PI TO PI.                                                   
C     ****       JMAXG LATITUDE POINTS, NOT NECESSARILY EQUALLY                 
C     ****         SPACED.  RANGE -PI/2 TO PI/2.                                
C     ****                                                                      
C     ****     GEOMAGNETIC GRID:                                                
C     ****       IMAXM EQUALLY SPACED LONGITUDE POINTS.  RANGE                  
C     ****         FLONM (IN RADIANS) TO FLONM+2*PI                             
C     ****       JMAXM LATITUDE POINTS, NOT NECESSARILY EQUALLY                 
C     ****         SPACED.  RANGE -PI/2 TO PI/2.  MUST INCLUDE                  
C     ****         POLES                                                        
C     ****                                                                      
C     ****     INPUT:                                                           
C     ****                                                                      
C     ****       ALATM(L,JMAXG) IS 2-DIM ARRAY GIVING THE MAGNETIC
C     ****         LATITUDE (IN RADIANS) CORRESPONDING TO EACH OF               
C     ****         THE IMAXG*JMAXG GEOGRAPHIC GRID POINTS.                      
C     ****         (L.GE.IMAXG)                                                 
C     ****                                                                      
C     ****       ALONM(L,JMAXG) IS 2-DIM ARRAY GIVING THE MAGNETIC
C     ****         LONGITUDE (IN RADIANS) CORRESPONDING TO EACH OF              
C     ****         THE IMAXG*JMAXG GEOGRAPHIC GRID POINTS.                      
C     ****         (L.GE.IMAXG)                                                 
C     ****                                                                      
C     ****       IMAXG, JMAXG  DIMENSIONS OF GEOGRAPHIC GRID                    
C     ****                                                                      
C     ****       IMAXM, JMAXM  DIMENSIONS OF GEOMAGNETIC GRID                   
C     ****                                                                      
C     ****       YLATM(JMAXM) IS 1-DIM ARRAY GIVING THE JMAXM LATITUDE
C     ****         VALUES (IN RADIANS AND ARRANGED MONOTONICALLY)               
C     ****                                                                      
C     ****       L IS FIRST DIMENSION OF ARRAYS ALATM, ALONM, IM,               
C     ****         JM, DIM, DJM IN CALLING PROGRAM                              
C     ****                                                                      
C     ****       FLONM IS THE FIRST MAGNETIC LONGITUDE IN RADIANS               
C     ****                                                                      
C     ****     OUTPUT:                                                          
C     ****                                                                      
C     ****       IM(L,JMAXG) GIVES ROUNDED DOWN GEOMAGNETIC LONGITUDE
C     ****         INDEX FOR EACH GEOGRAPHIC GRID POINT                         
C     ****                                                                      
C     ****       JM(L,JMAXG) GIVES ROUNDED DOWN GEOMAGNETIC LATITUDE
C     ****         INDEX FOR EACH GEOGRAPHIC GRID POINT                         
C     ****                                                                      
C     ****       DIM(L,JMAXG) IS FRACTIONAL GEOMAGNETIC LONGITUDE INDEX
C     ****         FOR EACH GEOGRAPHIC GRID POINT.  USED FOR BILINEAR           
C     ****         INTERPOLATION WITHIN GRID CELL.
C     ****                                                                      
C     ****       DJM(L,JMAXG) IS FRACTIONAL GEOMAGNETIC LATITUDE INDEX
C     ****         FOR EACH GEOGRAPHIC GRID POINT.  USED FOR BILINEAR           
C     ****         INTERPOLATION WITHIN GRID CELL.
C     ****                                                                      
      DIMENSION ALATM(L,JMAXG),ALONM(L,JMAXG),IM(L,JMAXG),JM(L,JMAXG),
     | DIM(L,JMAXG),DJM(L,JMAXG),YLATM(JMAXM)
      DATA E/1.E-8/                                                             
      PI = 4.*ATAN(1.)                                                          
c     rtd = 180./pi
c     write (6,"(1x,'grdset: ylatm(rad)='/(1x,10e12.4))") ylatm
c     write (6,"(1x,'grdset: alatm(rad) for j=34'/(1x,10e12.4))")
c    |  (alatm(i,34),i=1,l)
c     write (6,"(1x,'grdset: alonm(rad) for j=5'/(1x,10e12.4))")
c    |  (alonm(i,34),i=1,l)
C     ****                                                                      
C     ****     CALCULATE JM FOR EACH GEOGRAPHIC GRID POINT                      
C     ****                                                                      
      DO 1 JG = 1,JMAXG                                                         
        DO 2 IG = 1,IMAXG                                                       
          JM(IG,JG) = 0.                                                        
    2   CONTINUE                                                                
        DO 3 J = 1,JMAXM-1                                                      
          DO 3 IG = 1,IMAXG                                                     
C  CVMGP=(n1,n2,exp) means n1 if exp >= 0, n2 if exp < 0
c           JM(IG,JG) = JM(IG,JG)+CVMGP(J,0,(ALATM(IG,JG)-YLATM(J))*
c    1      (YLATM(J+1)-ALATM(IG,JG)))
	    exp = (alatm(ig,jg)-ylatm(j))*(ylatm(j+1)-alatm(ig,jg))
	    iadd = 0
	    if (exp .ge. 0.) iadd = j
C  The following statement meant that if alatm=ylatm, 2 numbers are
C   added to jm, so it would be wrong and maybe out of range.
c           jm(ig,jg) = jm(ig,jg) + iadd
	    jm(ig,jg) = max0(jm(ig,jg),iadd)
c       if (jg .eq. 34 .and. ig .eq. 48) then
c       write (6,"(1x,'j jg ig iadd ylatm alatm exp=',4i3,5e12.4)")
c    |    j,jg,ig,iadd,ylatm(j),ylatm(j)*rtd,alatm(ig,jg),
c    |    alatm(ig,jg)*rtd,exp
c       endif
    3   CONTINUE                                                                
        JM(IMAXG+1,JG) = JM(1,JG)                                               
c       write (6,"(1x,'jm at jg=',i2/(1x,10i3))") jg,(jm(i,jg),i=1,imaxg)
    1 CONTINUE                                                                  
C     ****                                                                      
C     ****     EVALUATE IM, DIM, DJM                                            
C     ****                                                                      
      DO 5 JG = 1,JMAXG                                                         
       DO 4 IG = 1,IMAXG                                                        
        DIM(IG,JG) = (ALONM(IG,JG)-FLONM)*FLOAT(IMAXM)/(2.*PI)+IMAXM            
        IM(IG,JG) = IFIX(AMOD(DIM(IG,JG)+E,FLOAT(IMAXM))) + 1                   
        DIM(IG,JG) = DIM(IG,JG)-IFIX(DIM(IG,JG))                                
        DJM(IG,JG) = (ALATM(IG,JG)-YLATM(JM(IG,JG)))/(YLATM(JM(IG,JG)           
     1   +1)-YLATM(JM(IG,JG)))                                                  
    4  CONTINUE                                                                 
       IM(IMAXG+1,JG) = IM(1,JG)                                                
c       write (6,"(1x,'im dim djm at jg=',i2/(1x,10i3))") jg,(im(i,jg),
c    |   i=1,imaxg)
c       write (6,"(1x,10f8.2)") (dim(i,jg),i=1,imaxg)
c       write (6,"(1x,10f8.2)") (djm(i,jg),i=1,imaxg)
    5 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC        
