source: trunk/SRC/ReadWrite/readoldoparestart.pro @ 226

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @categories
; For OPA
;
; @param UNIT
;
;
; @param PARAMS
;
; @param NUM
;
; @restrictions
; bug for etab and etan written on the same record???
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June 2002
;
; @version
; $Id$
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------


FUNCTION read2fromopa, unit, params, num
;
  compile_opt idl2, strictarrsubs
;
   offset=params.reclen*params.jpk*(num-1L)
   a=assoc(unit,dblarr(params.jpiglo,params.jpjglo,/nozero),offset)
   return, a[0]
end
;+
; @categories
; For OPA
;
; @param UNIT
;
;
; @param PARAMS
;
; @param NUM
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June 2002
;-
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)
   return, a[0]
end

;+
; @file_comments
; read the old restart files of OPA (before NetCDF)
; based on the OPA subroutine dtrlec included at the end of the file.
;
; @categories
; For OPA
;
; @param FILENAME {in}{required}
; 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
;
; @keyword UB
;
;
; @keyword VB
;
;
; @keyword TB
;
; @keyword SB
;
; @keyword ROTB
;
; @keyword HDIVB
;
; @keyword UN
;
; @keyword VN
;
; @keyword TN
;
; @keyword SN
;
; @keyword ROTN
;
; @keyword HDIVN
;
; @keyword GCX
;
; @keyword GCXB
;
; @keyword ETAB
;
; @keyword ETAN
;
; @keyword BSFB
;
; @keyword BSFN
;
; @keyword BSFD
;
; @keyword EN
; the variable we want to read.
;
; @returns
; According to the given keywords.
;
; @restrictions
; Bug for etab and etan written on the same record???
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      June 2002
;
; @version
; $Id$
;-

PRO readoldoparestart, filename, jpiglo, jpjglo, jpk, IBLOC = ibloc, JPBYT = jpbyt, NUMREC = numrec, ub = ub, vb = vb, tb = tb, sb = sb, rotb = rotb, hdivb = hdivb, un = un, vn = vn, tn = tn, sn = sn, rotn = rotn, hdivn = hdivn, gcx = gcx, gcxb = gcxb, etab = etab, etan = etan, bsfb = bsfb, bsfn = bsfn, bsfd = bsfd, en = en
;
;
  compile_opt idl2, strictarrsubs
;
   iname_file = findfile(filename)
   if iname_file[0] EQ '' then begin
      print, 'Bad file name'
      return
   ENDIF ELSE iname_file = iname_file[0]
; open the file
   openr,numrst , iname_file, /get_lun, /swap_if_little_endian
; check the size of the file
   filepamameters = fstat(numrst)
; parameter definition
   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)
   IF keyword_set(numrec) THEN numrec = long(numrec) ELSE numrec = 19L*jpk
   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
   endif
; first record: six 64-bit integer to read.
; default definition
   ino1 = long64(9999)
   it1 = long64(9999)
   isor1 = long64(9999)
   ipcg1 = long64(9999)
   itke1 = long64(9999)
   idast1 = long64(9999)
; read
   readu, numrst, ino1, it1, isor1, ipcg1, itke1, idast1
   print, ino1, it1, isor1, ipcg1, itke1, idast1
; other records
   params = {jpiglo:jpiglo, jpjglo:jpjglo, jpk:jpk, reclen:reclen}
;      CALL read3(numrst,ub   ,2 )
   IF arg_present(ub) THEN ub = read3fromopa(numrst, params, 2)
;      CALL read3(numrst,vb   ,3 )
   IF arg_present(vb) THEN vb = read3fromopa(numrst, params, 3)
;      CALL read3(numrst,tb   ,5 )
   IF arg_present(tb) THEN tb = read3fromopa(numrst, params, 5)
;      CALL read3(numrst,sb   ,6 )
   IF arg_present(sb) THEN sb = read3fromopa(numrst, params, 6)
;      CALL read3(numrst,rotb ,7 )
   IF arg_present(rotb) THEN rotb = read3fromopa(numrst, params, 7)
;      CALL read3(numrst,hdivb,8 )
   IF arg_present(hdivb) THEN hdivb = read3fromopa(numrst, params, 8)
;      CALL read3(numrst,un   ,9 )
   IF arg_present(un) THEN un = read3fromopa(numrst, params, 9)
;      CALL read3(numrst,vn   ,10)
   IF arg_present(vn) THEN vn = read3fromopa(numrst, params, 10)
;      CALL read3(numrst,tn   ,12)
   IF arg_present(tn) THEN tn = read3fromopa(numrst, params, 12)
;      CALL read3(numrst,sn   ,13)
   IF arg_present(sn) THEN sn = read3fromopa(numrst, params, 13)
;      CALL read3(numrst,rotn ,14)
   IF arg_present(rotn) THEN rotn = read3fromopa(numrst, params, 14)
;      CALL read3(numrst,hdivn,15)
   IF arg_present(hdivn) THEN hdivn = read3fromopa(numrst, params, 15)
;C
;C ... Read elliptic solver arrays
;C
;      CALL read2(numrst,gcx ,jpk,17)
   IF arg_present(gcx) THEN gcx = read2fromopa(numrst, params, 17)
;      CALL read2(numrst,gcxb,jpk,18)
   IF arg_present(gcxb) THEN gcxb = read2fromopa(numrst, params, 18)
;C
;#ifdef key_freesurf_cstvol
;C
;C ... free surface formulation (eta)
;C
;      CALL read2(numrst,etab ,jpk,4 )
   IF arg_present(etab) THEN etab = read2fromopa(numrst, params, 4)
;      CALL read2(numrst,etan ,jpk,4 )
   IF arg_present(etan) THEN etan = read2fromopa(numrst, params, 4)
;#  else
;C
;C ... Rigid-lid formulation (bsf)
;C
;      CALL read2(numrst,bsfb ,jpk,4 )
   IF arg_present(bsfb) THEN bsfb = read2fromopa(numrst, params, 4)
;      CALL read2(numrst,bsfn ,jpk,11)
   IF arg_present(bsfn) THEN bsfn = read2fromopa(numrst, params, 11)
;      CALL read2(numrst,bsfd ,jpk,16)
   IF arg_present(bsfd) THEN bsfd  = read2fromopa(numrst, params, 16)
;#endif
;#ifdef key_zdftke
;          CALL read3(numrst,en,19)
   IF arg_present(en) THEN en = read3fromopa(numrst, params, 19)



   close, numrst
   free_lun, numrst

   return
end



;CDIR$ LIST
;      SUBROUTINE dtrlec
;CCC---------------------------------------------------------------------
;CCC
;CCC                       ROUTINE dtrlec
;CCC                     ******************
;CCC
;CCC  Purpose :
;CCC  --------
;CCC     Read files for restart
;CCC
;CC   Method :
;CC   -------
;CC      Read the previous fields on the file numrst
;CC      the first record indicates previous characteristics
;CC      after control with the present run, we read :
;CC      - prognostic variables on the second record
;CC      - elliptic solver arrays
;CC     - barotropic stream function arrays (default option)
;CC       or free surface arrays ("key_freesurf_cstvol" defined)
;CC      - tke arrays ("key_zdftke" defined)
;CC      for this last three records,  the previous characteristics
;CC      could be different with those used in the present run.
;CC
;CC   Input :
;CC   ------
;CC      common
;CC            /comrst/          : restart parameter
;CC            /comctl/          : parameters for the control
;CC
;CC   Output :
;CC   -------
;CC      common
;CC            /combef/          : previous fields (before)
;CC            /comnow/          : present fields (now)
;CC            /combsf/          : barotropic stream function
;CC            /comspg/          : surface pressure
;CC            /comsol/          : diagonal preconditioned conjugate
;CC
;CC   Modifications :
;CC   --------------
;CC      original  : 91-03 ()
;CC      additions : 92-01 (M. Imbard)
;CC                : 92-06 correction restart file (M. Imbard)
;CC                : 98-02 (M. Guyon) FETI method
;CC      addition  : 98-05 (G. Roullet) free surface
;CC----------------------------------------------------------------------
;CC parameters and commons
;CC ======================
;CDIR$ NOLIST
;#include "parameter.h"
;#include "common.h"
;CDIR$ LIST
;CC----------------------------------------------------------------------
;CC local declarations
;CC ==================
;      INTEGER ji, jj, jk, jl
;      INTEGER ino0, it0, ipcg0, isor0, itke0
;      INTEGER ino1, it1, isor1, ipcg1, itke1, idast1
;CC----------------------------------------------------------------------
;CC statement functions
;CC ===================
;CDIR$ NOLIST
;#include "stafun.h"
;CDIR$ LIST
;CCC---------------------------------------------------------------------
;CCC  OPA8, LODYC (1997)
;CCC---------------------------------------------------------------------
;C
;C
;C 0. Initialisations
;C ------------------
;C
;      ino0  = no
;      it0   = nit000
;      ipcg0 = 0
;      isor0 = 0
;      itke0 = 0
;      isor0 = nsolv-1
;      ipcg0 = 2-nsolv
;#ifdef key_zdftke
;      itke0 = 1
;#endif
;C FETI method
;      IF (nsolv .EQ. 3) THEN
;          isor0=2
;          ipcg0=2
;      ENDIF
;C
;      IF(lwp) THEN
;          WRITE(numout,*) ' '
;          WRITE(numout,*) ' *** dtrlec:  beginning of restart'
;          WRITE(numout,*) ' '
;          WRITE(numout,*) ' the present run :'
;          WRITE(numout,*) '   job number : ', no
;          WRITE(numout,*) '   with nit000 : ', nit000
;          WRITE(numout,*) '   with pcg option ipcg0 : ', ipcg0
;          WRITE(numout,*) '   with sor option isor0 : ', isor0
;          WRITE(numout,*) '   with FETI solver option ipcg0 & isor0 : ',
;     $        ipcg0,' & ',isor0
;          WRITE(numout,*) '   with tke option itke0 : ', itke0
;      ENDIF
;C
;C 1. Read numrst
;C --------------
;C
;C ... First record
;C
;      READ(numrst,REC=1) ino1, it1, isor1, ipcg1, itke1, idast1
;C
;      IF(lwp) THEN
;          WRITE(numout,*) ' '
;          WRITE(numout,*) ' READ numrst with '
;          WRITE(numout,*) '   job number : ', ino1
;          WRITE(numout,*) '   with time step it : ', it1
;          WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
;          WRITE(numout,*) '   with sor option isor1 : ', isor1
;          WRITE(numout,*) '   with tke option itke1 : ', itke1
;          WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ',
;     $        ipcg1 + isor1
;          WRITE(numout,*) ' '
;      ENDIF
;C
;C ... Control of date
;C
;      IF ( (it0-it1).NE.1 .AND. abs(nrstdt).EQ.1 ) THEN
;          IF(lwp) THEN
;              WRITE(numout,*) ' ===>>>> : problem with nit000 for the',
;     $            ' restart'
;              WRITE(numout,*) ' =======                              ',
;     $            ' ======='
;              WRITE(numout,*) ' we stop. verify the file'
;              WRITE(numout,*) ' or rerun with the value  0 for the'
;              WRITE(numout,*) ' control of time parameter  nrstdt'
;              WRITE(numout,*) ' '
;          ENDIF
;          STOP 'dtrlec'
;      ENDIF
;      IF ( nrstdt.EQ.1 ) ndate0 = idast1
;C
;C ... Read prognostic variables
;C
;      CALL read3(numrst,ub   ,2 )
;      CALL read3(numrst,vb   ,3 )
;      CALL read3(numrst,tb   ,5 )
;      CALL read3(numrst,sb   ,6 )
;      CALL read3(numrst,rotb ,7 )
;      CALL read3(numrst,hdivb,8 )
;      CALL read3(numrst,un   ,9 )
;      CALL read3(numrst,vn   ,10)
;      CALL read3(numrst,tn   ,12)
;      CALL read3(numrst,sn   ,13)
;      CALL read3(numrst,rotn ,14)
;      CALL read3(numrst,hdivn,15)
;C
;C ... Read elliptic solver arrays
;C
;      CALL read2(numrst,gcx ,jpk,17)
;      CALL read2(numrst,gcxb,jpk,18)
;C
;#ifdef key_freesurf_cstvol
;C
;C ... free surface formulation (eta)
;C
;      CALL read2(numrst,etab ,jpk,4 )
;      CALL read2(numrst,etan ,jpk,4 )
;#  else
;C
;C ... Rigid-lid formulation (bsf)
;C
;      CALL read2(numrst,bsfb ,jpk,4 )
;      CALL read2(numrst,bsfn ,jpk,11)
;      CALL read2(numrst,bsfd ,jpk,16)
;#endif
;C
;#ifdef key_zdftke
;C
;C ... Read tke arrays
;C
;      IF(itke1.eq.1) THEN
;          CALL read3(numrst,en,19)
;      ELSE
;          IF(lwp) THEN
;              WRITE(numout,*) ' ===>>>> : the previous restart file',
;     $            ' did''nt used  tke scheme'
;              WRITE(numout,*) ' =======                ======='
;          ENDIF
;          nrstdt=2
;      ENDIF
;#endif
;C
;C
;      RETURN
;      END