source: trunk/SRC/Grid/checkperio.pro @ 287

;+
;
; @file_comments
; check ORCA2 and ORCA05 periodicity (east-west and noth fold) for files related to the T GRID ONLY
;
; @categories
; Grid
;
; @param file {in}{required}{type=scalar string}
; A string giving the name of the NetCdf file
;
; @keyword WRITE {default=0}{type=scalar: 0 or 1}
; Define to 1 to write the good periodicity in the file.
; By default print only message of periodicity errors
;
; @restrictions
; Works only for grid T points
; if the WRITE keyword is activated, we need the writing access on the file!
;
; @examples
;    IDL> checkperio, '/Users/smasson/OPA/modipsl/bin/runoff_1m_nomask.nc'
;    checkperio, '/Users/smasson/OPA/modipsl/bin/runoff_1m_nomask.nc'
;    nav_lon bad y periodicity (1)
;    nav_lon bad y periodicity (2)
;    nav_lon bad x periodicity (1)
;    nav_lon bad x periodicity (2)
;    nav_lat bad y periodicity (1)
;    nav_lat bad y periodicity (2)
;    nav_lat bad x periodicity (1)
;    nav_lat bad x periodicity (2)
;    socoefr bad y periodicity (1)
;    socoefr bad x periodicity (1)
;    socoefr bad x periodicity (2)
;    sorunoff bad x periodicity (1)
;    sorunoff bad x periodicity (2)
;    IDL> checkperio, '/Users/smasson/OPA/modipsl/bin/runoff_1m_nomask.nc', /WRITE
;    IDL> checkperio, '/Users/smasson/OPA/modipsl/bin/runoff_1m_nomask.nc'
;
; @history
; October 2006: Sebastien Masson (smasson\@locean-ipsl.upmc.fr)
;
; @version
; $Id$
;-
;
PRO checkperio, file, WRITE = write
;
  compile_opt idl2, strictarrsubs
;
  IF file_test(file) EQ 0 THEN BEGIN
    ras = report( 'file '+file+' not found')
    RETURN
  ENDIF

  cdfid = ncdf_open(file, WRITE = write)
  finq = ncdf_inquire(cdfid)
;
  dimsz = lonarr(finq.ndims)
  FOR did = 0, finq.ndims-1 DO BEGIN
    ncdf_diminq, cdfid, did, name, size
    dimsz[did] = size
  ENDFOR
;
  FOR vid = 0, finq.nvars-1 DO BEGIN
    vinq =  ncdf_varinq(cdfid, vid)
    IF vinq.ndims GE 2 THEN BEGIN
      vdimsz = dimsz[vinq.dim]
      jpi = vdimsz[0]
      jpj = vdimsz[1]
      CASE vinq.ndims OF
;------------------------------
; 2D VAR
;------------------------------
        2:BEGIN
; y periodicity (1)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1], offset = [1, jpj-1]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1], offset = [1, jpj-3]
              part2 = reverse(shift(part2, -1))
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1], offset = [1, jpj-1]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1], offset = [1, jpj-1]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1], offset = [1, jpj-2]
              part2 = reverse(part2)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1], offset = [1, jpj-1]
              ENDIF
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; y periodicity (2)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [(jpi-1)/2-2+1, 1], offset = [(jpi-1)/2+2, jpj-2]
              ncdf_varget, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1], offset = [          2, jpj-2]
              part2 = reverse(part2)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (2)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1], offset = [(jpi-1)/2+2, jpj-2]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; x periodicity (1)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj], offset = [0, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj], offset = [jpi-2, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (1)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj], offset = [0, 0]
          ENDIF
; x periodicity (2)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj], offset = [jpi-1, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj], offset = [1, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (2)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj], offset = [jpi-1, 0]
          ENDIF
        END
;------------------------------
; 3D VAR
;------------------------------
        3:BEGIN
          jpk = vdimsz[2]
; y periodicity (1)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1, jpk], offset = [1, jpj-1, 0]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1, jpk], offset = [1, jpj-3, 0]
              IF jpk EQ 1 THEN part2 = reform(part2, jpi-2, 1, jpk, /over)
              part2 = reverse(shift(part2, -1, 0, 0), 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1, jpk], offset = [1, jpj-1, 0]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1, jpk], offset = [1, jpj-1, 0]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1, jpk], offset = [1, jpj-2, 0]
              IF jpk EQ 1 THEN part2 = reform(part2, jpi-2, 1, jpk, /over)
              part2 = reverse(part2, 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1, jpk], offset = [1, jpj-1, 0]
              ENDIF
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; y periodicity (2)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [(jpi-1)/2-2+1, 1, jpk], offset = [(jpi-1)/2+2, jpj-2, 0]
              ncdf_varget, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1, jpk], offset = [          2, jpj-2, 0]
              part2 = reverse(part2, 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (2)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1, jpk], offset = [(jpi-1)/2+2, jpj-2, 0]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; x periodicity (1)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj, jpk], offset = [0, 0, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj, jpk], offset = [jpi-2, 0, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (1)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj, jpk], offset = [0, 0, 0]
          ENDIF
; x periodicity (2)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj, jpk], offset = [jpi-1, 0, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj, jpk], offset = [1, 0, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (2)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj, jpk], offset = [jpi-1, 0, 0]
          ENDIF
        END
;------------------------------
; 4D VAR
;------------------------------
        4:BEGIN
          jpk = vdimsz[2]
          jpt = vdimsz[3]
; y periodicity (1)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-1, 0, 0]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-3, 0, 0]
              IF jpt EQ 1 THEN part2 = reform(part2, jpi-2, 1, jpk, jpt, /over)
              part2 = reverse(shift(part2, -1, 0, 0, 0), 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-1, 0, 0]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
              ncdf_varget, cdfid, vid, part1, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-1, 0, 0]
              ncdf_varget, cdfid, vid, part2, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-2, 0, 0]
              IF jpt EQ 1 THEN part2 = reform(part2, jpi-2, 1, jpk, jpt, /over)
              part2 = reverse(part2, 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (1)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [jpi-2, 1, jpk, jpt], offset = [1, jpj-1, 0, 0]
              ENDIF
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; y periodicity (2)
          CASE 1 OF
            array_equal([182, 149], vdimsz[0:1]):BEGIN ; ORCA2
              ncdf_varget, cdfid, vid, part1, count = [(jpi-1)/2-2+1, 1, jpk, jpt], offset = [(jpi-1)/2+2, jpj-2, 0, 0]
              ncdf_varget, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1, jpk, jpt], offset = [          2, jpj-2, 0, 0]
              part2 = reverse(part2, 1)
              IF array_equal(part1, part2) NE 1 THEN BEGIN
                ras = report( vinq.name+' bad y periodicity (2)')
                IF keyword_set(write) THEN $
                   ncdf_varput, cdfid, vid, part2, count = [(jpi-1)/2-2+1, 1, jpk, jpt], offset = [(jpi-1)/2+2, jpj-2, 0, 0]
              ENDIF
            END
            array_equal([722, 511], vdimsz[0:1]):BEGIN ; ORCA05
            END
            ELSE:ras = report( vinq.name+' nothing to check')
          ENDCASE
; x periodicity (1)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj, jpk, jpt], offset = [0, 0, 0, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj, jpk, jpt], offset = [jpi-2, 0, 0, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (1)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj, jpk, jpt], offset = [0, 0, 0, 0]
          ENDIF
; x periodicity (2)
          ncdf_varget, cdfid, vid, part1, count = [1, jpj, jpk, jpt], offset = [jpi-1, 0, 0, 0]
          ncdf_varget, cdfid, vid, part2, count = [1, jpj, jpk, jpt], offset = [1, 0, 0, 0]
          IF array_equal(part1, part2) NE 1 THEN BEGIN
            ras = report( vinq.name+' bad x periodicity (2)')
            IF keyword_set(write) THEN $
               ncdf_varput, cdfid, vid, part2, count = [1, jpj, jpk, jpt], offset = [jpi-1, 0, 0, 0]
          ENDIF
        END
        ELSE:ras = report( vinq.name+' nothing to check')
      ENDCASE
    ENDIF
  ENDFOR

  ncdf_close, cdfid

  RETURN
END