source: branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/nclatlon.F90 @ 2893

PROGRAM nclatlon
   USE netcdf
   IMPLICIT NONE
   !
   ! Command line arguments for output file and input file
   !
#ifndef NOIARGCPROTO
   INTEGER,EXTERNAL :: iargc
#endif
   INTEGER :: nargs
   CHARACTER(len=256) :: cdinfile
   CHARACTER(len=64) :: cdlatvar
   CHARACTER(len=64) :: cdlonvar
   CHARACTER(len=256) :: cdoutfile
   CHARACTER(len=32) :: cdscale
   !
   ! Input data
   !
   REAL, ALLOCATABLE, DIMENSION(:) :: zlat, zlon
   REAL :: zscalat,zscalon
   INTEGER :: nopoint
   !
   ! NetCDF ids
   !
   INTEGER :: idfile, idlat, idlon
   INTEGER, ALLOCATABLE, DIMENSION(:) :: idimlats, idimlons, idimlens
   INTEGER :: ndimlat, ndimlon
   !
   ! Loop variables
   !
   INTEGER :: i, isize
   !
   ! Get number of command line arguments
   !
   nargs=IARGC()
   IF ((nargs<4).OR.(nargs>5)) THEN
      WRITE(*,'(A)')'Usage:'
      WRITE(*,'(A)')'nclatlon inputfile latvar lonvar outputfile [scale]'
      CALL abort()
   ENDIF
   CALL getarg(1,cdinfile)
   CALL getarg(2,cdlatvar)
   CALL getarg(3,cdlonvar)
   CALL getarg(4,cdoutfile)
   WRITE(*,*)'cdinfile  = ',TRIM(cdinfile)
   WRITE(*,*)'cdlatvar  = ',TRIM(cdlatvar)
   WRITE(*,*)'cdlonvar  = ',TRIM(cdlonvar)
   WRITE(*,*)'cdoutfile = ',TRIM(cdoutfile)
   !
   IF (nargs==5) THEN
      CALL getarg(5,cdscale)
      WRITE(*,*)'cdscale = ',TRIM(cdscale)
      READ(cdscale,*)zscalat
      zscalon=zscalat
   ELSE
      zscalat=1.0
      zscalon=1.0
   ENDIF
   !
   ! Open netCDF file
   ! 
   IF (nf90_open(TRIM(cdinfile),nf90_nowrite,idfile)/=nf90_noerr) THEN
      WRITE(*,*)'Error opening input file'
      CALL abort
   ENDIF
   !
   ! Get latitude id
   !
   IF (nf90_inq_varid(idfile,TRIM(cdlatvar),idlat)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting lat id'
      CALL abort
   ENDIF
   !
   ! Get number of latitude dimensions
   !
   IF (nf90_inquire_variable(idfile,idlat,ndims=ndimlat)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting number of latitude dimensions'
      CALL abort
   ENDIF
   WRITE(*,*)'ndimlat=',ndimlat
   !
   ! Get latitude dimensions
   !
   ALLOCATE(idimlats(ndimlat))
   IF (nf90_inquire_variable(idfile,idlat,dimids=idimlats)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting latitude dimensions'
      CALL abort
   ENDIF
   WRITE(*,*)'idimlats=',idimlats
   !
   ! Get longitude id
   !
   IF (nf90_inq_varid(idfile,TRIM(cdlonvar),idlon)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting lon id'
      CALL abort
   ENDIF
   !
   ! Get number of longitude dimensions
   !
   IF (nf90_inquire_variable(idfile,idlon,ndims=ndimlon)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting number of longitude dimensions'
      CALL abort
   ENDIF
   WRITE(*,*)'ndimlon=',ndimlon
   !
   ! Get longitude dimensions
   !
   ALLOCATE(idimlons(ndimlon))
   IF (nf90_inquire_variable(idfile,idlon,dimids=idimlons)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting longitude dimensions'
      CALL abort
   ENDIF
   WRITE(*,*)'idimlons=',idimlons
   !
   ! Check that latitude and longitude has the same shape
   !
   IF (ndimlon/=ndimlat) THEN
      WRITE(*,*)'Different number of dimensions for latitude and longitude'
      CALL abort
   ENDIF
   DO i=1,ndimlat
      IF (idimlons(i)/=idimlats(i)) THEN
         WRITE(*,*)'Different dimension for latitude and longitude'
         CALL abort
      ENDIF
   ENDDO
   !
   ! Get dimensions
   !
   ALLOCATE(idimlens(ndimlat))
   DO i=1,ndimlat
      IF (nf90_inquire_dimension(idfile,idimlats(i),len=idimlens(i))/=nf90_noerr) THEN
         WRITE(*,*)'Error getting dimension length'
         CALL abort
      ENDIF
   ENDDO
   WRITE(*,*)'idimlens=',idimlens
   !
   ! Get the data
   !
   isize=1
   DO i=1,ndimlat
      isize=isize*idimlens(i)
   ENDDO
   ALLOCATE(zlat(isize),zlon(isize))
   IF (nf90_get_var(idfile,idlat,zlat)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting latitude'
      CALL abort
   ENDIF
   IF (nf90_get_var(idfile,idlon,zlon)/=nf90_noerr) THEN
      WRITE(*,*)'Error getting longitude'
      CALL abort
   ENDIF
   ! 
   ! Get the scale_factor for latitude (if it exists)
   !
   IF (zscalat==1.0) THEN
      IF (nf90_inquire_attribute(idfile,idlat,&
         &                     "scale_factor")==nf90_noerr) THEN
         IF (nf90_get_att(idfile,idlat,"scale_factor",zscalat) &
            & /=nf90_noerr) THEN
            WRITE(*,*)'Error getting latitude scale_factor'
            CALL abort
         ENDIF
         WRITE(*,*)'Latitude scale factor',zscalat
      ENDIF
   ENDIF
   ! 
   ! Get the scale_factor for longitude (if it exists)
   !
   IF (zscalon==1.0) THEN
      IF (nf90_inquire_attribute(idfile,idlon,&
         &                     "scale_factor")==nf90_noerr) THEN
         IF (nf90_get_att(idfile,idlon,"scale_factor",zscalon) &
            & /=nf90_noerr) THEN
            WRITE(*,*)'Error getting longitude scale_factor'
            CALL abort
         ENDIF
         WRITE(*,*)'Longitude scale factor',zscalon
      ENDIF
   ENDIF
   !
   ! Close netCDF file
   ! 
   IF (nf90_close(idfile)/=nf90_noerr) THEN
      WRITE(*,*)'Error closing input file'
      CALL abort
   ENDIF
   !
   ! Write the data
   !
   OPEN(10,file=cdoutfile)
   DO i=1,isize
      WRITE(10,'(2F12.6)')zscalat*zlat(i),zscalon*zlon(i)
   ENDDO
   CLOSE(10)
   !
END PROGRAM nclatlon