       program inline2camx

!----------------------------------------------------------------------
!   CMAQ2ASC 1.0 Convert 3D CMAQ conc to gridded ascii file
!----------------------------------------------------------------------

c *** Variable declarations
      implicit none

c garnet:

c3 ifort -check -traceback -O2 -align dcommons -extend_source -convert big_endian -ipo -o inline2camx.x inline2camx.F -L/garnet/home/wdx/lib/x86_64/ifc/ioapi_3/Linux2_x86_64ifort -lioapi -L/usr/local/netcdf-4.0/lib -lnetcdf -lnetcdff
      
      include '/garnet/home/wdx/lib/x86_64/ifc/ioapi_3/ioapi/fixed_src/PARMS3.EXT'
      include '/garnet/home/wdx/lib/x86_64/ifc/ioapi_3/ioapi/fixed_src/IODECL3.EXT'
      include '/garnet/home/wdx/lib/x86_64/ifc/ioapi_3/ioapi/fixed_src/FDESC3.EXT'

c '

      INTEGER :: NCOLS, NROWS, NLAYS, NVARS, JDATE, JTIME, RUNLEN
      INTEGER :: TSTEP, NSTEPS, I,V, VMAX, UMAX, DMAX
      CHARACTER*16    VNAME( MXVARS3 ) !  list of vble names, from user
      CHARACTER*16    UNITS( MXVARS3 ) !  list of vble units
      CHARACTER*80    VDESC( MXVARS3 ) !  list of vble descs
      CHARACTER*256   MESG    !  buffer for m3exit(), etc
      CHARACTER*16    INAME   !  logical name of the input file
      LOGICAL         EFLAG !flag: error has happened

      REAL    :: rlon0,rlat0,tlat1,tlat2,xloc,yloc,rlat,rlon,xpos,ypos
      REAL    :: xorg,yorg,dx,rone,rzero,DUP(9)
      INTEGER :: xorig,yorig,nx,ny,IH,ione,izero,J,tempa,tempb,tempc,cdate
      integer    c, r, k, layer, hour, ihour, idate, n, nfiles,HOURS
      real cfac

      REAL, ALLOCATABLE :: TOTS( :,:,: ),EMIS(:,:),QA(:,:,:)
      REAL, ALLOCATABLE :: zero(:)
      REAL, ALLOCATABLE :: GRID( :,:,:),ALTGRID(:,:),xloca(:,:),yloca(:,:)
      REAL, ALLOCATABLE :: stkdm(:,:),stkht(:,:),stktk(:,:),stkve(:,:),stkflw(:,:)
      INTEGER, ALLOCATABLE :: row(:,:),col(:,:)

      CHARACTER*4 IFILE(10), NOTE(60), MSPEC(10,1000)
      CHARACTER*10 cspec

      CHARACTER(LEN=16) :: progname
      CHARACTER(LEN=16) :: met, outfile
      CHARACTER(LEN=180) :: ename(366),fname,land

      INTEGER :: istatus
      INTEGER :: TRIMLEN
      real       envreal

      EXTERNAL TRIMLEN, envreal

      data IFILE/'P','T','S','O','U','R','C','E',' ',' '/
      NOTE = "INLINE2CAMX VERSION 1.0" 
      eflag = .false.

c *** Initialize variables
      progname = 'INLINE2CAMX'
      cfac = (1./1000.) * 0.001102293 * 3600. !convert from g to tons and sec to hour
      ione = 1
      izero = 0
      rone = 1.
      rzero = 0.

      write(*,*) 'max vars',MXVARS3

      read(*,'(a)') fname
      read(fname,*) rlon0,rlat0,tlat1,tlat2,nx,ny,xorg,yorg,dx
      write(*,'(a,t20,3f10.0)') 'Projection:',rlon0,rlat0,tlat1,tlat2
      write(*,*) 'X and Y cells:',nx,ny
      xorg=xorg*1000.
      yorg=yorg*1000.
      dx=dx*1000.
      write(*,*) 'X orig, Y orig, cell size (m):',xorg,yorg,dx

      read(*,'(a)') fname
      open(11,file=fname,recl=1000)
      write(*,*)'Openend output text file:',fname

      read(*,'(a)') fname      
      open(9,file=fname,form='unformatted')
      write(*,*) 'Opened CAMx format point source file: ',fname

      read(*,'(a)') outfile
      write(*,*)'IOAPI output file:', outfile

      read(*,'(a)') fname
      read(fname,*) nfiles

      do n = 1, nfiles
      read(*,'(a)') ename(n)
      enddo

      do n = 1,nfiles
      write(*,*)n,ename(n)
      enddo

c---- open ioapi file: stack info

      land = ename(1)

      if ( .not. open3( land, FSREAD3, progname ) ) THEN
         MESG = 'Could not open file "' //
     &   land( 1: TRIMLEN(land))
     &   // '" for input'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

      IF ( .NOT. DESC3(land))THEN
         MESG = 'Could not get description info for file "' //
     &            land( 1: TRIMLEN( land) ) //'"'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      ENDIF

        NCOLS = NCOLS3D
        NROWS = NROWS3D
        NLAYS = NLAYS3D
        TSTEP = TSTEP3D

        write(*,*) NCOLS,NROWS,NLAYS,TSTEP

c---- allocate and clear arrays

      allocate ( TOTS(NCOLS, NROWS ,1000) )
      allocate ( QA(NX, NY ,1000) )
      allocate ( EMIS(1000, NROWS) )
      allocate ( zero(NROWS) )
      zero = 0.

C          DO V = 1,1000
C            DO R = 1, ny
C             DO C = 1, nx
C                 TOTS( C,R,V ) = 0.0
C                 QA(C,R,V) = 0.0
C             ENDDO
C            ENDDO
C          ENDDO

      tots = 0.
      qa   = 0.

      allocate ( GRID(NCOLS, NROWS, NLAYS) )
      allocate ( stkdm(NCOLS, NROWS))
      allocate ( stkht(NCOLS, NROWS))
      allocate ( stktk(NCOLS, NROWS))
      allocate ( stkve(NCOLS, NROWS))
      allocate ( stkflw(NCOLS, NROWS))
      allocate ( row(NCOLS, NROWS))
      allocate ( col(NCOLS, NROWS))
      allocate ( xloca(NCOLS, NROWS))
      allocate ( yloca(NCOLS, NROWS))

C.......   Get max string-lengths for use in variables-listing:

        VMAX = TRIMLEN( VNAME3D( 1 ) )
        UMAX = TRIMLEN( UNITS3D( 1 ) )
        DMAX = TRIMLEN( VDESC3D( 1 ) )
        DO  I = 1, NVARS3D
            VMAX = MAX( VMAX , TRIMLEN( VNAME3D( I ) ) )
            UMAX = MAX( UMAX , TRIMLEN( UNITS3D( I ) ) )
            DMAX = MAX( DMAX , TRIMLEN( VDESC3D( I ) ) )
        END DO

        WRITE( *,92000 )
     &  ' ', 'The list of variables in this file is:', ' ',
     &  ( VNAME3D( I )( 1:VMAX ) // ' (' //
     &    UNITS3D( I )( 1:UMAX ) // '): ' //
     &    VDESC3D( I )( 1:DMAX ), I = 1, NVARS3D )

        write(*,*) NVARS3D

                NVARS = NVARS3D
                DO  V = 1, NVARS3D
                 VNAME( V ) = VNAME3D( V )
                 UNITS( V ) = UNITS3D( V )
                 VDESC( V ) = VDESC3D( V )
                 EFLAG = EFLAG .OR. ( VTYPE3D( V ) .NE. M3REAL )
                END DO

c---- read in variables and sum over all hours

       JDATE = SDATE3D
       JTIME = STIME3D

        DO  V = 1, NVARS

            IF ( .NOT. READ3( land, VNAME( V ), ALLAYS3,
     &                        JDATE, JTIME, GRID ) ) THEN
                MESG = 'Read failure:  file ' // land //
     &                 ' variable ' // VNAME( V )
                CALL M3EXIT( 'VERTOT:VERSTEP', JDATE, JTIME,
     &                       MESG, 2 )
            END IF      !  if read3() failed

             if(VNAME(V).eq.'STKDM') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              stkdm( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'STKHT') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              stkht( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'STKTK') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              stktk( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'STKVE') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              stkve( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'STKFLW') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              stkflw( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'XLOCA') then   
           DO R = 1, NROWS
            DO C = 1, NCOLS
              xloca( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'YLOCA') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              yloca( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'ROW') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              row( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             elseif(VNAME(V).eq.'COL') then
           DO R = 1, NROWS
            DO C = 1, NCOLS
              col( C,R) =  GRID( C,R,1 )
            ENDDO
           ENDDO
             endif

        ENDDO !end loop over variables on file

       write(*,*) 'got here'

      if ( .not. close3 ( land ) ) THEN
         MESG = 'Could not close file'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if



c---- open ioapi file: emissions

      land = ename(2)

      if ( .not. open3( land, FSREAD3, progname ) ) THEN
         MESG = 'Could not open file "' //
     &   land( 1: TRIMLEN(land))
     &   // '" for input'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

      IF ( .NOT. DESC3(land))THEN
         MESG = 'Could not get description info for file "' //
     &            land( 1: TRIMLEN( land) ) //'"'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      ENDIF

        NCOLS = NCOLS3D
        NROWS = NROWS3D
        NLAYS = NLAYS3D
        TSTEP = TSTEP3D

        write(*,*) NCOLS,NROWS,NLAYS,TSTEP

C.......   Get max string-lengths for use in variables-listing:

        VMAX = TRIMLEN( VNAME3D( 1 ) )
        UMAX = TRIMLEN( UNITS3D( 1 ) )
        DMAX = TRIMLEN( VDESC3D( 1 ) )
        DO  I = 1, NVARS3D
            VMAX = MAX( VMAX , TRIMLEN( VNAME3D( I ) ) )
            UMAX = MAX( UMAX , TRIMLEN( UNITS3D( I ) ) )
            DMAX = MAX( DMAX , TRIMLEN( VDESC3D( I ) ) )
        END DO

        WRITE( *,92000 )
     &  ' ', 'The list of variables in this file is:', ' ',
     &  ( VNAME3D( I )( 1:VMAX ) // ' (' //
     &    UNITS3D( I )( 1:UMAX ) // '): ' //
     &    VDESC3D( I )( 1:DMAX ), I = 1, NVARS3D )

        write(*,*) NVARS3D

                NVARS = NVARS3D
                DO  V = 1, NVARS3D
                 VNAME( V ) = VNAME3D( V )
                 UNITS( V ) = UNITS3D( V )
                 VDESC( V ) = VDESC3D( V )
                 EFLAG = EFLAG .OR. ( VTYPE3D( V ) .NE. M3REAL )
                END DO

c---- camx file header

        write(*,*)'Number of stacks:',NROWS

         do V = 1,NVARS  ! convert species names
          cspec = VNAME(V)
           if(cspec.eq.'POC') cspec = 'POA'
           if(cspec.eq.'PMFINE') cspec = 'FPRM'
           if(cspec.eq.'PMC') cspec = 'CPRM'
           if(cspec.eq.'SESQ') cspec = 'SQT'
           if(cspec.eq.'UNR') cspec = 'NR'
           if(cspec.eq.'HGNRVA') cspec = 'HG0'
           if(cspec.eq.'PHGI') cspec = 'HGP'
           if(cspec.eq.'HGIIGAS') cspec = 'HG2'
           if(cspec.eq.'CL2') cspec = 'CHLORINE'
           if(cspec.eq.'HCL') cspec = 'HCL_GAS'
           if(cspec.eq.'TOL') dup(1) = V
           if(cspec.eq.'XYL') dup(2) = V
           if(cspec.eq.'ISOP') dup(3) = V
           if(cspec.eq.'TERP') dup(4) = V
          do i=1,10
          write(MSPEC(i,V),'(a)') cspec(i:i)
          enddo
         enddo

          cspec = 'TOLA'
          do i=1,10
          write(MSPEC(i,(NVARS+1)),'(a)') cspec(i:i)
          enddo
          cspec = 'XYLA'
          do i=1,10
          write(MSPEC(i,(NVARS+2)),'(a)') cspec(i:i)
          enddo
          cspec = 'ISP'
          do i=1,10
          write(MSPEC(i,(NVARS+3)),'(a)') cspec(i:i)
          enddo
          cspec = 'TRP'
          do i=1,10
          write(MSPEC(i,(NVARS+4)),'(a)') cspec(i:i)
          enddo

         tempa = int(SDATE3D/1000)*1000 !convert to camx date formate YYDDD
         tempb = SDATE3D - tempa !jday
         tempc = tempa - 2000000
         cdate = tempc + tempb

      WRITE(9) IFILE,NOTE,ione,(NVARS+4),cdate,real(0),cdate,real(23)
      WRITE(9) rzero,rzero,izero,XORG,YORG,DX,DX,NX,NY,
     $ ione,izero,izero,rzero,rzero,rzero
      WRITE (*,1007) cdate, real(0), cdate, real(23)
      WRITE (9) izero,izero,NX,NY
      WRITE (9) ((MSPEC(I,J),I=1,10),J=1,(NVARS+4))
      WRITE (*,1013) ((MSPEC(I,J),I=1,10),J=1,(NVARS+4))
      WRITE(9) ione,NROWS
      WRITE(9) (xloca(1,R),yloca(1,R),stkht(1,R),stkdm(1,R),stktk(1,R),
     &   stkve(1,R),R=1,NROWS)


c---- read in variables and sum over all hours

       JDATE = SDATE3D
       JTIME = STIME3D
       do HOURS = 0, 23
       JTIME = HOURS*10000
       write(*,*) 'Reading file at',JDATE,JTIME

        DO  V = 1, NVARS

            IF ( .NOT. READ3( land, VNAME( V ), ALLAYS3,
     &                        JDATE, JTIME, GRID ) ) THEN
                MESG = 'Read failure:  file ' // land //
     &                 ' variable ' // VNAME( V )
                CALL M3EXIT( 'VERTOT:VERSTEP', JDATE, JTIME,
     &                       MESG, 2 )
            END IF      !  if read3() failed

          IH = HOURS+1
           DO R = 1, NROWS
            DO C = 1, NCOLS
                EMIS( V,R) = GRID(C,R,1)
                TOTS( C,R,V ) = TOTS(C,R,V) + GRID( C,R,1 )
            ENDDO
           ENDDO

        ENDDO !end loop over variables on file

c---- write out to camx file for this hour

	WRITE (9) cdate,real(IH),cdate,real(IH+1)
        WRITE(*,1007) cdate,real(IH),cdate,real(IH+1)

             WRITE(9) ione,NROWS
             WRITE(9) (col(1,R),row(1,R),ione,(stkflw(1,R)*3600.),
     &          ione,R=1,NROWS)

              DO V=1,NVARS

              WRITE(9) ione,(MSPEC(I,V),I=1,10),((EMIS(V,R)*3600.),
     &             R=1,NROWS)

              ENDDO !end nspecs loop

              DO K=1,4 !duplicate species
                 if ( dup(k) .ne. 0 ) then
         WRITE(9) ione,(MSPEC(I,(NVARS+K)),I=1,10),((EMIS(DUP(K),R)*3600.),
     &             R=1,NROWS)
              else
         WRITE(9) ione,(MSPEC(I,(NVARS+K)),I=1,10),((zero(R)),
     &             R=1,NROWS)
              endif

              enddo !end loop over duplicate species for SOA chemistry

c----- done writing to camx file

       write(*,*) 'got to this point',HOURS
      enddo !end loop over hours

      if ( .not. close3 ( land ) ) THEN
         MESG = 'Could not close file'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

c-----create ascii domain total output

           DO R = 1, NROWS !rows and columns are empty on the IOAPI file
            DO C = 1, NCOLS
c              write(*,*) xloca(C,R),yloca(C,R),xorg,yorg,dx
              row(C,R) = int((yloca(C,R)-yorg)/dx)
              col(C,R) = int((xloca(C,R)-xorg)/dx)               
            ENDDO
           ENDDO

        write(11,1002)
 1002   format('specie,date,domain daily total,units')

            DO V = 1, NVARS
           cfac = 0.
           DO R = 1, NROWS
            DO C = 1, NCOLS
              cfac = cfac + TOTS( C,R,V )
               K = col(C,R)
               N = row(C,R)
               if ( k .ge. 1 .and. k .le. nx .and.
     &              n .ge. 1 .and. n .le. ny )
     &        QA(K,N,V) = QA(K,N,V) + TOTS( C,R,V )
            ENDDO
           ENDDO
            write(11,8080) VNAME(V),JDATE,cfac,UNITS(V)
            ENDDO
 8080       format(a20,1x,i10,1x,f20.10,1x,a10)        

 1007 FORMAT(2(I10,F10.2),I10)
 1013 FORMAT(1X,10A1)

c------ create output file
                 
c        IF ( TSTEP .EQ. 0 ) THEN
            JDATE  = 0
            JTIME  = 0
            NSTEPS = 1
c        ELSE
       
            NCOLS3D = NX
            NROWS3D = NY
            XORIG3D = XORG
            YORIG3D = YORG
            XCELL3D = DX
            YCELL3D = DX
            GDTYP3D = 2
            P_ALP3D = 33.
            P_BET3D = 45.
            P_GAM3D = -97.
            XCENT3D = -97.
            YCENT3D = 40.
            SDATE3D = JDATE
            STIME3D = JTIME
            NLAYS3D = 1
            NVARS3D = NVARS
            DO I = 1, NVARS
                VNAME3D( I ) = VNAME( I )
                UNITS3D( I ) = UNITS( I )
                VDESC3D( I ) = VDESC( I )
                VTYPE3D( I ) = M3REAL  
            END DO     
                
      if ( .not. open3( outfile, FSUNKN3, progname ) ) THEN
          MESG = 'Could not open file "' //
     &     outfile( 1: TRIMLEN(outfile))
     &     // '" for output'
           CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if    
                
          do V = 1,NVARS
           write(*,*)VNAME3D(V)

      IF ( .NOT. WRITE3( OUTFILE,VNAME(V),JDATE, JTIME, QA(1,1,V) )) THEN
                 MESG = 'Could not write totals to ' // VNAME(V)
                  CALL M3EXIT( 'VERTOT:VERSTEP', JDATE, JTIME,
     &                          MESG, 2 )
                END IF      !  if write3() failed
      
          enddo !end loop over variables



C******************  FORMAT  STATEMENTS   ******************************

C...........   Error and warning message formats..... 91xxx

91000   FORMAT ( //5X , '*** ERROR ABORT in program VERTOT ***',
     &            /5X , A ,
     &           // )        !  generic error message format

C...........   Informational (LOG) message formats... 92xxx

92000   FORMAT ( 5X , A )

92999   FORMAT ( //5X , A, // )

C...........   Formatted file I/O formats............ 93xxx

93000   FORMAT ( A16 )

C...........   Internal buffering formats............ 94xxx

C...........   Miscellaneous formats................. 95xxx

95000   FORMAT ( /5X , A , $ )          !  generic prompt format.


      stop

      END


c-----Start subroutines
      subroutine lcpgeo(iway,phic,xlonc,truelat1,truelat2,xloc,yloc,
     &                  xlon,ylat)
c
c     LCPGEO performs Lambert Conformal to geodetic (lat/lon) translation
c
c     Code based on the TERRAIN preprocessor for MM5 v2.0,
c     developed by Yong-Run Guo and Sue Chen, National Center for
c     Atmospheric Research, and Pennsylvania State University
c     10/21/1993
c
c     Input arguments:
c        iway                Conversion type
c                            0 = geodetic to Lambert Conformal
c                            1 = Lambert Conformal to geodetic
c        phic                Central latitude (deg, neg for southern hem)
c        xlonc               Central longitude (deg, neg for western hem)
c        truelat1            First true latitute (deg, neg for southern hem)
c        truelat2            Second true latitute (deg, neg for southern hem)
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
c     Output arguments:
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
      data conv/57.29578/, a/6370./
c
c-----Entry Point
c
      if (phic.lt.0) then
        sign = -1.
      else
        sign = 1.
      endif
      pole = 90.
      if (abs(truelat1).gt.90.) then
        truelat1 = 60.
        truelat2 = 30.
        truelat1 = sign*truelat1
        truelat2 = sign*truelat2
      endif
      xn = alog10(cos(truelat1/conv)) - alog10(cos(truelat2/conv))
      xn = xn/(alog10(tan((45. - sign*truelat1/2.)/conv)) -
     &         alog10(tan((45. - sign*truelat2/2.)/conv)))           
      psi1 = 90. - sign*truelat1
      psi1 = psi1/conv
      if (phic.lt.0.) then
        psi1 = -psi1
        pole = -pole
      endif
      psi0 = (pole - phic)/conv
      xc = 0.
      yc = -a/xn*sin(psi1)*(tan(psi0/2.)/tan(psi1/2.))**xn
c
c-----Calculate lat/lon of the point (xloc,yloc)
c
      if (iway.eq.1) then
        xloc = xloc + xc
        yloc = yloc + yc
        if (yloc.eq.0.) then
          if (xloc.ge.0.) flp = 90./conv
          if (xloc.lt.0.) flp = -90./conv
        else
          if (phic.lt.0.) then
            flp = atan2(xloc,yloc)
          else
            flp = atan2(xloc,-yloc)
          endif
        endif
        flpp = (flp/xn)*conv + xlonc
        if (flpp.lt.-180.) flpp = flpp + 360.
        if (flpp.gt. 180.) flpp = flpp - 360. 
        xlon = flpp 
c
        r = sqrt(xloc*xloc + yloc*yloc)
        if (phic.lt.0.) r = -r
        cell = (r*xn)/(a*sin(psi1))
        rxn  = 1.0/xn
        cel1 = tan(psi1/2.)*cell**rxn
        cel2 = atan(cel1)
        psx  = 2.*cel2*conv
        ylat = pole - psx
c
c-----Calculate x/y from lat/lon
c
      else
        ylon = xlon - xlonc
        if (ylon.gt. 180.) ylon = ylon - 360.
        if (ylon.lt.-180.) ylon = ylon + 360.
        flp = xn*ylon/conv
        psx = (pole - ylat)/conv
        r = -a/xn*sin(psi1)*(tan(psx/2.)/tan(psi1/2.))**xn
        if (phic.lt.0.) then
          xloc = r*sin(flp)
          yloc = r*cos(flp)
        else
          xloc = -r*sin(flp)
          yloc =  r*cos(flp)
        endif
      endif
c
      xloc = xloc - xc
      yloc = yloc - yc
c
      return
      end

      subroutine caldate(idate)

      integer idate
      dimension nday(12)
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/

c
c-----Entry point
c
c-----If it is already in calender date, return
c

      icent = int(idate/100000)
      iyear = int((idate - icent*100000)/1000)
      jday = idate - icent*100000 - iyear*1000

      nday(2) = 28
      if (mod(iyear,4).eq.0) nday(2) = 29
      mday = 0
      do imonth = 1,12
        mday = mday + nday(imonth)
        if (mday.ge.jday) go to 20
      enddo
 20   iday = jday - (mday - nday(imonth))
      idate = icent*1000000 + iyear*10000 + imonth*100 + iday
c
c      write(*,'(4i6)') icent,iyear,imonth,iday

      return
      end
