MODULE flincom
!-
!$Id: flincom.f90 2281 2010-10-15 14:21:13Z smasson $
!-
! This software is governed by the CeCILL license
! See IOIPSL/IOIPSL_License_CeCILL.txt
!---------------------------------------------------------------------
  USE netcdf
!-
  USE calendar,  ONLY : ju2ymds, ymds2ju, ioconf_calendar
  USE errioipsl, ONLY : histerr
  USE stringop,  ONLY : strlowercase
!-
  IMPLICIT NONE
!-
  PRIVATE
  PUBLIC :: flinput, flincre, flinget, flinclo, &
            flinopen, flininfo, flininspect, flinquery_var
!-
  INTERFACE flinopen
!---------------------------------------------------------------------
!- The "flinopen" routines will open an input file
!-
!- INPUT
!-
!- filename  : Name of the netCDF file to be opened
!-
!- iideb     : index i for zoom     !
!- iilen     : length  of  zoom     !   for
!- jjdeb     : index j for zoom     !   zoom
!- jjlen     : length  of  zoom     !
!-
!- do_test   : A flag that enables the testing of the content
!-             of the file against the input from the model
!-
!- INPUT if do_test=TRUE      OUTPUT else
!-
!- iim       : size in the x direction in the file (longitude)
!- jjm       : size in the y direction
!- llm       : number of levels
!-             (llm = 0 means no axis to be expected)
!- lon       : array of (iilen,jjlen) (zoom), or (iim,jjm) (no zoom),
!-             that contains the longitude of each point
!- lat       : same for latitude
!- lev       : An array of llm for the latitude
!-
!- WARNING :
!- In the case of do_test=FALSE it is for the user to check
!- that the dimensions of lon lat and lev are correct when passed to
!- flinopen. This can be done after the call when iim and jjm have
!- been retrieved from the netCDF file. In F90 this problem will
!- be solved with an internal assign
!- IF iim, jjm, llm or ttm are parameters in the calling program and
!- you use the option do_test=FALSE it will create a segmentation fault
!-
!-  OUTPUT
!-
!-  ttm       : size of time axis
!-  itaus     : Time steps within this file
!-  date0     : Julian date at which itau = 0
!-  dt        : length of the time steps of the data
!-  fid       : returned file ID which is later used to read the data
!---------------------------------------------------------------------
    MODULE PROCEDURE flinopen_zoom2d, flinopen_nozoom
  END INTERFACE
!-
  INTERFACE flinput
!---------------------------------------------------------------------
!- The "flinput" routines will put a variable
!- on the netCDF file created by flincre.
!- If the sizes of the axis do not match the one of the IDs
!- then a new axis is created.
!- That is we loose the possibility of writting hyperslabs of data.
!-
!- Again here if iim = jjm = llm = ttm = 0
!- then a global attribute is added to the file.
!-
!- INPUT
!-
!- fid      : Identification of the file in which we will write
!- varname  : Name of variable to be written
!- iim      : size in x of variable
!- nlonid   : ID of x axis which could fit for this axis
!- jjm      : size in y of variable
!- nlatid   : ID of y axis which could fit for this axis
!- llm      : size in z of variable
!- zdimid   : ID of z axis which could fit for this axis
!- ttm      : size in t of variable
!- tdimid   : ID of t axis which could fit for this axis
!-
!- OUTPUT
!-
!- NONE
!---------------------------------------------------------------------
    MODULE PROCEDURE flinput_r4d, flinput_r3d, flinput_r2d, &
                     flinput_r1d, flinput_scal
  END INTERFACE
!-
  INTERFACE flinget
    MODULE PROCEDURE flinget_r4d, flinget_r3d, flinget_r2d, &
                     flinget_r1d, flinget_scal, &
                     flinget_r4d_zoom2d, flinget_r3d_zoom2d, &
                     flinget_r2d_zoom2d
  END INTERFACE
!-
! This is the data we keep on each file we open
!-
  INTEGER, PARAMETER :: nbfile_max = 200
  INTEGER, SAVE :: nbfiles = 0
  INTEGER, SAVE :: ncids(nbfile_max), ncnbd(nbfile_max), &
                   ncfunli(nbfile_max), ncnba(nbfile_max)
  INTEGER, SAVE :: ncnbva(nbfile_max), ncdims(nbfile_max,4)
  LOGICAL, SAVE :: ncfileopen(nbfile_max)=.FALSE.
!-
  INTEGER, SAVE :: cind_vid, cind_fid, cind_len
  INTEGER,DIMENSION(:),ALLOCATABLE,SAVE :: cindex
!-
  INTEGER,DIMENSION(4) :: w_sta, w_len, w_dim
!-
CONTAINS
!-
!===
!-
SUBROUTINE flincre &
  (filename, iim1, jjm1, lon1, lat1, llm1, lev1, ttm1, itaus, &
   time0, dt, fid_out, nlonid1, nlatid1, zdimid1, tdimid1)
!---------------------------------------------------------------------
!- This is a "low level" subroutine for opening netCDF files wich
!- contain the major coordinate system of the model.
!- Other coordinates needed for other variables
!- will be added as they are needed.
!-
!- INPUT
!-
!- filename    : Name of the file to be created
!- iim1, jjm1  : Horizontal size of the grid
!-               which will be stored in the file
!- lon1, lat1  : Horizontal grids
!- llm1        : Size of the vertical grid
!- lev1        : Vertical grid
!- ttm1        : Size of time axis
!- itaus       : time steps on the time axis
!- time0       : Time in julian days at which itau = 0
!- dt          : time step in seconds between itaus
!-               (one step of itau)
!-
!- OUTPUT
!-
!- fid         : File identification
!- nlonid1     : Identification of longitudinal axis
!- nlatid1     : Identification of latitudinal axis
!- zdimid1     : ID of vertical axis
!- tdimid1     : ID of time axis
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  INTEGER :: iim1, jjm1, llm1, ttm1
  REAL :: lon1(iim1,jjm1)
  REAL :: lat1(iim1,jjm1)
  REAL :: lev1(llm1)
  INTEGER :: itaus(ttm1)
  REAL :: time0
  REAL :: dt
  INTEGER :: fid_out, zdimid1, nlonid1, nlatid1, tdimid1
!-
! LOCAL
!-
  INTEGER :: iret, lll, fid
  INTEGER :: lonid, latid, levid, timeid
  INTEGER :: year, month, day
  REAL :: sec
  CHARACTER(LEN=250):: name
!-
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  lll = LEN_TRIM(filename)
  IF (filename(lll-2:lll) /= '.nc') THEN
    name=filename(1:lll)//'.nc'
  ELSE
    name=filename(1:lll)
  ENDIF
!-
  iret = NF90_CREATE (name, NF90_CLOBBER, fid)
!-
  iret = NF90_DEF_DIM (fid, 'x',     iim1, nlonid1)
  iret = NF90_DEF_DIM (fid, 'y',     jjm1, nlatid1)
  iret = NF90_DEF_DIM (fid, 'lev',   llm1, zdimid1)
  iret = NF90_DEF_DIM (fid, 'tstep', ttm1, tdimid1)
!-
! Vertical axis
!-
  IF (check) WRITE(*,*) 'flincre Vertical axis'
!-
  iret = NF90_DEF_VAR (fid, 'lev', NF90_FLOAT, zdimid1, levid)
  iret = NF90_PUT_ATT (fid, levid, 'units',     '-')
  iret = NF90_PUT_ATT (fid, levid, 'title',     'levels')
  iret = NF90_PUT_ATT (fid, levid, 'long_name', 'Sigma Levels')
!-
! Time axis
!-
  IF (check) WRITE(*,*) 'flincre time axis'
!-
  iret = NF90_DEF_VAR (fid, 'tstep', NF90_FLOAT, tdimid1, timeid)
  iret = NF90_PUT_ATT (fid, timeid, 'units',     '-')
  iret = NF90_PUT_ATT (fid, timeid, 'title',     'time')
  iret = NF90_PUT_ATT (fid, timeid, 'long_name', 'time steps')
!-
! The longitude
!-
  IF (check) WRITE(*,*) 'flincre Longitude axis'
!-
  iret = NF90_DEF_VAR (fid, "nav_lon", NF90_FLOAT, &
                       (/ nlonid1, nlatid1 /), lonid)
  iret = NF90_PUT_ATT (fid, lonid, 'units', "degrees_east")
  iret = NF90_PUT_ATT (fid, lonid, 'title', "Longitude")
  iret = NF90_PUT_ATT (fid, lonid, 'nav_model', &
                       "Lambert projection of PROMES")
  iret = NF90_PUT_ATT (fid, lonid, 'valid_min', &
                       REAL(MINVAL(lon1),KIND=4))
  iret = NF90_PUT_ATT (fid, lonid, 'valid_max', &
                       REAL(MAXVAL(lon1),KIND=4))
!-
! The Latitude
!-
  IF (check) WRITE(*,*) 'flincre Latitude axis'
!-
  iret = NF90_DEF_VAR (fid, "nav_lat", NF90_FLOAT, &
                     (/ nlonid1, nlatid1 /), latid)
  iret = NF90_PUT_ATT (fid, latid, 'units', "degrees_north")
  iret = NF90_PUT_ATT (fid, latid, 'title', "Latitude")
  iret = NF90_PUT_ATT (fid, latid, 'nav_model', &
                       "Lambert projection of PROMES")
  iret = NF90_PUT_ATT (fid, latid, 'valid_min', &
                       REAL(MINVAL(lat1),KIND=4))
  iret = NF90_PUT_ATT (fid, latid, 'valid_max', &
                       REAL(MAXVAL(lat1),KIND=4))
!-
! The time coordinates
!-
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, 'delta_tstep_sec', &
                       REAL(dt,KIND=4))
!-
  CALL ju2ymds (time0, year, month, day, sec)
!-
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, 'year0',  REAL(year,KIND=4))
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, 'month0', REAL(month,KIND=4))
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, 'day0',   REAL(day,KIND=4))
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, 'sec0',   REAL(sec,KIND=4))
!-
  iret = NF90_ENDDEF (fid)
!-
  IF (check) WRITE(*,*) 'flincre Variable'
!-
  iret = NF90_PUT_VAR (fid, levid, lev1(1:llm1))
!-
  IF (check) WRITE(*,*) 'flincre Time Variable'
!-
  iret = NF90_PUT_VAR (fid, timeid, REAL(itaus(1:ttm1)))
!-
  IF (check) WRITE(*,*) 'flincre Longitude'
!-
  iret = NF90_PUT_VAR (fid, lonid, lon1(1:iim1,1:jjm1))
!-
  IF (check) WRITE(*,*) 'flincre Latitude'
!-
  iret = NF90_PUT_VAR (fid, latid, lat1(1:iim1,1:jjm1))
!-
! Keep all this information
!-
  nbfiles = nbfiles+1
!-
  IF (nbfiles > nbfile_max) THEN
    CALL histerr (3,'flincre', &
     'Too many files. Please increase nbfil_max', &
     'in program flincom.F90.',' ')
  ENDIF
!-
  ncids(nbfiles) = fid
  ncnbd(nbfiles) = 4
!-
  ncdims(nbfiles,1:4) = (/ iim1, jjm1, llm1, ttm1 /)
!-
  ncfunli(nbfiles) = -1
  ncnba(nbfiles)   =  4
  ncnbva(nbfiles)  =  0
  ncfileopen(nbfiles) = .TRUE.
!-
  fid_out = nbfiles
!---------------------
END SUBROUTINE flincre
!-
!===
!-
SUBROUTINE flinopen_zoom2d &
  (filename, iideb, iilen, jjdeb, jjlen, do_test, &
  iim, jjm, llm, lon, lat, lev, ttm, itaus, date0, dt, fid_out)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  LOGICAL :: do_test
  INTEGER :: iim, jjm, llm, ttm, iideb, iilen, jjdeb, jjlen
  REAL :: lon(iilen,jjlen), lat(iilen,jjlen), lev(llm)
  INTEGER :: itaus(ttm)
  REAL :: date0, dt
  INTEGER :: fid_out
!-
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE (*,*) ' iideb, iilen, jjdeb, jjlen, iim, jjm ', &
                           iideb, iilen, jjdeb, jjlen, iim, jjm
  IF (check) WRITE (*,*) ' lon ', lon(1,1), lon(iilen,jjlen)
  IF (check) WRITE (*,*) ' lat ', lat(1,1), lat(iilen,jjlen)
!-
  CALL flinopen_work &
    (filename, iideb, iilen, jjdeb, jjlen, do_test, &
     iim, jjm, llm, lon, lat, lev, ttm, itaus, date0, dt, fid_out)
!-----------------------------
END SUBROUTINE flinopen_zoom2d
!-
!===
!-
SUBROUTINE flinopen_nozoom &
  (filename, do_test, iim, jjm, llm, lon, lat, lev, ttm, &
   itaus, date0, dt, fid_out)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  LOGICAL :: do_test
  INTEGER :: iim, jjm, llm, ttm
  REAL :: lon(iim,jjm), lat(iim,jjm), lev(llm)
  INTEGER :: itaus(ttm)
  REAL :: date0, dt
  INTEGER :: fid_out
!---------------------------------------------------------------------
  CALL flinopen_work &
    (filename, 1, iim, 1, jjm, do_test, &
     iim, jjm, llm, lon, lat, lev, ttm, itaus, date0, dt, fid_out)
!-------------------------
END SUBROUTINE flinopen_nozoom
!-
!===
!-
SUBROUTINE flinopen_work &
  (filename, iideb, iilen, jjdeb, jjlen, do_test, &
   iim, jjm, llm, lon, lat, lev, ttm, itaus, date0, dt, fid_out)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  LOGICAL :: do_test
  INTEGER :: iim, jjm, llm, ttm, iideb, iilen, jjdeb, jjlen
  REAL :: lon(iilen,jjlen), lat(iilen,jjlen), lev(llm)
  INTEGER :: itaus(ttm)
  REAL :: date0, dt
  INTEGER :: fid_out
!-
! LOCAL
!-
  REAL, PARAMETER :: eps = 1.e-4
!-
  INTEGER :: iret, vid, fid, nbdim, i, iilast, jjlast
  INTEGER :: gdtt_id, old_id, iv, gdtmaf_id
  CHARACTER(LEN=250) :: name
  CHARACTER(LEN=80) :: units, calendar
  INTEGER :: tmp_iim, tmp_jjm, tmp_llm, tmp_ttm
  REAL :: x_first, x_last
  INTEGER :: year, month, day
  REAL :: r_year, r_month, r_day
  INTEGER :: year0, month0, day0, hours0, minutes0, seci
  REAL :: sec, sec0
  CHARACTER :: strc
!-
  REAL,DIMENSION(:),ALLOCATABLE :: vec_tmp
!-
  LOGICAL :: open_file
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  iilast = iideb+iilen-1
  jjlast = jjdeb+jjlen-1
  IF (check) WRITE (*,*) &
    ' flinopen_work zoom 2D information '// &
    ' iideb, iilen, iilast, jjdeb, jjlen, jjlast ', &
      iideb, iilen, iilast, jjdeb, jjlen, jjlast
!-
! 1.0 get all infos on the file
!-
! Either the fid_out has not been initialized (0 or very large)
! then we have to open anyway. Else we only need to open the file
! if it has not been opened before.
!-
  IF ( (fid_out < 1).OR.(fid_out > nbfile_max) ) THEN
    open_file = .TRUE.
  ELSE IF (.NOT.ncfileopen(fid_out)) THEN
    open_file = .TRUE.
  ELSE
    open_file = .FALSE.
  ENDIF
!-
  IF (open_file) THEN
    CALL flininfo (filename,tmp_iim,tmp_jjm,tmp_llm,tmp_ttm,fid_out)
  ELSE
!-- The user has already opened the file
!-- and we trust that he knows the dimensions
    tmp_iim = iim
    tmp_jjm = jjm
    tmp_llm = llm
    tmp_ttm = ttm
  ENDIF
!-
  IF (check) &
    WRITE(*,*) 'OUT OF flininfo :',tmp_iim,tmp_jjm,tmp_llm,tmp_ttm
!-
  fid = ncids(fid_out)
!-
! 2.0 get the sizes and names of the different coordinates
!     and do a first set of verification.
!-
! 2.2 We test the axis if we have to.
!-
  IF (check) &
    WRITE(*,*) 'flininfo 2.2 We test if we have to test : ',do_test
!-
  IF (do_test) THEN
    IF      (iim /= tmp_iim) THEN
      CALL histerr (3,'flinopen', &
        'file '//filename//' does not have the ', &
        'required dimension in x direction (longitude)',' ')
    ELSE IF (jjm /= tmp_jjm) THEN
      CALL histerr (3,'flinopen', &
        'file '//filename//' does not have the ', &
        'required dimension in y direction (latitude)',' ')
    ELSE IF ( llm /= tmp_llm .AND. llm > 0 ) THEN
      CALL histerr (3,'flinopen', &
        'file '//filename//' does not have the ', &
        'required dimension in the vertical',' ')
    ENDIF
  ELSE
!---
!-- 2.3 Else the sizes of the axes are returned to the user
!---
    IF (check) WRITE(*,*) 'flinopen 2.3 Else sizes are returned'
!---
    iim = tmp_iim
    jjm = tmp_jjm
    llm = tmp_llm
  ENDIF
!-
  ttm = tmp_ttm
!-
! 3.0 Check if we are realy talking about the same coodinate system
!     if not then we get the lon, lat and lev variables from the file
!-
  IF (check) WRITE(*,*) 'flinopen 3.0 we are realy talking'
!-
  IF (do_test) THEN
!---
    CALL flinfindcood (fid_out, 'lon', vid, nbdim)
    iret = NF90_GET_VAR (fid, vid, x_first, start=(/ iideb, jjdeb /))
    iret = NF90_GET_VAR (fid, vid, x_last, start=(/ iilast, jjlast /))
!---
    IF (check) &
      WRITE(*,*) 'from file lon first and last, modulo 360. ', &
        x_first, x_last, MODULO(x_first,360.), MODULO(x_last,360.)
    IF (check) &
      WRITE(*,*) 'from model lon first and last, modulo 360. ', &
        lon(1,1),lon(iilen,jjlen), &
        MODULO(lon(1,1),360.), MODULO(lon(iilen,jjlen),360.)
    IF (    (ABS( MODULO(x_first,360.) &
                 -MODULO(lon(1,1),360.)) > eps) &
        .OR.(ABS( MODULO(x_last,360.) &
                 -MODULO(lon(iilen ,jjlen),360.)) > eps ) ) THEN
      CALL histerr (3,'flinopen', &
        'file '//filename//' and the model do not', &
        'share the same longitude coordinate', &
        'Obtained by comparing the first and last values ')
    ENDIF
!---
    CALL flinfindcood (fid_out, 'lat', vid, nbdim)
    iret = NF90_GET_VAR (fid, vid, x_first, start=(/ iideb, jjdeb /))
    iret = NF90_GET_VAR (fid, vid, x_last, start=(/ iilast, jjlast /))
!---
    IF (check) WRITE(*,*) &
      'from file lat first and last ',x_first,x_last
    IF (check) WRITE(*,*) &
      'from model lat first and last ',lat(1,1),lat(iilen,jjlen)
!---
    IF (    (ABS(x_first-lat(1,1)) > eps) &
        .OR.(ABS(x_last-lat(iilen,jjlen)) > eps) ) THEN
      CALL histerr (3,'flinopen', &
        'file '//filename//' and the model do not', &
        'share the same latitude coordinate', &
        'Obtained by comparing the first and last values ')
    ENDIF
!---
    IF (llm > 0) THEN
      CALL flinfindcood (fid_out, 'lev', vid, nbdim)
      iret = NF90_GET_VAR (fid, vid, x_first, start=(/ 1 /))
      iret = NF90_GET_VAR (fid, vid, x_last, start=(/ llm /))
!-----
      IF (check) WRITE(*,*) &
        'from file lev first and last ',x_first ,x_last
      IF (check) WRITE(*,*) &
        'from model lev first and last ',lev(1),lev(llm)
!-----
      IF (    (ABS(x_first-lev(1)) > eps) &
          .OR.(ABS(x_last-lev(llm)) > eps) ) THEN
        CALL histerr (3,'flinopen', &
          'file '//filename//' and the model do not', &
          'share the same vertical coordinate', &
          'Obtained by comparing the first and last values')
      ENDIF
    ENDIF
!---
  ELSE
!---
!-- 4.0 extracting the coordinates if we do not check
!---
    IF (check) WRITE(*,*) 'flinopen 4.0 extracting the coordinates'
!---
    CALL flinfindcood (fid_out, 'lon', vid, nbdim)
    IF (nbdim == 2) THEN
      iret = NF90_GET_VAR (fid, vid, lon, &
               start=(/ iideb, jjdeb /), count=(/ iilen, jjlen /))
    ELSE
      ALLOCATE(vec_tmp(iilen))
      iret = NF90_GET_VAR (fid, vid, vec_tmp, &
               start=(/ iideb /), count=(/ iilen /))
      DO i=1,jjlen
        lon(:,i) = vec_tmp(:)
      ENDDO
      DEALLOCATE(vec_tmp)
    ENDIF
!---
    CALL flinfindcood (fid_out, 'lat', vid, nbdim)
    IF (nbdim == 2) THEN
      iret = NF90_GET_VAR (fid, vid, lat, &
               start=(/ iideb, jjdeb /), count=(/ iilen, jjlen /))
    ELSE
      ALLOCATE(vec_tmp(jjlen))
      iret = NF90_GET_VAR (fid, vid, vec_tmp, &
               start=(/ jjdeb /), count=(/ jjlen /))
      DO i=1,iilen
        lat(i,:) = vec_tmp(:)
      ENDDO
      DEALLOCATE(vec_tmp)
    ENDIF
!---
    IF (llm > 0) THEN
      CALL flinfindcood (fid_out, 'lev', vid, nbdim)
      IF (nbdim == 1) THEN
        iret = NF90_GET_VAR (fid, vid, lev, &
                 start=(/ 1 /), count=(/ llm /))
      ELSE
        CALL histerr (3,'flinopen', &
          'Can not handle vertical coordinates that have more',&
          'than 1 dimension',' ')
      ENDIF
    ENDIF
  ENDIF
!-
! 5.0 Get all the details for the time if possible needed
!-
  IF (check) WRITE(*,*) 'flinopen 5.0 Get time'
!-
  IF (ttm > 0) THEN
!---
!-- 5.1 Find the time axis. Prefered method is the 'timestep since'
!---
    gdtmaf_id = -1
    gdtt_id = -1
    old_id = -1
    DO iv=1,ncnbva(fid_out)
      name=''
      iret = NF90_INQUIRE_VARIABLE (fid, iv, name=name)
      units=''
      iret = NF90_GET_ATT (fid, iv, 'units', units)
      IF (INDEX(units,'seconds since') > 0) gdtmaf_id = iv
      IF (INDEX(units,'timesteps since') > 0) gdtt_id = iv
      IF (INDEX(name, 'tstep') > 0) old_id = iv
    ENDDO
!---
    IF (gdtt_id > 0) THEN
      vid = gdtt_id
    ELSE IF (gdtmaf_id > 0) THEN
      vid = gdtmaf_id
    ELSE IF (old_id > 0) THEN
      vid = old_id
    ELSE
      CALL histerr (3, 'flinopen', 'No time axis found',' ',' ')
    ENDIF
!---
    ALLOCATE(vec_tmp(ttm))
    iret = NF90_GET_VAR (fid,vid,vec_tmp,start=(/ 1 /),count=(/ ttm /))
    itaus(1:ttm) = NINT(vec_tmp(1:ttm))
    DEALLOCATE(vec_tmp)
!---
    IF (check) WRITE(*,*) 'flinopen 5.1 Times ',itaus
!---
!-- Getting all the details for the time axis
!---
!-- Find the calendar
    calendar = ''
    iret = NF90_GET_ATT (fid,gdtmaf_id,'calendar',calendar)
    IF (iret == NF90_NOERR) THEN
      CALL ioconf_calendar(calendar)
    ENDIF
!--
    units = ''
    iret = NF90_GET_ATT (fid,vid,'units',units)
    IF (gdtt_id > 0) THEN
      units = units(INDEX(units,'since')+6:LEN_TRIM(units))
      READ (units,'(I4.4,5(a,I2.2))') &
        year0, strc, month0, strc, day0, &
               strc, hours0, strc, minutes0, strc, seci
      sec0 = hours0*3600. + minutes0*60. + seci
      CALL ymds2ju (year0, month0, day0, sec0, date0)
      IF (check) &
        WRITE(*,*) 'flinopen 5.1 gdtt_id year0 ... date0 ', &
                   year0, month0, day0, sec0, date0
!-----
      iret = NF90_GET_ATT (fid, gdtt_id, 'tstep_sec', dt)
    ELSE IF (gdtmaf_id > 0) THEN
      units = units(INDEX(units,'since')+6:LEN_TRIM(units))
      READ (units,'(I4.4,5(a,I2.2))') &
        year0, strc, month0, strc, day0, &
               strc, hours0, strc, minutes0, strc, seci
      sec0 = hours0*3600. + minutes0*60. + seci
      CALL ymds2ju (year0, month0, day0, sec0, date0)
!-----
      IF (check) &
        WRITE(*,*) 'flinopen 5.1 gdtmaf_id year0 ... date0 ', &
                   year0, month0, day0, sec0, date0
    ELSE IF (old_id > 0) THEN
      iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'delta_tstep_sec', dt)
      iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'day0', r_day)
      iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'sec0', sec)
      iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'year0', r_year)
      iret = NF90_GET_ATT (fid, NF90_GLOBAL, 'month0', r_month)
!-----
      day = INT(r_day)
      month = INT(r_month)
      year = INT(r_year)
!-----
      CALL ymds2ju (year, month, day, sec, date0)
    ENDIF
  ENDIF
!-
  IF (check) WRITE(*,*) 'flinopen 6.0 File opened', date0, dt
!---------------------------
END SUBROUTINE flinopen_work
!-
!===
!-
SUBROUTINE flininfo (filename, iim, jjm, llm, ttm, fid_out)
!---------------------------------------------------------------------
!- This subroutine allows to get some information.
!- It is usualy done within flinopen but the user may want to call
!- it before in order to allocate the space needed to extract the
!- data from the file.
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  INTEGER :: iim, jjm, llm, ttm, fid_out
!-
! LOCAL
!-
  INTEGER :: iret, fid, ndims, nvars, nb_atts, id_unlim
  INTEGER :: iv, lll
  INTEGER :: xid, yid, zid, tid
  CHARACTER(LEN=80) :: name
  CHARACTER(LEN=30) :: axname
!-
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  lll = LEN_TRIM(filename)
  IF (filename(lll-2:lll) /= '.nc') THEN
    name = filename(1:lll)//'.nc'
  ELSE
    name = filename(1:lll)
  ENDIF
!-
  iret = NF90_OPEN (name, NF90_NOWRITE, fid)
  IF (iret /= NF90_NOERR) THEN
    CALL histerr(3, 'flininfo','Could not open file :',TRIM(name),' ')
  ENDIF
!-
  iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, &
                      nAttributes=nb_atts, unlimitedDimId=id_unlim)
!-
  xid = -1; iim = 0;
  yid = -1; jjm = 0;
  zid = -1; llm = 0;
  tid = -1; ttm = 0;
!-
  DO iv=1,ndims
!---
    iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll)
    CALL strlowercase (axname)
    axname = ADJUSTL(axname)
!---
    IF (check) WRITE(*,*) &
      'flininfo - getting axname',iv,axname,lll
!---
    IF      (    (INDEX(axname,'x') == 1) &
             .OR.(INDEX(axname,'lon') == 1) ) THEN
      xid = iv; iim = lll;
    ELSE IF (    (INDEX(axname,'y') == 1) &
             .OR.(INDEX(axname,'lat') == 1) ) THEN
      yid = iv; jjm = lll;
    ELSE IF (    (INDEX(axname,'lev') == 1) &
             .OR.(INDEX(axname,'plev') == 1) &
             .OR.(INDEX(axname,'z') == 1) &
             .OR.(INDEX(axname,'depth') == 1) ) THEN
      zid = iv; llm = lll;
    ELSE IF (    (INDEX(axname,'tstep') == 1) &
             .OR.(INDEX(axname,'time_counter') == 1) ) THEN
!---- For the time we certainly need to allow for other names
      tid = iv; ttm = lll;
    ELSE IF (ndims == 1) THEN
!---- Nothing was found and ndims=1 then we have a vector of data
      xid = 1; iim = lll;
    ENDIF
!---
  ENDDO
!-
! Keep all this information
!-
  nbfiles = nbfiles+1
!-
  IF (nbfiles > nbfile_max) THEN
    CALL histerr (3,'flininfo', &
      'Too many files. Please increase nbfil_max', &
      'in program flincom.F90.',' ')
  ENDIF
!-
  ncids(nbfiles) = fid
  ncnbd(nbfiles) = ndims
!-
  ncdims(nbfiles,1:4) = (/ iim, jjm, llm, ttm /)
!-
  ncfunli(nbfiles) = id_unlim
  ncnba(nbfiles)   = nb_atts
  ncnbva(nbfiles)  = nvars
  ncfileopen(nbfiles) = .TRUE.
!-
  fid_out = nbfiles
!----------------------
END SUBROUTINE flininfo
!-
!===
!-
SUBROUTINE flinput_r1d &
  (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  REAL :: var(:)
!-
  INTEGER :: fid, ncvarid, ndim, iret
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE(*,*) &
     "flinput_r1d : SIZE(var) = ",SIZE(var)
!-
  CALL flinput_mat &
    (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid, &
     fid,ncvarid,ndim)
!-
  iret = NF90_PUT_VAR (fid, ncvarid, var, &
           start=w_sta(1:ndim), count=w_len(1:ndim))
!-------------------------
END SUBROUTINE flinput_r1d
!-
!===
!-
SUBROUTINE flinput_r2d &
  (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  REAL :: var(:,:)
!-
  INTEGER :: fid, ncvarid, ndim, iret
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE(*,*) &
     "flinput_r2d : SIZE(var) = ",SIZE(var)
!-
  CALL flinput_mat &
    (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid, &
     fid,ncvarid,ndim)
!-
  iret = NF90_PUT_VAR (fid, ncvarid, var, &
           start=w_sta(1:ndim), count=w_len(1:ndim))
!-------------------------
END SUBROUTINE flinput_r2d
!-
!===
!-
SUBROUTINE flinput_r3d &
  (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  REAL :: var(:,:,:)
!-
  INTEGER :: fid, ncvarid, ndim, iret
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE(*,*) &
     "flinput_r3d : SIZE(var) = ",SIZE(var)
!-
  CALL flinput_mat &
    (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid, &
     fid,ncvarid,ndim)
!-
  iret = NF90_PUT_VAR (fid, ncvarid, var, &
           start=w_sta(1:ndim), count=w_len(1:ndim))
!-------------------------
END SUBROUTINE flinput_r3d
!-
!===
!-
SUBROUTINE flinput_r4d &
  (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  REAL :: var(:,:,:,:)
!-
  INTEGER :: fid, ncvarid, ndim, iret
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE(*,*) &
     "flinput_r4d : SIZE(var) = ",SIZE(var)
!-
  CALL flinput_mat &
    (fid_in,varname,iim,nlonid,jjm,nlatid,llm,zdimid,ttm,tdimid, &
     fid,ncvarid,ndim)
!-
  iret = NF90_PUT_VAR (fid, ncvarid, var, &
           start=w_sta(1:ndim), count=w_len(1:ndim))
!-------------------------
END SUBROUTINE flinput_r4d
!-
!===
!-
SUBROUTINE flinput_mat &
  (fid_in,varname,iim,nlonid,jjm,nlatid, &
                  llm,zdimid,ttm,tdimid,fid,ncvarid,ndim)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  INTEGER :: fid, ncvarid, ndim
!-
! LOCAL
!-
  INTEGER :: iret
!---------------------------------------------------------------------
  fid = ncids(fid_in)
!-
  w_sta(1:4) = (/      1,      1,  1,  1 /)
  w_len(1:2) = (/    iim,    jjm /)
  w_dim(1:2) = (/ nlonid, nlatid /)
!-
  IF ( (llm > 0).AND.(ttm > 0) ) THEN
    ndim = 4
    w_len(3:4) = (/    llm,    ttm /)
    w_dim(3:4) = (/ zdimid, tdimid /)
  ELSE IF (llm > 0) THEN
    ndim = 3
    w_dim(3) = zdimid
    w_len(3) = llm
  ELSE IF (ttm > 0) THEN
    ndim = 3
    w_dim(3) = tdimid
    w_len(3) = ttm
  ELSE
    ndim = 2
  ENDIF
!-
  iret = NF90_REDEF   (fid)
  iret = NF90_DEF_VAR (fid,varname,NF90_FLOAT,w_dim(1:ndim),ncvarid)
  iret = NF90_PUT_ATT (fid,ncvarid,'short_name',TRIM(varname))
  iret = NF90_ENDDEF  (fid)
!--------------------------
END  SUBROUTINE flinput_mat
!-
!===
!-
SUBROUTINE flinput_scal &
  (fid_in, varname, iim, nlonid, jjm, nlatid, &
                    llm, zdimid, ttm, tdimid, var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, nlonid, jjm, nlatid, llm, zdimid, ttm, tdimid
  REAL :: var
!-
! LOCAL
!-
  INTEGER :: fid, iret
!---------------------------------------------------------------------
  fid = ncids(fid_in)
!-
  iret = NF90_REDEF   (fid)
  iret = NF90_PUT_ATT (fid, NF90_GLOBAL, varname, REAL(var,KIND=4))
  iret = NF90_ENDDEF  (fid)
!---------------------------
END  SUBROUTINE flinput_scal
!-
!===
!-
SUBROUTINE flinget_r1d &
  (fid_in,varname,iim,jjm,llm,ttm,itau_dep,itau_fin,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm, itau_dep, itau_fin
  REAL :: var(:)
!-
  INTEGER :: jl, ji
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r1d : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r1d : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,1,iim,1,jjm,buff_tmp)
!-
  jl=0
  DO ji=1,SIZE(var,1)
    jl=jl+1
    var(ji) = buff_tmp(jl)
  ENDDO
!-------------------------
END SUBROUTINE flinget_r1d
!-
!===
!-
SUBROUTINE flinget_r2d &
  (fid_in,varname,iim,jjm,llm,ttm,itau_dep,itau_fin,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm, itau_dep, itau_fin
  REAL :: var(:,:)
!-
  INTEGER :: jl, jj, ji
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r2d : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r2d : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,1,iim,1,jjm,buff_tmp)
!-
  jl=0
  DO jj=1,SIZE(var,2)
    DO ji=1,SIZE(var,1)
      jl=jl+1
      var(ji,jj) = buff_tmp(jl)
    ENDDO
  ENDDO
!-------------------------
END SUBROUTINE flinget_r2d
!-
!===
!-
SUBROUTINE flinget_r2d_zoom2d &
  (fid_in,varname,iim,jjm,llm,ttm, &
   itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim,jjm,llm,ttm,itau_dep,itau_fin,iideb,jjdeb,iilen,jjlen
  REAL :: var(:,:)
!-
  INTEGER :: jl, jj, ji
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r2d_zoom : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r2d_zoom : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,buff_tmp)
!-
  jl=0
  DO jj=1,SIZE(var,2)
    DO ji=1,SIZE(var,1)
      jl=jl+1
      var(ji,jj) = buff_tmp(jl)
    ENDDO
  ENDDO
!--------------------------------
END SUBROUTINE flinget_r2d_zoom2d
!-
!===
!-
SUBROUTINE flinget_r3d &
  (fid_in,varname,iim,jjm,llm,ttm,itau_dep,itau_fin,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm, itau_dep, itau_fin
  REAL :: var(:,:,:)
!-
  INTEGER :: jl, jk, jj, ji
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r3d : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r3d : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,1,iim,1,jjm,buff_tmp)
!-
  jl=0
  DO jk=1,SIZE(var,3)
    DO jj=1,SIZE(var,2)
      DO ji=1,SIZE(var,1)
        jl=jl+1
        var(ji,jj,jk) = buff_tmp(jl)
      ENDDO
    ENDDO
  ENDDO
!-------------------------
END SUBROUTINE flinget_r3d
!-
!===
!-
SUBROUTINE flinget_r3d_zoom2d &
  (fid_in,varname,iim,jjm,llm,ttm, &
   itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim,jjm,llm,ttm,itau_dep,itau_fin,iideb,jjdeb,iilen,jjlen
  REAL :: var(:,:,:)
!-
  INTEGER :: jl, jk, jj, ji
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r3d_zoom : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r3d_zoom : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,buff_tmp)
!-
  jl=0
  DO jk=1,SIZE(var,3)
    DO jj=1,SIZE(var,2)
      DO ji=1,SIZE(var,1)
        jl=jl+1
        var(ji,jj,jk) = buff_tmp(jl)
      ENDDO
    ENDDO
  ENDDO
!--------------------------------
END SUBROUTINE flinget_r3d_zoom2d
!-
!===
!-
SUBROUTINE flinget_r4d &
  (fid_in,varname,iim,jjm,llm,ttm,itau_dep,itau_fin,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm, itau_dep, itau_fin
  REAL :: var(:,:,:,:)
!-
  INTEGER :: jl, jk, jj, ji, jm
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r4d : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r4d : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,1,iim,1,jjm,buff_tmp)
!-
  jl=0
  DO jm=1,SIZE(var,4)
    DO jk=1,SIZE(var,3)
      DO jj=1,SIZE(var,2)
        DO ji=1,SIZE(var,1)
          jl=jl+1
          var(ji,jj,jk,jm) = buff_tmp(jl)
        ENDDO
      ENDDO
    ENDDO
  ENDDO
!-------------------------
END SUBROUTINE flinget_r4d
!-
!===
!-
SUBROUTINE flinget_r4d_zoom2d &
  (fid_in,varname,iim,jjm,llm,ttm, &
   itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,var)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim,jjm,llm,ttm,itau_dep,itau_fin,iideb,jjdeb,iilen,jjlen
  REAL :: var(:,:,:,:)
!-
  INTEGER :: jl, jk, jj, ji, jm
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: buff_tmp
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (.NOT.ALLOCATED(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r4d_zoom : allocate buff_tmp for buff_sz = ",SIZE(var)
    ALLOCATE (buff_tmp(SIZE(var)))
  ELSE IF (SIZE(var) > SIZE(buff_tmp)) THEN
    IF (check) WRITE(*,*) &
      "flinget_r4d_zoom : re-allocate buff_tmp for buff_sz = ",SIZE(var)
    DEALLOCATE (buff_tmp)
    ALLOCATE (buff_tmp(SIZE(var)))
  ENDIF
!-
  CALL flinget_mat (fid_in,varname,iim,jjm,llm,ttm, &
                    itau_dep,itau_fin,iideb,iilen,jjdeb,jjlen,buff_tmp)
!-
  jl=0
  DO jm=1,SIZE(var,4)
    DO jk=1,SIZE(var,3)
      DO jj=1,SIZE(var,2)
        DO ji=1,SIZE(var,1)
          jl=jl+1
          var(ji,jj,jk,jm) = buff_tmp(jl)
        ENDDO
      ENDDO
    ENDDO
  ENDDO
!--------------------------------
END SUBROUTINE flinget_r4d_zoom2d
!-
!===
!-
SUBROUTINE flinget_mat &
  (fid_in, varname, iim, jjm, llm, ttm, itau_dep, &
   itau_fin, iideb, iilen, jjdeb, jjlen, var)
!---------------------------------------------------------------------
!- This subroutine will read the variable named varname from
!- the file previously opened by flinopen and identified by fid
!-
!- It is checked that the dimensions of the variable to be read
!- correspond to what the user requested when he specified
!- iim, jjm and llm. The only exception which is allowed is
!- for compressed data where the horizontal grid is not expected
!- to be iim x jjm.
!-
!- If variable is of size zero a global attribute is read.
!- This global attribute will be of type real
!-
!- INPUT
!-
!- fid      : File ID returned by flinopen
!- varname  : Name of the variable to be read from the file
!- iim      : | These three variables give the size of the variables
!- jjm      : | to be read. It will be verified that the variables
!- llm      : | fits in there.
!- ttm      : |
!- itau_dep : Time step at which we will start to read
!- itau_fin : Time step until which we are going to read
!-            For the moment this is done on indexes
!-            but it should be in the physical space.
!-            If there is no time-axis in the file then use a
!-            itau_fin < itau_dep, this will tell flinget not to
!-            expect a time-axis in the file.
!- iideb    : index i for zoom
!- iilen    : length of zoom
!- jjdeb    : index j for zoom
!- jjlen    : length of zoom
!-
!- OUTPUT
!-
!- var      : array that will contain the data
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm
  INTEGER :: itau_dep, itau_fin, iideb, iilen, jjdeb, jjlen
  REAL :: var(:)
!-
! LOCAL
!-
  INTEGER :: iret, fid
  INTEGER :: vid, cvid, clen
  CHARACTER(LEN=70) :: str1
  CHARACTER(LEN=250) :: att_n, tmp_n
  CHARACTER(LEN=5) :: axs_l
  INTEGER :: tmp_i
  REAL,SAVE :: mis_v=0.
  REAL :: tmp_r
  INTEGER :: ndims, x_typ, nb_atts
  INTEGER,DIMENSION(NF90_MAX_VAR_DIMS) :: dimids
  INTEGER :: i, nvars, i2d, cnd
  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: var_tmp
  LOGICAL :: uncompress = .FALSE.
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  fid = ncids(fid_in)
!-
  IF (check) THEN
    WRITE(*,*) &
    'flinget_mat : fid_in, fid, varname :', fid_in, fid, TRIM(varname)
    WRITE(*,*) &
    'flinget_mat : iim, jjm, llm, ttm, itau_dep, itau_fin :', &
    iim, jjm, llm, ttm, itau_dep, itau_fin
    WRITE(*,*) &
    'flinget_mat : iideb, iilen, jjdeb, jjlen :', &
    iideb, iilen, jjdeb, jjlen
  ENDIF
!-
  uncompress = .FALSE.
!-
! 1.0 We get first all the details on this variable from the file
!-
  nvars = ncnbva(fid_in)
!-
  vid = -1
  iret = NF90_INQ_VARID (fid, varname, vid)
!-
  IF (vid < 0 .OR. iret /= NF90_NOERR) THEN
    CALL histerr (3,'flinget', &
      'Variable '//TRIM(varname)//' not found in file',' ',' ')
  ENDIF
!-
  iret = NF90_INQUIRE_VARIABLE (fid, vid, &
           ndims=ndims, dimids=dimids, nAtts=nb_atts)
  IF (check) THEN
    WRITE(*,*) &
    'flinget_mat : fid, vid :', fid, vid
    WRITE(*,*) &
    'flinget_mat : ndims, dimids(1:ndims), nb_atts :', &
    ndims, dimids(1:ndims), nb_atts
  ENDIF
!-
  w_dim(:) = 0
  DO i=1,ndims
    iret  = NF90_INQUIRE_DIMENSION (fid, dimids(i), len=w_dim(i))
  ENDDO
  IF (check) WRITE(*,*) &
    'flinget_mat : w_dim :', w_dim(1:ndims)
!-
  mis_v = 0.0; axs_l = ' ';
!-
  IF (nb_atts > 0) THEN
    IF (check) THEN
      WRITE(*,*) 'flinget_mat : attributes for variable :'
    ENDIF
  ENDIF
  DO i=1,nb_atts
    iret = NF90_INQ_ATTNAME (fid, vid, i, att_n)
    iret = NF90_INQUIRE_ATTRIBUTE (fid, vid, att_n, xtype=x_typ)
    CALL strlowercase (att_n)
    IF      (    (x_typ == NF90_INT).OR.(x_typ == NF90_SHORT) &
             .OR.(x_typ == NF90_BYTE) ) THEN
      iret = NF90_GET_ATT (fid, vid, att_n, tmp_i)
      IF (check) THEN
        WRITE(*,*) '   ',TRIM(att_n),' : ',tmp_i
      ENDIF
    ELSE IF ( (x_typ == NF90_FLOAT).OR.(x_typ == NF90_DOUBLE) ) THEN
      iret = NF90_GET_ATT (fid, vid, att_n, tmp_r)
      IF (check) THEN
        WRITE(*,*) '   ',TRIM(att_n),' : ',tmp_r
      ENDIF
      IF (index(att_n,'missing_value') > 0) THEN
        mis_v = tmp_r
      ENDIF
    ELSE
      tmp_n = ''
      iret = NF90_GET_ATT (fid, vid, att_n, tmp_n)
      IF (check) THEN
        WRITE(*,*) '   ',TRIM(att_n),' : ',TRIM(tmp_n)
      ENDIF
      IF (index(att_n,'axis') > 0) THEN
        axs_l = tmp_n
      ENDIF
    ENDIF
  ENDDO
!?
!!!!!!!!!! We will need a verification on the type of the variable
!?
!-
! 2.0 The dimensions are analysed to determine what is to be read
!-
! 2.1 the longitudes
!-
  IF ( w_dim(1) /= iim .OR. w_dim(2) /= jjm) THEN
!---
!-- There is a possibility that we have to deal with a compressed axis !
!---
    iret = NF90_INQUIRE_DIMENSION (fid, dimids(1), &
             name=tmp_n, len=clen)
    iret = NF90_INQ_VARID (fid, tmp_n, cvid)
!---
    IF (check) WRITE(*,*) &
      'Dimname, iret , NF90_NOERR : ',TRIM(tmp_n),iret,NF90_NOERR
!---
!-- If we have an axis which has the same name
!-- as the dimension we can see if it is compressed
!---
!-- TODO TODO for zoom2d
!---
    IF (iret == NF90_NOERR) THEN
      iret = NF90_GET_ATT (fid, cvid, 'compress', str1)
!-----
      IF (iret == NF90_NOERR) THEN
        iret = NF90_INQUIRE_VARIABLE (fid,cvid,xtype=x_typ,ndims=cnd)
!-------
        IF ( cnd /= 1 .AND. x_typ /= NF90_INT) THEN
          CALL histerr (3,'flinget', &
            'Variable '//TRIM(tmp_n)//' can not be a compressed axis', &
            'Either it has too many dimensions'// &
            ' or it is not of type integer', ' ')
        ELSE
!---------
!-------- Let us see if we already have that index table
!---------
          IF (    (cind_len /= clen).OR.(cind_vid /= cvid) &
              .OR.(cind_fid /= fid) ) THEN
            IF (ALLOCATED(cindex))   DEALLOCATE(cindex)
            ALLOCATE(cindex(clen))
            cind_len = clen
            cind_vid = cvid
            cind_fid = fid
            iret = NF90_GET_VAR (fid, cvid, cindex)
          ENDIF
!---------
!-------- In any case we need to set the slab of data to be read
!---------
          uncompress = .TRUE.
          w_sta(1) = 1
          w_len(1) = clen
          i2d = 1
        ENDIF
      ELSE
        str1 = 'The horizontal dimensions of '//varname
        CALL histerr (3,'flinget',str1, &
          'is not compressed and does not'// &
          ' correspond to the requested size',' ')
      ENDIF
    ELSE
      IF (w_dim(1) /= iim) THEN
        str1 = 'The longitude dimension of '//varname
        CALL histerr (3,'flinget',str1, &
          'in the file is not equal to the dimension', &
          'that should be read')
      ENDIF
      IF (w_dim(2) /= jjm) THEN
        str1 = 'The latitude dimension of '//varname
        CALL histerr (3,'flinget',str1, &
          'in the file is not equal to the dimension', &
          'that should be read')
      ENDIF
    ENDIF
  ELSE
    w_sta(1:2) = (/ iideb, jjdeb /)
    w_len(1:2) = (/ iilen, jjlen /)
    i2d = 2
  ENDIF
!-
! 2.3 Now the difficult part, the 3rd dimension which can be
! time or levels.
!-
! Priority is given to the time axis if only three axes are present.
!-
  IF (ndims > i2d) THEN
!---
!-- 2.3.1 We have a vertical axis
!---
    IF (llm == 1 .AND. ndims == i2d+2 .OR. llm == w_dim(i2d+1)) THEN
!-----
      IF (w_dim(i2d+1) /= llm) THEN
        CALL histerr (3,'flinget', &
          'The vertical dimension of '//varname, &
          'in the file is not equal to the dimension', &
          'that should be read')
      ELSE
        w_sta(i2d+1) = 1
        IF (llm > 0) THEN
          w_len(i2d+1) = llm
        ELSE
          w_len(i2d+1) = w_sta(i2d+1)
        ENDIF
      ENDIF
!-----
      IF ((itau_fin-itau_dep) >= 0) THEN
        IF      (ndims /= i2d+2) THEN
          CALL histerr (3,'flinget', &
            'You attempt to read a time slab', &
            'but there is no time axis on this variable', varname)
        ELSE IF ((itau_fin - itau_dep) <= w_dim(i2d+2)) THEN
          w_sta(i2d+2) = itau_dep
          w_len(i2d+2) = itau_fin-itau_dep+1
        ELSE
          CALL histerr (3,'flinget', &
            'The time step you try to read is not', &
            'in the file (1)', varname)
        ENDIF
      ELSE IF (ndims == i2d+2 .AND. w_dim(i2d+2) > 1) THEN
        CALL histerr (3,'flinget', &
          'There is a time axis in the file but no', &
          'time step give in the call', varname)
      ELSE
        w_sta(i2d+2) = 1
        w_len(i2d+2) = 1
      ENDIF
    ELSE
!-----
!---- 2.3.2 We do not have any vertical axis
!-----
      IF (ndims == i2d+2) THEN
        CALL histerr (3,'flinget', &
          'The file contains 4 dimensions', &
          'but only 3 are requestes for variable ', varname)
      ENDIF
      IF ((itau_fin-itau_dep) >= 0) THEN
        IF (ndims == i2d+1) THEN
          IF ((itau_fin-itau_dep) < w_dim(i2d+1) ) THEN
            w_sta(i2d+1) = itau_dep
            w_len(i2d+1) = itau_fin-itau_dep+1
          ELSE
            CALL histerr (3,'flinget', &
              'The time step you try to read is not', &
              'in the file (2)', varname)
          ENDIF
        ELSE
          CALL histerr (3,'flinget', &
            'From your input you sould have 3 dimensions', &
            'in the file but there are 4', varname)
        ENDIF
      ELSE
        IF (ndims == i2d+1 .AND. w_dim(i2d+1) > 1) THEN
          CALL histerr (3,'flinget', &
            'There is a time axis in the file but no', &
            'time step given in the call', varname)
        ELSE
          w_sta(i2d+1) = 1
          w_len(i2d+1) = 1
        ENDIF
      ENDIF
    ENDIF
  ELSE
!---
!-- 2.3.3 We do not have any vertical axis
!---
    w_sta(i2d+1:i2d+2) = (/ 0, 0 /)
    w_len(i2d+1:i2d+2) = (/ 0, 0 /)
  ENDIF
!-
! 3.0 Reading the data
!-
  IF (check) WRITE(*,*) &
    'flinget_mat 3.0 : ', uncompress, w_sta, w_len
!---
  IF (uncompress) THEN
!---
    IF (ALLOCATED(var_tmp)) THEN
      IF (SIZE(var_tmp) < clen) THEN
        DEALLOCATE(var_tmp)
        ALLOCATE(var_tmp(clen))
      ENDIF
    ELSE
      ALLOCATE(var_tmp(clen))
    ENDIF
!---
    iret = NF90_GET_VAR (fid, vid, var_tmp, &
             start=w_sta(:), count=w_len(:))
!---
    var(:) = mis_v
    var(cindex(:)) = var_tmp(:)
!---
  ELSE
    iret = NF90_GET_VAR (fid, vid, var, &
             start=w_sta(:), count=w_len(:))
  ENDIF
!-
  IF (check) WRITE(*,*) 'flinget_mat 3.1 : ',NF90_STRERROR (iret)
!--------------------------
END  SUBROUTINE flinget_mat
!-
!===
!-
SUBROUTINE flinget_scal &
  (fid_in, varname, iim, jjm, llm, ttm, itau_dep, itau_fin, var)
!---------------------------------------------------------------------
!- This subroutine will read the variable named varname from
!- the file previously opened by flinopen and identified by fid
!-
!- If variable is of size zero a global attribute is read. This
!- global attribute will be of type real
!-
!- INPUT
!-
!- fid      : File ID returned by flinopen
!- varname  : Name of the variable to be read from the file
!- iim      : | These three variables give the size of the variables
!- jjm      : | to be read. It will be verified that the variables
!- llm      : | fits in there.
!- ttm      : |
!- itau_dep : Time step at which we will start to read
!- itau_fin : Time step until which we are going to read
!-           For the moment this is done on indeces but it should be
!-           in the physical space
!-           If there is no time-axis in the file then use a
!-           itau_fin < itau_dep, this will tell flinget not to
!-           expect a time-axis in the file.
!-
!- OUTPUT
!-
!- var      : scalar that will contain the data
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) :: varname
  INTEGER :: iim, jjm, llm, ttm, itau_dep, itau_fin
  REAL :: var
!-
! LOCAL
!-
  INTEGER :: iret, fid
!-
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) THEN
    WRITE (*,*) 'flinget_scal in file with id ',fid_in
  ENDIF
!-
  fid = ncids(fid_in)
!-
! 1.0 Reading a global attribute
!-
  iret = NF90_GET_ATT (fid, NF90_GLOBAL, varname, var)
!---------------------------
END  SUBROUTINE flinget_scal
!-
!===
!-
SUBROUTINE flinfindcood (fid_in, axtype, vid, ndim)
!---------------------------------------------------------------------
!- This subroutine explores the file in order to find
!- the coordinate according to a number of rules
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! ARGUMENTS
!-
  INTEGER :: fid_in, vid, ndim
  CHARACTER(LEN=3) :: axtype
!-
! LOCAL
!-
  INTEGER :: iv, iret, dimnb
  CHARACTER(LEN=40) :: dimname, dimuni1, dimuni2, dimuni3
  CHARACTER(LEN=80) :: str1
  LOGICAL :: found_rule = .FALSE.
!---------------------------------------------------------------------
  vid = -1
!-
! Make sure all strings are invalid
!-
  dimname = '?-?'
  dimuni1 = '?-?'
  dimuni2 = '?-?'
  dimuni3 = '?-?'
!-
! First rule : we look for the correct units
! lon : east
! lat : north
! We make an exact check as it would be too easy to mistake
! some units by just comparing the substrings.
!-
  SELECTCASE(axtype)
  CASE ('lon')
    dimuni1 = 'degree_e'
    dimuni2 = 'degrees_e'
    found_rule = .TRUE.
  CASE('lat')
    dimuni1 = 'degree_n'
    dimuni2 = 'degrees_n'
    found_rule = .TRUE.
  CASE('lev')
    dimuni1 = 'm'
    dimuni2 = 'km'
    dimuni3 = 'hpa'
    found_rule = .TRUE.
  CASE DEFAULT
    found_rule = .FALSE.
  END SELECT
!-
  IF (found_rule) THEN
    iv = 0
    DO WHILE ( (vid < 0).AND.(iv < ncnbva(fid_in)) )
      iv = iv+1
      str1 = ''
      iret = NF90_GET_ATT (ncids(fid_in), iv, 'units', str1)
      IF (iret == NF90_NOERR) THEN
        CALL strlowercase (str1)
        IF (    (INDEX(str1, TRIM(dimuni1)) == 1) &
            .OR.(INDEX(str1, TRIM(dimuni2)) == 1) &
            .OR.(INDEX(str1, TRIM(dimuni3)) == 1) ) THEN
          vid = iv
          iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, ndims=ndim)
        ENDIF
      ENDIF
    ENDDO
  ENDIF
!-
! Second rule : we find specific names :
! lon : nav_lon
! lat : nav_lat
! Here we can check if we find the substring as the
! names are more specific.
!-
  SELECTCASE(axtype)
  CASE ('lon')
    dimname = 'nav_lon lon longitude'
    found_rule = .TRUE.
  CASE('lat')
    dimname = 'nav_lat lat latitude'
    found_rule = .TRUE.
  CASE('lev')
    dimname = 'plev level depth deptht'
    found_rule = .TRUE.
  CASE DEFAULT
    found_rule = .FALSE.
  END SELECT
!-
  IF (found_rule) THEN
    iv = 0
    DO WHILE ( (vid < 0).AND.(iv < ncnbva(fid_in)) )
      iv = iv+1
      str1=''
      iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, &
                                    name=str1, ndims=ndim)
      IF (INDEX(dimname,TRIM(str1)) >= 1) THEN
        vid = iv
      ENDIF
    ENDDO
  ENDIF
!-
! Third rule : we find a variable with the same name as the dimension
! lon = 1
! lat = 2
! lev = 3
!-
  IF (vid < 0) THEN
    SELECTCASE(axtype)
    CASE ('lon')
      dimnb = 1
      found_rule = .TRUE.
    CASE('lat')
      dimnb = 2
      found_rule = .TRUE.
    CASE('lev')
      dimnb = 3
      found_rule = .TRUE.
    CASE DEFAULT
      found_rule = .FALSE.
    END SELECT
!---
    IF (found_rule) THEN
      iret = NF90_INQUIRE_DIMENSION (ncids(fid_in), dimnb, name=dimname)
      iv = 0
      DO WHILE ( (vid < 0).AND.(iv < ncnbva(fid_in)) )
        iv = iv+1
        str1=''
        iret = NF90_INQUIRE_VARIABLE (ncids(fid_in), iv, &
                                      name=str1, ndims=ndim)
        IF (INDEX(dimname,TRIM(str1)) == 1) THEN
          vid = iv
        ENDIF
      ENDDO
    ENDIF
  ENDIF
!-
! Stop the program if no coordinate was found
!-
  IF (vid < 0) THEN
    CALL histerr (3,'flinfindcood', &
           'No coordinate axis was found in the file', &
           'The data in this file can not be used', axtype)
  ENDIF
!--------------------------
END SUBROUTINE flinfindcood
!-
!===
!-
SUBROUTINE flinclo (fid_in)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
!-
  INTEGER :: iret
!---------------------------------------------------------------------
  iret = NF90_CLOSE (ncids(fid_in))
  ncfileopen(fid_in) = .FALSE.
!---------------------
END SUBROUTINE flinclo
!-
!===
!-
SUBROUTINE flinquery_var(fid_in, varname, exists)
!---------------------------------------------------------------------
!- Queries the existance of a variable in the file.
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: fid_in
  CHARACTER(LEN=*) varname
  LOGICAL :: exists
!-
  INTEGER :: iret, fid, vid
!---------------------------------------------------------------------
  fid = ncids(fid_in)
  vid = -1
  iret = NF90_INQ_VARID (fid, varname, vid)
!-
  exists = ( (vid >= 0).AND.(iret == NF90_NOERR) )
!---------------------------
END SUBROUTINE flinquery_var
!-
!===
!-
SUBROUTINE flininspect (fid_in)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
! fid : File id to inspect
!-
  INTEGER :: fid_in
!-
!- LOCAL
!-
  INTEGER :: iim, jjm, llm, ttm, fid_out
  INTEGER :: iret, fid, ndims, nvars, nb_atts, id_unlim
  INTEGER :: iv, in, lll
  INTEGER :: xid, yid, zid, tid
  INTEGER,DIMENSION(NF90_MAX_VAR_DIMS) :: idimid
  CHARACTER(LEN=80) :: name
  CHARACTER(LEN=30) :: axname
!---------------------------------------------------------------------
  fid = ncids(fid_in)
!-
  iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, &
                       nAttributes=nb_atts, unlimitedDimId=id_unlim)
!-
  WRITE (*,*) 'IOIPSL ID                   : ',fid_in
  WRITE (*,*) 'NetCDF ID                   : ',fid
  WRITE (*,*) 'Number of dimensions        : ',ndims
  WRITE (*,*) 'Number of variables         : ',nvars
  WRITE (*,*) 'Number of global attributes : ',nb_atts
  WRITE (*,*) 'ID unlimited                : ',id_unlim
!-
  xid = -1; iim = 0;
  yid = -1; jjm = 0;
  zid = -1; llm = 0;
  tid = -1; ttm = 0;
!-
  DO iv=1,ndims
!---
    iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll)
    CALL strlowercase (axname)
    axname = ADJUSTL(axname)
!---
    WRITE (*,*) 'Dimension number : ',iv
    WRITE (*,*) 'Dimension name   : ',TRIM(axname)
!---
    IF      (    (INDEX(axname,'x') == 1) &
             .OR.(INDEX(axname,'lon') == 1)) THEN
      xid = iv; iim = lll;
      WRITE (*,*) 'Dimension X size   : ',iim
    ELSE IF (    (INDEX(axname,'y') == 1) &
             .OR.(INDEX(axname,'lat') == 1)) THEN
      yid = iv; jjm = lll;
      WRITE (*,*) 'Dimension Y size   : ',jjm
    ELSE IF (    (INDEX(axname,'lev') == 1) &
             .OR.(INDEX(axname,'plev') == 1) &
             .OR.(INDEX(axname,'z') == 1) &
             .OR.(INDEX(axname,'depth') == 1)) THEN
      zid = iv; llm = lll;
      WRITE (*,*) 'Dimension Z size   : ',llm
    ELSE IF (    (INDEX(axname,'tstep') == 1) &
             .OR.(INDEX(axname,'time_counter') == 1)) THEN
!---- For the time we certainly need to allow for other names
      tid = iv; ttm = lll;
    ELSE IF (ndims == 1) THEN
!---- Nothing was found and ndims=1 then we have a vector of data
      xid = 1; iim = lll;
    ENDIF
!---
  ENDDO
!-
! Keep all this information
!-
  nbfiles = nbfiles+1
!-
  IF (nbfiles > nbfile_max) THEN
    CALL histerr(3,'flininspect', &
      'Too many files. Please increase nbfil_max', &
      'in program flincom.F90.',' ')
  ENDIF
!-
  ncids(nbfiles) = fid
  ncnbd(nbfiles) = ndims
!-
  ncdims(nbfiles,1:4) = (/ iim, jjm, llm, ttm /)
!-
  ncfunli(nbfiles) = id_unlim
  ncnba(nbfiles)   = nb_atts
  ncnbva(nbfiles)  = nvars
  ncfileopen(nbfiles) = .TRUE.
!-
  fid_out = nbfiles
!-
  DO in=1,nvars
    iret = NF90_INQUIRE_VARIABLE (fid, in, &
             name=name, ndims=ndims, dimids=idimid, nAtts=nb_atts)
    WRITE (*,*) 'Variable number  ------------ > ', in
    WRITE (*,*) 'Variable name        : ', TRIM(name)
    WRITE (*,*) 'Number of dimensions : ', ndims
    WRITE (*,*) 'Dimensions ID''s     : ', idimid(1:ndims)
    WRITE (*,*) 'Number of attributes : ', nb_atts
  ENDDO
!-------------------------
END SUBROUTINE flininspect
!-
!===
!-
END MODULE flincom


text:
! POP_GridTypeMappings.f90
!
! Author : Joseph Schoonover
! E-mail : jschoonover@lanl.gov, schoonover.numerics@gmail.com
!
! Copyright 2017 Joseph Schoonover, Los Alamos National Laboratory
!
! Redistribution and use in source and binary forms, with or without
! modification,
! are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors
! may be used to endorse or promote products derived from this
!  software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY,  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
! BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! ////////////////////////////////////////////////////////////////////////////////////////////////

MODULE POP_GridTypeMappings

! src/common/
USE ConstantsDictionary

IMPLICIT NONE

CONTAINS

 SUBROUTINE GetTrueIJ( meshtype, i, j, nX, nY, true_i, true_j )

   IMPLICIT NONE
   INTEGER, INTENT(in)    :: meshtype, i, j, nX, nY
   INTEGER, INTENT(inout) :: true_i, true_j

      IF( meshtype == PeriodicTripole )THEN

         CALL GetTrueIJ_PeriodicTripole( i, j, nX, nY, true_i, true_j )

      ELSE

         PRINT*, 'Module POP_GridTypeMappings : S/R GetTrueIJ : Unknown Grid Type'
         STOP

      ENDIF

 END SUBROUTINE GetTrueIJ
!
 SUBROUTINE GetTrueIJ_PeriodicTripole( i, j, nX, nY, true_i, true_j )
   IMPLICIT NONE
   INTEGER, INTENT(in)  :: i, j, nX, nY
   INTEGER, INTENT(out) :: true_i, true_j
   ! Local
   INTEGER :: dj

      IF( j > nY )THEN

         dj = j - nY
         true_j = nY - (dj-1)
         IF( i == 0 .OR. i == nX+1 )THEN
            true_i = i
         ELSE
            true_i = nX-(i-1)
         ENDIF

      ELSE

         true_j = j

         ! Taking care of the x-periodicity in the map
         IF( i < 1 )THEN
            true_i = nX + i
         ELSEIF( i > nX )THEN
            true_i = i-nX
         ELSE
            true_i = i
         ENDIF

      ENDIF

 END SUBROUTINE GetTrueIJ_PeriodicTripole

END MODULE POP_GridTypeMappings

max_stars_repo_name: lanl/FEOTS

max_issues_repo_name: lanl/FEOTS

max_forks_repo_name: lanl/FEOTS


text:
      SUBROUTINE KSINT(I1)
*
*
*       Regularized integration.
*       ------------------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
      COMMON/MODES/  EB0(NTMAX),ZJ0(NTMAX),ECRIT(NTMAX),AR(NTMAX),
     &               BR(NTMAX),EOSC(4,NTMAX),EDEC(NTMAX),TOSC(NTMAX),
     &               ZP(NTMAX),ES(NTMAX),CZ(2,NTMAX),IOSC(NTMAX),
     &               NAMEC(NTMAX)
      COMMON/SLOW0/  RANGE,ISLOW(10)
      COMMON/GAMDOT/  DGAM
      REAL*8  UI(4),UIDOT(4),XI(6),VI(6),FP(6),FD(6)
      LOGICAL IQ
*
*
*       Set second component, pair index & c.m. index.
      I2 = I1 + 1
      IPAIR = KVEC(I1)
      I = N + IPAIR
*
*       Define perturber membership & inverse c.m. mass.
      NNB0 = LIST(1,I1)
      BODYIN = 1.0/BODY(I)
*
*       Check for further unperturbed motion.
      IF (NNB0.EQ.0.AND.H(IPAIR).LT.0.0) THEN
*       Include possible eccentricity modulation of hierarchical binary.
          IF (NAME(I).LT.0) THEN
              IM = 0
              DO 1 K = 1,NMERGE
                  IF (NAMEM(K).EQ.NAME(I)) IM = K
    1         CONTINUE
              IF (IM.GT.0.AND.TIME.GT.TMDIS(IM)) THEN
                  SEMI = -0.5*BODY(I)/H(IPAIR)
                  ECC2 = (1.0 - R(IPAIR)/SEMI)**2 +
     &                                   TDOT2(IPAIR)**2/(BODY(I)*SEMI)
                  RP = SEMI*(1.0 - SQRT(ECC2))*(1.0 - 2.0*GAMMA(IPAIR))
                  IF (RP.LT.R0(IPAIR)) THEN
                      if(rank.eq.0)
     &                WRITE (77,4)  NAME(I1), SQRT(ECC2), RP, R0(IPAIR)
    4                 FORMAT (' TMDIS TERM    NM E RP R0',I7,1P,3E10.2)
                      CALL FLUSH(77)
                      INSTAB = INSTAB + 1
                      GO TO 90
                  ELSE IF (KZ(27).EQ.2) THEN
                      CALL ECCMOD(I,ITERM)
*       Update time on termination to prevent looping.
                      IF (ITERM.GT.0) THEN
                          T0(I1) = TIME
                          GO TO 90
                      END IF
                  END IF
*       Check any inner and outer circularizing binary using NAMEM & NAMEG.
                  DO 3 K = 1,NCHAOS
                      IF (NAMEC(K).EQ.NZERO - NAMEM(IM)) THEN
*       Update unperturbed binary if T - TOSC > 10 Myr (cf. IMPACT & DECIDE).
                          IF ((TIME - TOSC(K))*TSTAR.GT.10.0) THEN
                              T0(I1) = TIME
                              GO TO 90
                          END IF
                      END IF
                      IF (NAMEC(K).EQ.NAMEG(IM)) THEN
                          IF ((TIME - TOSC(K))*TSTAR.GT.10.0) THEN
                              T0(I1) = TIME
                              GO TO 90
                          END IF
                      END IF
    3             CONTINUE
              END IF
          END IF
          CALL UNPERT(IPAIR)
          GO TO 100
      END IF
*
*       Perform KS prediction of U, UDOT & H.
      CALL KSPRED(IPAIR,I1,I,BODYIN,UI,UIDOT,XI,VI)
*
*       Obtain the perturbing force & derivative.
      CALL KSPERT(I1,NNB0,XI,VI,FP,FD)
*
*       Save old radial velocity & relative perturbation and set new GAMMA.
      RDOT = TDOT2(IPAIR)
*     G0 = GAMMA(IPAIR)
      GI = SQRT(FP(1)**2 + FP(2)**2 + FP(3)**2)*R(IPAIR)**2*BODYIN
      GAMMA(IPAIR) = GI
*
*       Apply the Hermite corrector.
      CALL KSCORR(IPAIR,UI,UIDOT,FP,FD,TD2,TDOT4,TDOT5,TDOT6)
*
*       Increase regularization time-step counter and update the time.
      NSTEPU = NSTEPU + 1
      T0(I1) = TIME
*       Check for early return during termination (called from KSTERM).
      IF (IPHASE.NE.0) GO TO 100
*
*       Define useful scalars.
      RI = R(IPAIR)
      HI = H(IPAIR)
*
*       Initialize termination indicator and check for large perturbation.
      IQ = .FALSE.
      IF (GI.GT.0.25) GO TO 2
      IF (GI.LT.0.1) THEN
          JCOMP = 0
          GO TO 20
      END IF
      CALL FLYBY(I,ITERM)
      IF (ITERM.EQ.0.AND.KSTAR(I).LT.0) THEN
          IQ = .TRUE.
      END IF
      IF (ITERM.EQ.1) THEN
*       Delay chain regularization search until near end of block-step.
          IF (TIME + STEP(I1).GT.TBLOCK) THEN
              CALL IMPACT(I)
              IF (IPHASE.GT.0) GO TO 100
          END IF
      ELSE IF (ITERM.EQ.2) THEN
          IQ = .TRUE.
          GO TO 20
      ELSE
          GO TO 20
      END IF
*
*       Find the dominant body for large perturbations.
    2 S = 4.0*STEP(I)
      FMAX = BODY(I)/RI**2
*       Initialize JCOMP for prediction and optional diagnostics in KSTERM.
      JCOMP = 0
      DO 10 L = 2,NNB0+1
          J = LIST(L,I1)
*       Only search bodies within twice the c.m. time-step.
          IF (STEP(J).GT.S) GO TO 10
*       Compare strong perturber and either component with current pair.
          DO 5 K = I1,I2
              RIJ2 = (X(1,J) - X(1,K))**2 + (X(2,J) - X(2,K))**2 +
     &                                      (X(3,J) - X(3,K))**2
              IF (BODY(J) + BODY(K).GT.RIJ2*FMAX) JCOMP = J
    5     CONTINUE
   10 CONTINUE
*
*       Set termination if strong perturber <= N forms dominant pair.
      IF (JCOMP.GT.0.OR.GI.GT.1.0) THEN
*       Check optional binary search.
*         IF (KZ(4).GT.0) THEN
*             DGAM = GI - G0
*             K = KZ(4)
*             CALL EVOLVE(IPAIR,K)
*         END IF
          IF (JCOMP.LE.N.OR.GI.GT.2.0) IQ = .TRUE.
      END IF
*
*       Check termination of hyperbolic encounter (R > R0 or R > 2*RMIN).
   20 IF (HI.GT.0.0D0.AND.NAME(I).GT.0) THEN
          IF ((RI.GT.R0(IPAIR).AND.GI.GT.GMAX).OR.RI.GT.2.0*RMIN.OR.
     &        (GI.GT.0.5.AND.TD2.GT.0.0)) THEN
*       Skip termination delay in case of velocity kick (cf. routine KSTERM).
              IF (HI.LT.100.0.OR.GI.GT.0.1.OR.RI.GT.5.0*RMIN) THEN
                  IQ = .TRUE.
              END IF
          END IF
      END IF
*
*       Choose basic regularized step using binding energy or lower limit.
      IF (ABS(HI).GT.ECLOSE) THEN
          W1 = 0.5/ABS(HI)
      ELSE
          W1 = R0(IPAIR)*BODYIN
          W2 = 0.5/ABS(HI)
          W1 = MIN(W1,W2)
          IF (RI.GT.R0(IPAIR)) W1 = W1*R0(IPAIR)/RI
*       Maximum square step for soft binaries & weak hyperbolic pairs.
          IF (HI.LT.0.0D0) THEN
*       Include case of hard binary with massive components or merger.
              W2 = -0.5/HI
              IF (NAME(I).LT.0) THEN
                  W1 = RI*BODYIN
              END IF
              W1 = MIN(W2,W1)
          END IF
      END IF
*
*       Include perturbation factor in predicted step.
      IF (GI.LT.0.0005) THEN
*       Use second-order expansion of cube root for small perturbations.
          W3 = 333.3*GI
          W2 = SQRT(W1)/(1.0 + W3*(1.0 - W3))
      ELSE
          W3 = 1.0 + 1000.0*GI
          W2 = SQRT(W1)/W3**0.3333
      END IF
*
*       Form new regularized step (include inertial factor).
      DTU = ETAU*W2
      DTU = MIN(1.2*DTAU(IPAIR),DTU)
*
*       Include convergence criterion DH = H'*DTU + H''*DTU**2/2 = 0.001*|H|.
      IF (GI.GT.1.0D-04) THEN
          DH = 1.0E-03*MAX(ABS(H(IPAIR)),ECLOSE)
          XF = 2.0*DH/ABS(HDOT2(IPAIR))
          YF = HDOT(IPAIR)/HDOT2(IPAIR)
          DTU1 = SQRT(XF + YF**2) - ABS(YF)
          DTU = MIN(DTU1,DTU)
      END IF
*
*       Check pericentre step reduction for perturbed spiral.
      IF (KSTAR(I).EQ.-2.AND.TD2.LT.0.0) THEN
          SEMI = -0.5*BODY(I)/HI
          IF (RI.LT.SEMI) THEN
*       Predict radial velocity for step DTU (note scaled coefficients).
              RD = ((ONE3*TDOT5*DTU + TDOT4)*DTU +
     &                                TDOT3(IPAIR))*DTU + 2.0*TD2
*       Adopt pericentre step with 1% safety factor (small TDOT4 is OK).
              IF (RD.GT.0.0) THEN
                  A2 = 0.5*TDOT3(IPAIR)/TDOT4
                  DTU1 = SQRT(A2**2 - 2.0*TD2/TDOT4) - A2
                  DTU1 = 1.01*MAX(DTU1,1.0D-10)
                  DTU = MIN(DTU1,DTU)
              END IF
          END IF
      END IF
*
*       Reset reference energy and generate new Stumpff coefficients.
      H0(IPAIR) = H(IPAIR)
   30 Z = -0.5D0*H(IPAIR)*DTU**2
      CALL STUMPF(IPAIR,Z)
      Z5 = SF(6,IPAIR)
      Z6 = SF(7,IPAIR)
      DT12 = ONE12*DTU*Z6
*
*       Convert to physical time units modified by Stumpff coefficients.
      STEP(I1) = (((((TDOT6*DT12 + TDOT5*Z5)*0.2*DTU + 0.5D0*TDOT4)*DTU
     &                     + TDOT3(IPAIR))*ONE6*DTU + TD2)*DTU + RI)*DTU
*
*       Ensure that regularized step is smaller than the c.m. step.
      IF (STEP(I1).GT.STEP(I).AND.HI.LT.0) THEN
          DTU = 0.5D0*DTU
          GO TO 30
      END IF
      DTAU(IPAIR) = DTU
*
*       See whether the KS slow-down procedure is activated.
      IMOD = KSLOW(IPAIR)
      IF (IMOD.GT.1) THEN
          ZMOD = FLOAT(ISLOW(IMOD))
          STEP(I1) = ZMOD*STEP(I1)
      END IF
*
*       Check diagnostic print option.
      IF (KZ(10).GE.3.and.rank.eq.0) THEN
          WRITE (6,40)  IPAIR, TIME+TOFF, H(IPAIR), RI, DTAU(IPAIR),
     &                  GI, STEP(I1), LIST(1,I1), IMOD
   40     FORMAT (3X,'KS MOTION',I6,2F10.4,1P,4E10.2,2I4)
      END IF
*
*       Employ special termination criterion in merger case.
      IF (NAME(I).LT.0) THEN
*       Terminate if apocentre perturbation > 0.25 (R > SEMI) or GI > 0.25.
          IF (HI.LT.0.0) THEN
              SEMI = -0.5*BODY(I)/HI
              ECC2 = (1.0 - RI/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
              A0 = SEMI*(1.0 + SQRT(ECC2))/RI
              GA = GI*A0*A0*A0
              IF (GA.GT.0.25.AND.RI.GT.SEMI) IQ = .TRUE.
              IF (GI.GT.0.1.AND.RI.GT.RMIN) IQ = .TRUE.
              IF (GI.GT.0.25) IQ = .TRUE.
              IF (MIN(BODY(I1),BODY(I2)).LT.0.05*BODYM) THEN
                  IF (GI.GT.2.0D-04) IQ = .TRUE.
              END IF
          ELSE
              IF (TD2.GT.0.0.AND.(GI.GT.GMAX.OR.RI.GT.RMIN)) IQ = .TRUE.
          END IF
          IF (.NOT.IQ) GO TO 60
      END IF
*
*       Delay termination until end of block for large perturbations.
      IF (IQ) THEN
          DTR = TBLOCK - TIME
*         if(rank.eq.0)
*    &    WRITE (6,45)  IPAIR, TTOT, GI, RI, DTR, STEP(I1)
*  45     FORMAT (' TERM TEST    KS T G R DTR DT  ',
*    &                           I4,F10.4,F7.3,1P,E10.2,2E9.1)
          IF (DTR.LT.STEP(I1)) GO TO 90
      END IF
*
*       Check standard termination criterion (suppress on IQ = .true.).
*     IF (RI.GT.R0(IPAIR)) THEN
      IF (RI.GT.R0(IPAIR).AND.RI.GT.2.0*RMIN.AND..NOT.IQ) THEN
*
*       See whether termination can be delayed for intermediate energies.
          IF (RI.GT.RMIN) THEN
              A3 = RMIN*ABS(HI)*BODYIN
              GIMAX = A3*A3*A3
              IF (GI.LT.GIMAX) GO TO 60
*       Check updating of R0 for newly hardened binary orbit.
              IF (HI.LT.-ECLOSE) THEN
                  SEMI = -0.5*BODY(I)/HI
                  R0(IPAIR) = MAX(RMIN,2.0D0*SEMI)
                  GO TO 70
              END IF
          END IF
*
*         IF (KZ(4).GT.0.AND.GI.GT.GPRINT(1)) THEN
*             DGAM = GI - G0
*             K = KZ(4)
*             DO 55 L = 2,K
*                 IF (GI.LT.GPRINT(L)) THEN
*                     CALL EVOLVE(IPAIR,L-1)
*                     GO TO 90
*                 END IF
*  55         CONTINUE
*             CALL EVOLVE(IPAIR,K)
*         END IF
          GO TO 90
      END IF
*
*       End integration cycle for hyperbolic motion.
   60 IF (HI.GE.0.0D0) THEN
          IF (RDOT*TD2.LT.0.0D0) THEN
*       Determine pericentre for hyperbolic two-body motion.
              SEMI = -0.5D0*BODY(I)/HI
              ECC2 = (1.0 - RI/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
              QPERI = SEMI*(1.0D0 - SQRT(ECC2))
              DMIN2 = MIN(DMIN2,QPERI)
*
*       Check optional tidal interaction or stellar collision.
              IF (KZ(19).GE.3.AND.KSTAR(I).LT.10) THEN
                  VINF = SQRT(2.0*HI)*VSTAR
                  KS1 = KSTAR(I1)
                  KS2 = KSTAR(I2)
                  RX = MAX(RADIUS(I1),RADIUS(I2))
*       Determine maximum periastron factor for capture (VINF in km/sec).
                  IF (KZ(27).LE.2) THEN
                      RFAC = RPMAX2(RADIUS(I1),RADIUS(I2),BODY(I1),
     &                              BODY(I2),KS1,KS2,VINF)
                      RCAP = RFAC*RX
                  ELSE
                      DV = SQRT(2.0*HI)
*       Note that Quinlan & Shapiro function returns actual distance.
                      RCAP = RPMAX(BODY(I1),BODY(I2),VSTAR,DV,QPERI)
                  END IF
                  IF (QPERI.LT.5.0*RX) THEN
                      if(rank.eq.0)
     &                WRITE (54,54)  TTOT, NAME(I1), NAME(I2), KS1,
     &                               KS2, VINF, RCAP*SU, RX*SU, QPERI*SU
   54                 FORMAT (' CLOSE   T NAM K* VINF RCAP RX QP  ',
     &                                  F7.1,2I6,2I4,F6.2,3F6.1)
                  END IF
                  IF (QPERI.LT.RCAP) THEN
                      J1 = I1
                      IF (RADIUS(I2).GT.RADIUS(I1)) J1 = I2
                      FAC = 0.5*BODY(I)/BODY(J1)
*       Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
                      IF (KZ(27).LE.2) THEN
                          RCOLL = 1.7*FAC**0.3333*RADIUS(J1)
                      ELSE
                          RCOLL = 6.0*BODY(I)/CLIGHT**2
                      END IF
                      if(rank.eq.0)then
                      WRITE (55,58)  TTOT, IPAIR, NAME(I1), NAME(I2),
     &                               KS1, KS2, KSTAR(I), VINF
   58                 FORMAT (' RPMAX:    T KS NAM K* VINF ',
     &                                    F7.1,I5,2I6,3I4,F7.2)
                      WRITE (55,59)  SQRT(ECC2), HI, R(IPAIR), SEMI,
     &                               QPERI, BODY(I1), BODY(I2),
     &                               BODY(I)*ZMBAR
   59                 FORMAT (' E H R A QP BODY MT ',
     &                          F9.5,1P,6E10.2,0P,F6.1)
                      end if
                      RI2 = 0.0
                      VI2 = 0.0
                      DO 61 K = 1,3
                          RI2 = RI2 + (X(K,I) - CMR(K))**2
                          VI2 = VI2 + XDOT(K,I)**2
   61                 CONTINUE
                      if(rank.eq.0)
     &                WRITE (55,62)  SQRT(RI2)/RC, SQRT(VI2)*VSTAR,
     &                               RHOD, RADIUS(I1)*SU, RADIUS(I2)*SU,
     &                               RCAP, RADIUS(J1)/QPERI, RCOLL/QPERI
   62                 FORMAT (' r/RC V* <C> R* RCAP R1/QP RCOLL/QP ',
     &                          2F5.1,3F6.1,3F5.1)
                      CALL FLUSH(55)
                      IF (QPERI.LT.RCOLL) THEN
*       Obtain KS variables at pericentre before merging into one body.
                          CALL KSPERI(IPAIR)
                          KSPAIR = IPAIR
                          IQCOLL = -2
                          CALL CMBODY(QPERI,2)
                      ELSE IF (KSTAR(I).GE.0.AND.KZ(27).GT.0) THEN
                          CALL KSTIDE(IPAIR,QPERI)
                      END IF
                  END IF
*       Check options for artificial collisions.
              ELSE IF (KZ(27).EQ.-1.AND.KZ(13).LT.0) THEN
                  RFAC = 2.0
                  IF (QPERI.LT.RFAC*MAX(RADIUS(I1),RADIUS(I2))) THEN
                      J1 = I1
                      IF (RADIUS(I2).GT.RADIUS(I1)) J1 = I2
                      FAC = 0.5*BODY(I)/BODY(J1)
*       Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
                      RCOLL = 1.7*FAC**0.3333*RADIUS(J1)
                      IF (QPERI.LT.RCOLL) THEN
                          CALL TOUCH(IPAIR,I1,I2,RCOLL)
                      END IF
                  END IF
              END IF
          END IF
          GO TO 100
      END IF
*
*       Check perturbation threshold (H < 0 & GAMMA > GMAX).
*     IF (KZ(4).EQ.0.OR.G0.LT.GMAX) GO TO 70
*
*     K = KZ(4)
*     DO 65 L = 1,K
*         IF ((G0 - GPRINT(L))*(GI - GPRINT(L)).LE.0.0) THEN
*             IF (L.EQ.1) THEN
*                 W1 = -0.5*BODY(I)/HI
*                 W2 = W1*BODYIN
*                 TK = TWOPI*W1*SQRT(W2)
*             END IF
*
*       Estimate smallest permitted output interval at new level.
*             DTCRIT = TK*ORBITS(L)
*             IF (TIME - TLASTB(L).LT.DTCRIT) GO TO 70
*             DGAM = GI - G0
*             CALL EVOLVE(IPAIR,L)
*             GO TO 70
*         END IF
*  65 CONTINUE
*
*       Check for partial reflection during approach (NB! only IMOD = 1).
*  70 IF (GI.LT.GMIN.AND.TD2.LT.0.0D0) THEN
*       Skip apocentre position itself.
*         IF (RDOT.LT.0.0D0.AND.IMOD.EQ.1) THEN
*             IF (KZ(25).GT.0) CALL FREEZE(IPAIR)
*             GO TO 100
*         END IF
*     END IF
*
*       Determine new perturbers for binary at apocentre turning point.
   70 IF (RDOT*TD2.GE.0.0D0) GO TO 100
      SEMI = -0.5D0*BODY(I)/HI
*
*       Check minimum two-body separation just after pericentre.
      IF (RDOT.LT.0.0D0) THEN
*       Obtain pericentre by Mikkola's algorithm (GAMMA < 0.001).
          IF (GI.LT.0.001) THEN
              CALL PERI(UI,UIDOT,RI,BODY(I1),BODY(I2),QPERI)
          ELSE
              QPERI = RI
          END IF
          DMIN2 = MIN(DMIN2,QPERI)
*
*       Check optional tidal interaction or stellar collision (skip merger).
          IF (KZ(19).GE.3.AND.KSTAR(I).LE.10.AND.NAME(I).GT.0) THEN
              RFAC = 5.0
              IF (KZ(27).LE.2) THEN
                  IF (KZ(27).EQ.1) RFAC = 4.0
                  RX = RFAC*MAX(RADIUS(I1),RADIUS(I2))
              ELSE
                  RX = RPMIN(BODY(I1),BODY(I2),VSTAR,HI,QPERI)
              END IF
              IF (QPERI.LT.RX) THEN
                  J1 = I1
                  IF (RADIUS(I2).GT.RADIUS(I1)) J1 = I2
                  FAC = 0.5*BODY(I)/BODY(J1)
*       Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
                  IF (KZ(27).LE.2) THEN
                      RCOLL = 1.7*FAC**0.3333*RADIUS(J1)
                  ELSE
                      RCOLL = 6.0*BODY(I)/CLIGHT**2
                  END IF
                  IF (QPERI.LT.RCOLL) THEN
*       Obtain KS variables at pericentre before merging into one body.
                      CALL KSPERI(IPAIR)
                      KSPAIR = IPAIR
                      IQCOLL = -2
                      CALL CMBODY(QPERI,2)
                  ELSE IF (KSTAR(I).GE.0) THEN
*       Distinguish between sequential, standard and GR circularization.
                      IF (KZ(27).EQ.1) THEN
                          ICIRC = 1
                          TC = 0.0
                      ELSE IF (KZ(27).EQ.2) THEN
                          ECC2 = (1.0 - RI/SEMI)**2 +
     &                                    TDOT2(IPAIR)**2/(BODY(I)*SEMI)
                          ECC = SQRT(ECC2)
                          ICIRC = 0
                          CALL TCIRC(QPERI,ECC,I1,I2,ICIRC,TC)
                      ELSE
                          ICIRC = 1
                          TC = 0.0
                      END IF
                      IF (KSTAR(I).GE.10) ICIRC = 0
*       Skip tidal effects for circularization time above 100 Myr (07/08).
                      IF (ICIRC.GT.0.AND.KZ(27).GT.0.AND.
     &                    TC.LT.100.0) THEN
                          CALL KSTIDE(IPAIR,QPERI)
                      END IF
                  END IF
              END IF
*       Check for perturbed spiral or chaos case (skip collision).
              IF (KSTAR(I).EQ.-2.AND.IPHASE.EQ.0) THEN
                  CALL SPIRAL(IPAIR)
              ELSE IF (KSTAR(I).EQ.-1.AND.IPHASE.EQ.0) THEN
                  CALL KSTIDE(IPAIR,QPERI)
              END IF
*       Check options for artificial collisions.
          ELSE IF (KZ(27).EQ.-1.AND.KZ(13).LT.0) THEN
              RFAC = 2.0
              IF (QPERI.LT.RFAC*MAX(RADIUS(I1),RADIUS(I2))) THEN
                  J1 = I1
                  IF (RADIUS(I2).GT.RADIUS(I1)) J1 = I2
                  FAC = 0.5*BODY(I)/BODY(J1)
*       Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
                  RCOLL = 1.7*FAC**0.3333*RADIUS(J1)
                  IF (QPERI.LT.RCOLL) THEN
                      CALL TOUCH(IPAIR,I1,I2,RCOLL)
                  END IF
              END IF
          END IF
          GO TO 100
      END IF
*
*       Save maximum separation of persistent binary.
      RMAX = MAX(RMAX,RI)
*
*       Check binary reference radius or merger stability criterion.
      IF (NAME(I).GT.0) THEN
*       Update termination length scale in case of initial soft binary.
          EB = BODY(I1)*BODY(I2)*HI*BODYIN
          IF (EB.LT.EBH) R0(IPAIR) = MAX(RMIN,2.0*SEMI)
      ELSE
          ECC2 = (1.0 - RI/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
          ECC = SQRT(ECC2)
          RP = SEMI*(1.0 - ECC)*(1.0 - 2.0*GI)
*       Find merger index.
          IM = 0
          DO 72 K = 1,NMERGE
              IF (NAMEM(K).EQ.NAME(I)) IM = K
   72     CONTINUE
*       Exclude inner planets from the general stability test.
          IF (MIN(BODY(I1),BODY(I2)).LT.0.05*BODYM) THEN
              IF (RP.LT.R0(IPAIR)) GO TO 90
          END IF
*       Assess the stability inside critical pericentre (safety factor 1.04).
          IF (RP.LT.1.04*R0(IPAIR)) THEN
              EOUT = ECC
*       Increase tolerance near sensitive stability boundary (RM 10/2008).
              IF (EOUT.GT.0.90) THEN
                  DE = 0.5*(1.0 - EOUT)
                  DE = MIN(DE,0.01D0)
*       Add extra amount 0.011 to avoid switching.
                  IF (ECC.GT.0.9) DE = DE + 0.011
                  EOUT = EOUT - DE
              END IF
*       Note: assessment needs to use same eccentricity as for acceptance.
              CALL ASSESS(IPAIR,IM,EOUT,SEMI,ITERM)
              IF (ITERM.GT.0) THEN
                  INSTAB = INSTAB + 1
                  GO TO 90
              END IF
          END IF
*       Check possible eccentricity modulation or t_circ update.
          IF (IM.GT.0.AND.(TIME.GT.TMDIS(IM).OR.
     &        TMDIS(IM).GT.1.0D+06)) THEN
              IF (KZ(27).EQ.2) THEN
                  CALL ECCMOD(I,ITERM)
                  IF (ITERM.GT.0) THEN
*                     if(rank.eq.0)
*    &                WRITE (6,76)  RP, R0(IPAIR)
*  76                 FORMAT (' ECCMOD TERM    RP R0 ',1P,2E10.2)
                      GO TO 90
                  END IF
*       Consider both inner and possible outer circularizing binary.
                  DO 78 K = 1,NCHAOS
                      IF (NAMEC(K).EQ.NZERO - NAMEM(IM).AND.
     &                    KSTARM(IM).EQ.-2) THEN
*       Update unperturbed binary if T - TOSC > 10 Myr (cf. IMPACT & DECIDE).
                          IF ((TIME - TOSC(K))*TSTAR.GT.10.0) GO TO 90
                      END IF
                      IF (NAMEC(K).EQ.NAMEG(IM).AND.
     &                    KSTARM(IM).EQ.-2) THEN
                          IF ((TIME - TOSC(K))*TSTAR.GT.10.0) GO TO 90
                      END IF
*       Note: perturbed binary is treated if pericentre before next IMPACT.
   78             CONTINUE
              END IF
          END IF
      END IF
*
*       Produce diagnostics for any circularizing perturbed binary.
      IF (KSTAR(I).EQ.-2.AND.GI.GT.0.01) THEN
          ECC = RI/SEMI - 1.0
          QPS = SEMI*(1.0 - ECC)/MAX(RADIUS(I1),RADIUS(I2))
          ZM = BODY(I1)/BODY(I2)
          if(rank.eq.0)
     &    WRITE (21,80)  TTOT, NAME(I1), NAME(I2), KSTAR(I1), KSTAR(I2),
     &                   LIST(1,I1), QPS, ZM, GI, ECC, SEMI
   80     FORMAT (' PERT SPIRAL    T NAM K* NP QP/S M1/M2 G E A ',
     &                             F11.4,2I6,3I3,2F5.1,2F8.4,1P,E10.2)
          CALL FLUSH(21)
      END IF
*
*       See whether KS slow-down procedure should be (re)-checked (no Chaos).
      IF (KZ(26).GT.0.AND.KSTAR(I).GE.0) THEN
          KMOD = RANGE*GMIN/MAX(GI,1.0D-10)
          IF (KMOD.GT.1.OR.IMOD.GT.1) THEN
              CALL KSMOD(IPAIR,KMOD)
              IF (KMOD.LT.0) GO TO 100
          END IF
      END IF
*
*       Set approximate value of next period with perturbation included.
      TK = TWOPI*SEMI*SQRT(SEMI*BODYIN)*(1.0 + GI)
      IF (IMOD.GT.1) THEN
          TK = ZMOD*TK
      END IF
*
*       Use old perturber list if next apocentre is before the c.m. step.
      IF (TIME + TK.LT.T0(I) + STEP(I)) THEN
          GO TO 100
      END IF
*
*       Select new perturbers (J = N adopted for unperturbed Chaos).
      CALL KSLIST(IPAIR)
*
*       Check rectification of chaotic spiral at start of unperturbed motion.
      IF (KSTAR(I).EQ.-2.AND.LIST(1,I1).EQ.0) THEN
          DMR = 0.D0
          CALL CHRECT(IPAIR,DMR)
          IF (IPHASE.LT.0) GO TO 100
      ELSE
          CALL KSRECT(IPAIR)
      END IF
*
*       Check optional search criterion for multiple encounter or merger.
      IF (KZ(15).GT.0.AND.STEP(I).LT.DTMIN) THEN
          CALL IMPACT(I)
      END IF
      GO TO 100
*
*       Terminate regularization of current pair (IPAIR set in KSPAIR).
   90 KSPAIR = IPAIR
*       Set indicator for calling KSTERM in MAIN (permits phase overlay).
      IPHASE = 2
*       Check case of hierarchical binary.
      IF (NAME(I).LT.0) IPHASE = 7
*
  100 RETURN
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
C  =====================================================================
C  pgm: SHSAVP .. Buffer and save/restore parameter code
C
C  use:      CALL SHSAVP(CMD,ITEM,PARCOD)
C
C   in: CMD ....... routine command to save or retrieve p-code - CHAR*1
C   in:               'I' ... Initialize buffer parameter code to blank
C   in:               'P' ... Put parameter code into the buffer
C   in:               'G' ... Get parameter code from the buffer
C   in:               'S' ... Save buffer parameter code in array
C   in:               'R' ... Restore buffer parameter code from array
C   in: ITEM ...... field number of p-code to save or retrieve - INT
C  i/o: PARCOD .... parameter code - CHAR*8
C
C  cmt: Hardwired for 200 fields (or items).
C  =====================================================================
      SUBROUTINE SHSAVP(CMD,ITEM,PARCOD)

      CHARACTER*1   CMD
      INTEGER       ITEM,II
      CHARACTER*8   PARCOD,SVPCOD,SVARRP(200)

      SAVE      SVPCOD,SVARRP
C
C    ================================= RCS keyword statements ==========
      CHARACTER*68     RCSKW1,RCSKW2
      DATA             RCSKW1,RCSKW2 /                                 '
     .$Source: /fs/hseb/ob5/rfc/ofs/src/shefpars_driv/RCS/shsavp.f,v $
     . $',                                                             '
     .$Id: shsavp.f,v 1.2 1997/12/31 20:51:01 page Exp $
     . $' /
C    ===================================================================
C

      DATA      SVPCOD / ' ' /

        IF (CMD .EQ. 'I') THEN
            SVPCOD = ' '

        ELSEIF (CMD .EQ. 'P') THEN
            SVPCOD = PARCOD

        ELSEIF (CMD .EQ. 'G') THEN
            PARCOD = SVPCOD

        ELSEIF (CMD .EQ. 'S') THEN
             II = ITEM
            SVARRP(II) = SVPCOD

        ELSEIF (CMD .EQ. 'R') THEN
             II = ITEM
            SVPCOD = SVARRP(II)

        ENDIF

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
The Address(Eichhorn Family House, 38.540383717063, 121.744092458215) is located at the east end of campus, tucked behind the Center for Child & Family Studies. It is part of the Human Development Department of Human and Community Development. It was funded by private donation.

It houses several offices used by graduate students and faculty and well as research space. Some of these rooms are designed with oneway mirrors and inconspicuous cameras to record and monitor research participants.

The facility was designed to look and feel like a home rather than a more typical research space. It includes a fully functional kitchen and dining area, a living room with large windows, and previously had working laundry facilities. It was designed this way in order to make participants, in particular families with small children, feel more comfortable than they might in a more traditional research facility.

The building is not publicly accessible, and requires keycards to access it; these are only available to people who are involved in research there.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*----------------------------------------------------------------------*
      subroutine r12_form_post(form,nterms,
     &     idxtbar,idxcbar,idxham,idxtop,idxc12, iprint,
     &     op_info)
*----------------------------------------------------------------------*
*     resort terms according to tbar blocks and count terms
*     print some statistics
*----------------------------------------------------------------------*

      implicit none

      include 'stdunit.h'
      include 'opdim.h'
      include 'mdef_operator_info.h'
      include 'def_contraction.h'
      include 'def_formula_item.h'
      include 'def_formula_item_list.h'

      type(formula_item), intent(inout), target ::
     &     form
      integer, intent(in) ::
     &     idxtbar, idxcbar, idxham, idxtop, idxc12, iprint
      type(operator_info), intent(in) ::
     &     op_info
      integer, intent(out) ::
     &     nterms

      integer ::
     &     n_commu(1:15)
      type(formula_item_list), target ::
     &     fpl_reo
      type(formula_item_list), pointer ::
     &     fpl_pnt
      type(formula_item), pointer ::
     &     form_pnt
      type(operator), pointer ::
     &     op_tbar, op_cbar
      integer ::
     &     nblk_tbar, nblk_cbar, iblk

      nterms = 0
      ! initialize reodering list
      call init_formula_plist(fpl_reo)
      fpl_pnt => fpl_reo
      ! point to "initialize" item
      if (form%command.ne.command_set_target_init)
     &     call quit(1,'r12_form_post','wrong assumptions? (1)')
      fpl_pnt%item => form
      ! and advance to first "add"
      form_pnt => form%next
      if (form_pnt%command.ne.command_add_contribution)
     &     call quit(1,'r12_form_post','wrong assumptions? (2)')

      op_tbar => op_info%op_arr(abs(idxtbar))%op
      nblk_tbar = op_tbar%n_occ_cls
      op_cbar => op_info%op_arr(abs(idxcbar))%op
      nblk_cbar = op_cbar%n_occ_cls

      if (iprint.gt.0) then
        write(lulog,'(x,78("-"))')
        write(lulog,'(x,a)')
     &  ' L    number  n-fold commutators'
        write(lulog,'(x,a)')
     &  'class  of      0   1       2           3                4'
        write(lulog,'(x,a)')
     &  '      terms        T   C  TT  TC  CC TTT TTC TCC CCC TTTT '//
     &       'TTTC TTCC TCCC CCCC'
        write(lulog,'(x,78("-"))')
      end if

      ! loop over Tbar blocks (0 is the unit-projection)
      do iblk = 0, nblk_tbar
        if (iblk.eq.0) then
          call collect_terms_w_op(fpl_pnt,form_pnt,
     &         2,(/idxtbar,idxcbar/),-1,0)
          if (associated(fpl_pnt%next)) fpl_pnt => fpl_pnt%next
        else
          call collect_terms_w_op(fpl_pnt,form_pnt,1,idxtbar,iblk,1)
          if (associated(fpl_pnt%next)) fpl_pnt => fpl_pnt%next
        end if

        call r12_count_terms(fpl_pnt,idxtop,idxc12,n_commu,op_info)

        do while(associated(fpl_pnt).and.associated(fpl_pnt%next))
          fpl_pnt => fpl_pnt%next
        end do

        if (iprint.ge.10) then
c          write(lulog,'("??",x,"L",i2,3x,i6,x,15(i4))')
          write(lulog,'("??",x,"L",i1,x,i6,x,10(i4),5(i5))')
     &       iblk,sum(n_commu(1:15)),n_commu(1:15)
        else if (iprint.gt.0) then
          write(lulog,'(3x,"L",i1,x,i6,x,10(i4),5(i5))')
     &       iblk,sum(n_commu(1:15)),n_commu(1:15)
        end if
        nterms = nterms + sum(n_commu(1:15))

      end do
      ! loop over Cbar blocks
      do iblk = 1, nblk_cbar
        call collect_terms_w_op(fpl_pnt,form_pnt,1,idxcbar,iblk,1)
        if (associated(fpl_pnt%next)) then
          fpl_pnt => fpl_pnt%next
          call r12_count_terms(fpl_pnt,idxtop,idxc12,n_commu,op_info)
        else
          n_commu = 0
        end if

        do while(associated(fpl_pnt%next))
          fpl_pnt => fpl_pnt%next
        end do

        if (iprint.ge.10) then
          write(lulog,'("??",x,"C",i1,x,i6,x,10(i4),5(i5))')
c          write(lulog,'("??",x,"C",i2,3x,i6,x,15(i4))')
     &       iblk+1,sum(n_commu(1:15)),n_commu(1:15)
        else if (iprint.gt.0) then
          write(lulog,'(3x,"C",i1,x,i6,x,10(i4),5(i5))')
     &       iblk+1,sum(n_commu(1:15)),n_commu(1:15)
        end if
        nterms = nterms + sum(n_commu(1:15))

      end do

      if (iprint.gt.0)
     &     write(lulog,'(x,78("-"))')

      call relink_formula_list(form,fpl_reo)

      call dealloc_formula_plist(fpl_reo)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
module model_fluid_primary_fields_regularized
#include <messenger.h>
    use mod_kinds,          only: rk
    use mod_constants,      only: HALF, ONE, ZERO
    use mod_fluid,          only: Rgas, gam
    use type_model,         only: model_t
    use type_chidg_worker,  only: chidg_worker_t
    use DNAD_D
    implicit none

    !> This model computes regularized values for fluid primary field variables by
    !! enforcing positivity of density and energy.
    !!
    !! @author  Eric M. Wolf
    !! @date    08/05/2019
    !!
    !--------------------------------------------------------------------------------
    type, extends(model_t)  :: fluid_primary_fields_regularized_t

        real(rk)    :: window_size_density      = 1.0e-2_rk
        real(rk)    :: minimum_value_density    = 1.0e-11_rk
        logical     :: regularize_density       = .true.

        real(rk)    :: window_size_energy       = 1.0e-2_rk
        real(rk)    :: minimum_value_energy     = 1.0e-11_rk
        logical     :: regularize_energy        = .true.

    contains

        procedure   :: init
        procedure   :: compute

    end type fluid_primary_fields_regularized_t
    !***************************************************************************************

contains

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    08/05/2019
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(fluid_primary_fields_regularized_t), intent(inout)   :: self

        integer             :: unit, msg
        logical             :: file_exists
        real(rk)    :: window_size_density, minimum_value_density, window_size_energy, minimum_value_energy
        logical     :: regularize_density, regularize_energy

        namelist /regularize_fluid_primary_fields/ window_size_density, minimum_value_density, regularize_density, &
            window_size_energy, minimum_value_energy, regularize_energy

        call self%set_name('Regularized Fluid Primary Fields')
        call self%set_dependency('f(Q-)')

        call self%add_model_field('Density')
        call self%add_model_field('Energy')

        inquire(file='regularize_fluid_primary_fields.nml', exist=file_exists)
        if (file_exists) then
            open(newunit=unit,form='formatted',file='regularize_fluid_primary_fields.nml')
            read(unit,nml=regularize_fluid_primary_fields,iostat=msg)
            if (msg == 0) self%window_size_density = window_size_density
            if (msg == 0) self%minimum_value_density = minimum_value_density
            if (msg == 0) self%regularize_density = regularize_density
            if (msg == 0) self%window_size_energy = window_size_energy
            if (msg == 0) self%minimum_value_energy = minimum_value_energy
            if (msg == 0) self%regularize_energy = regularize_energy
            close(unit)
        end if

    end subroutine init
    !***************************************************************************************

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    08/05/2019
    !!
    !--------------------------------------------------------------------------------
    subroutine compute(self,worker)
        class(fluid_primary_fields_regularized_t),     intent(in)      :: self
        type(chidg_worker_t),   intent(inout)   :: worker

        type(AD_D), dimension(:),   allocatable :: &
            density, energy

        !
        ! Interpolate solution to quadrature nodes
        !
        density = worker%get_field('Density',    'value')
        energy  = worker%get_field('Energy',     'value')

        if (self%regularize_density) density = density*sin_ramp(density, ZERO, self%window_size_density) + self%minimum_value_density
        if (self%regularize_energy) energy = energy*sin_ramp(energy, ZERO, self%window_size_energy) + self%minimum_value_energy

        call worker%store_model_field('Density',    'value', density)
        call worker%store_model_field('Energy', 'value', energy)

    end subroutine compute
    !***************************************************************************************

end module model_fluid_primary_fields_regularized

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG