source: trunk/SRC/ReadWrite/readoldopadistcoast.pro @ 106

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; @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 before NetCDF
;
;
; @returns  a structure that contains two elements: tdistcoast (the
; distance for the t-points) and fdiscoast (the distance for the
; f-points).
;
; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June, 2002
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------

FUNCTION read3fromopa, unit, params, num
   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

;+
;     @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 ibloc size, default: ibloc = 4096L
;     @keyword JPBYT jpbyt size, default: jpbyt = 8L
;     @keyword NUMREC number of records in the file. default: numrec = 19L*jpk
; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June, 2002
;-
FUNCTION readoldopadistcoast, filename, jpiglo, jpjglo, jpk, IBLOC = ibloc, JPBYT = jpbyt, NUMREC = numrec
;
   iname_file = findfile(filename)
   if iname_file[0] EQ '' then begin
      print, '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
   filepamameters = fstat(numcost)
; defaut 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 contened
; array.
   toomuch = reclen-jpiglo*jpjglo*jpbyt
; expected size computation
   size = numrec*reclen-toomuch
   if size NE filepamameters.size then begin
      print, 'The size of the file is not the expected one!'
      print, 'Check your file or the values of ibloc, jpiglo,'
      print, '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
      print, 'iimlu = '+strtrim(iimlu, 1)+' differs from jpiglo ='+strtrim(jpiglo, 1) 
      return, -1
   endif
   if ijmlu NE jpjglo then begin
      print, 'ijmlu = '+strtrim(ijmlu, 1)+' differs from jpjglo ='+strtrim(jpjglo, 1) 
      return, -1
   endif
   if ikmlu NE jpk then begin
      print, '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