source: trunk/SRC/ReadWrite/readoldopadistcoast.pro

;+
;
; @file_comments
; read the old binary direct access file that contains the
; distance to the coast in OPA.
;
; based on the OPA subroutines dtacof and parctl
;
; @categories
; For OPA
;
; @param UNIT
;
; @param PARAMS
;
; @param NUM
;
; @returns
; a structure that contains two elements: tdistcoast (the
; distance for the t-points) and fdistcoast (the distance for the
; f-points).
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June, 2002
; @version
; $Id$
;
;-
FUNCTION read3fromopa, unit, params, num
;
  compile_opt idl2, strictarrsubs
;
   offset=params.reclen*params.jpk*(num-1L)
   a=assoc(unit,dblarr(params.jpiglo,params.jpjglo,params.jpk,/nozero),offset)
   b = a[0]
   return, b
end
;
;+
;
; @file_comments
;
; @categories
; For OPA
;
; @param FILENAME {in}{required}
; Filename (with the whole path if necessary)
;
; @param JPIGLO {in}{required}
;
; @param JPJGLO {in}{required}
;
; @param JPK {in}{required}
; Dimensions of the opa grid
;
; @keyword IBLOC {default=4096L}
; Ibloc size
;
; @keyword JPBYT {default=8L}
; Jpbyt size
;
; @keyword NUMREC {default=19L*jpk}
; Number of records in the file.
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June, 2002
;
; @version
; $Id$
;
;-
FUNCTION readoldopadistcoast, filename, jpiglo, jpjglo, jpk $
                            , IBLOC=ibloc, JPBYT=jpbyt, NUMREC=numrec
;
  compile_opt idl2, strictarrsubs
;
   iname_file = findfile(filename)
   if iname_file[0] EQ '' then begin
      ras = report( 'Bad file name')
      return, -1
   ENDIF ELSE iname_file = iname_file[0]
; open the file
   openr,numcost , iname_file, /get_lun, /swap_if_little_endian
; check the size of the file
   fileparameters = fstat(numcost)
; default parameter definition for ORCA2
   IF keyword_set(ibloc) THEN ibloc = long(ibloc) ELSE ibloc = 4096L
   jpiglo = long(jpiglo)
   jpjglo = long(jpjglo)
   jpk = long(jpk)
   IF keyword_set(jpbyt) THEN jpbyt = long(jpbyt) ELSE jpbyt = 8L
; record length computation
   reclen = ibloc*((jpiglo*jpjglo*jpbyt-1 )/ibloc+1)
; number of records
   IF keyword_set(numrec) THEN numrec = long(numrec) ELSE numrec = 3L*jpk
; difference between the record length and the size of the contained
; array.
   toomuch = reclen-jpiglo*jpjglo*jpbyt
; expected size computation
   size = numrec*reclen-toomuch
   if size NE fileparameters.size then begin
      ras = report(['The size of the file is not the expected one!', $
      'Check your file or the values of ibloc, jpiglo,', $
      'jpjglo, jpk, jpbyt, numrec in this program'])
      return, -1
   endif
; first record: six 64-bit integer to read.
; default definition
   iimlu = long64(999)
   ijmlu = long64(999)
   ikmlu = long64(999)
; read
   readu, numcost, iimlu, ijmlu, ikmlu
   if iimlu NE jpiglo then begin
      ras = report( 'iimlu = '+strtrim(iimlu, 1)+' differs from jpiglo ='+strtrim(jpiglo, 1))
      return, -1
   endif
   if ijmlu NE jpjglo then begin
      ras = report( 'ijmlu = '+strtrim(ijmlu, 1)+' differs from jpjglo ='+strtrim(jpjglo, 1))
      return, -1
   endif
   if ikmlu NE jpk then begin
      ras = report( 'ikmlu = '+strtrim(ikmlu, 1)+' differs from jpk ='+strtrim(jpk, 1))
      return, -1
   endif
; other records
   params = {jpiglo:jpiglo, jpjglo:jpjglo, jpk:jpk, reclen:reclen}
   tdistcoast = read3fromopa(numcost, params, 2)
   fdistcoast = read3fromopa(numcost, params, 3)
;
   close, numcost
   free_lun, numcost

   return, {tdistcoast:tdistcoast, fdistcoast:fdistcoast}
end