subroutine physdem1pal(filename,lonfi,latfi,nsoil,nq,
     .          phystep,day_ini,
     .          time,tsurf,tsoil,emis,q2,qsurf,
     .          airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe,
     .          oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal)



      use planet_h         
      use comgeomfi_h 
      use comsoil_h
      implicit none
c-------------------------------------------------------------
c
c create physics (re-)start data file "restartfi.nc"
c
c
c
#include "dimensions.h"
#include "paramet.h"
c-----------------------------------------------------------------------
#include "comvert.h"
#include "comgeom2.h"
#include "control.h"
#include "comdissnew.h"
#include "logic.h"
#include "ener.h"
#include "netcdf.inc"
#include "dimphys.h"
#include"advtrac.h"
#include"callkeys.h"

      INTEGER nid,iq
      INTEGER, parameter :: ivap=1
      REAL, parameter :: qsolmax= 150.0
      character (len=*) :: filename
      character (len=7) :: str7

      REAL day_ini
      REAL oblipal,peri_daypal,eccpal   !TB16 : change of obliquity/periday/eccentricity for paleo run
      REAL peripal,aphepal   !TB16 : change of perihelion/aphelion from eccentricity
      REAL tpalnew   !TB16 : change of obliquity for paleo run
      REAL adjustnew   !TB21 : change in N2 ice albedo
      REAL phisfipal(ngridmx)   !TB16 : change of phisfi for paleo run
      INTEGER nsoil,nq
      integer ierr,idim1,idim2,idim3,idim4,idim5,nvarid

c
      REAL phystep,time
      REAL latfi(ngridmx), lonfi(ngridmx)
!      REAL champhys(ngridmx)
      REAL tsurf(ngridmx)
      INTEGER length
      PARAMETER (length=100)
      REAL tab_cntrl(length)

c

!      EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
!      EXTERNAL exner_hyb , SSUM
c
#include "serre.h"
#include "clesph0.h"
#include "fxyprim.h"
#include "surfdat.h"
!#include "comsoil.h"
!#include "planete.h"
#include "comcstfi.h"

      real tsoil(ngridmx,nsoil),emis(ngridmx)
      real q2(ngridmx, llm+1),qsurf(ngridmx,nq)
      real airefi(ngridmx)
      real alb(ngridmx),ith(ngridmx,nsoil)
      real pzmea(ngridmx),pzstd(ngridmx)
      real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
      integer ig
      INTEGER lnblnk
      EXTERNAL lnblnk

! flag which identifies if we are using old tracer names (qsurf01,...)
      integer :: count
      character(len=30) :: txt ! to store some text
c-----------------------------------------------------------------------

      ! copy airefi(:) to area(:)
      CALL SCOPY(ngridmx,airefi,1,area,1)
      ! note: area() is defined in comgeomfi.h

      DO ig=1,ngridmx
         albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
         zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
         zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat.h
         zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat.h
         zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat.h
         zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat.h
         inertiedat(ig,:)=ith(ig,:) ! note inertiedat() is defined in comsoil.h
      ENDDO

c
c  things to store in the physics start file:
c
      ierr = NF_CREATE(adjustl(filename),NF_CLOBBER, nid)
      IF (ierr.NE.NF_NOERR) THEN
        WRITE(6,*)'physdem1: Problem creating file ',adjustl(filename)
        write(6,*) NF_STRERROR(ierr)
        CALL ABORT
      ENDIF
c
      print*,'physdem1d creating restart'

      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
     .                       "Physics start file")
c
      ierr = NF_DEF_DIM (nid,"index",length,idim1)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining index dimension'
	write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"physical_points",ngridmx,idim2)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining physical_points dimension'
	write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
      if (ierr.ne.NF_NOERR) then
      WRITE(6,*)'physdem1: Problem defining subsurface_layers dimension'
	write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
!      ierr = NF_DEF_DIM (nid,"nlayer+1",llm+1,idim4)
      ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining nlayer+1 dimension'
        write(6,*) NF_STRERROR(ierr)
	call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining advected fields dimension'
        WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
        write(6,*) NF_STRERROR(ierr)
      endif

      ierr = NF_ENDDEF(nid) ! exit NetCDF define mode

c clear tab_cntrl(:) array
      DO ierr = 1, length
         tab_cntrl(ierr) = 0.0
      ENDDO

      write(*,*) "physdem1: ngridmx: ",ngridmx
      write(*,*) "change of eccentricity, ecc=",eccpal
      peripal=(periheli+aphelie)/2.*(1-eccpal)
      aphepal=(periheli+aphelie)/2.*(1+eccpal)
      write(*,*) "check peri old/new=",periheli,peripal
      write(*,*) "check aphe old/new=",aphelie,aphepal
      write(*,*) "check a old/new=",aphelie+periheli,aphepal+peripal

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Fill control array tab_cntrl(:) with paramleters for this run
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Informations on the physics grid
      tab_cntrl(1) = float(ngridmx)  ! number of nodes on physics grid
      tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
      tab_cntrl(3) = day_ini + int(time)         ! initial day 
      tab_cntrl(4) = time -int(time)            ! initiale time of day

c Informations about Pluto, used by dynamics and physics
      tab_cntrl(5) = rad      ! radius (m)
      tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
      tab_cntrl(7) = g        ! gravity (m.s-2) 
      tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) 
      tab_cntrl(9) = rcp      !  = r/cp (=kappa dans dynamique)
      tab_cntrl(10) = daysec  ! length of a sol (s) 

      tab_cntrl(11) = phystep  ! time step in the physics
      tab_cntrl(12) = 0.
      tab_cntrl(13) = 0.

c Informations about Pluto, only for physics
      tab_cntrl(14) = year_day  ! length of year (sols) 
      tab_cntrl(15) = peripal  ! min. Sun-Pluto distance (Mkm) 
      tab_cntrl(16) = aphepal   ! max. SUn-Pluto distance (Mkm) 
      tab_cntrl(17) = peri_daypal  ! date of perihelion (sols since N. spring)
      tab_cntrl(18) = oblipal  ! Obliquity of the dwarf planet (deg) 

c Boundary layer and turbulence
      tab_cntrl(19) = z0        ! surface roughness (m) 

c Paleoclimate run
      tab_cntrl(20) = tpalnew  ! Time to start paleoclimate run (e.g -10)
      tab_cntrl(21) = adjustnew  ! Change in N2 ice albedo 
c aerosol properties

      tab_cntrl(28) = 0. 
      tab_cntrl(29) = 0.
      tab_cntrl(30) = 0.

! Soil properties:
      tab_cntrl(35) = volcapa ! soil volumetric heat capacity
      

!      write(*,*) "physdem1: tab_cntrl():",tab_cntrl
      
      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
         CALL abort
      ENDIF
      ! define variable
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to define controle'
         CALL abort
      ENDIF
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
     .                        "Control parameters")
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to define controle title attribute'
         CALL abort
      ENDIF
      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to swich out of NetCDF define mode'
         CALL abort
      ENDIF
      ! write variable
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to write controle data'
         CALL abort
      ENDIF

! write mid-layer depths mlayer() !known from comsoil.h

      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
      ! define variable
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
#else
      ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
     .                        "Soil mid-layer depth")
      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
      ! write variable
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
     .                        "Longitudes of physics grid")
      ierr = NF_ENDDEF(nid)

#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
     .                        "Latitudes of physics grid")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
     .                        "Mesh area")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
     .                        "Geopotential at the surface")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfipal)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfipal)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
     .                        "Albedo of bare ground")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
#endif

c
c   some data for Francois Lott programs
c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
     .                        "Relief: mean relief")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
     .                        "Relief: standard deviation")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: sigma parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: gamma parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: theta parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
#endif

c Write the physical fields

! Soil Thermal inertia

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
     &                   2,(/idim2,idim3/),nvarid)
#else
      ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
     &                   2,(/idim2,idim3/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
     .                        "Soil thermal inertia")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
#endif

!  Surface temperature

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
     .                        "Surface temperature")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsurf)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsurf)
#endif

! Soil temperature

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"tsoil",NF_DOUBLE,2,(/idim2,idim3/),nvarid)
#else
!      ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 2, idim2,nvarid)
      ierr = NF_DEF_VAR (nid,"tsoil",NF_FLOAT,2,(/idim2,idim3/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
     .                        "Soil temperature")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil)
#endif

c emissivity

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "emis", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
     .                        "Surface emissivity")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,emis)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
#endif

c planetary boundary layer

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR(nid,"q2",NF_DOUBLE, 2,(/idim2,idim4/),nvarid)
#else
      ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 2,(/idim2,idim4/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
     .                        "pbl wind variance")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q2)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,q2)
#endif
      IF (ierr.NE.NF_NOERR) THEN
        PRINT*, 'physdem1: Failed to write q2'
        CALL abort
      ENDIF

c tracers

! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
!                    qsurf02, ...)
      count=0
      do iq=1,nqmx
        txt= " "
        write(txt,'(a1,i2.2)')'q',iq
        if (txt.ne.tnom(iq)) then ! use tracer names stored in dynamics
          ! did not find old tracer name
          exit ! might as well stop here
        else
          ! found old tracer name
          count=count+1
        endif
      enddo
      if (count.eq.nqmx) then
        write(*,*) "physdem1:tracers seem to follow old naming ",
     &             "convention (qsurf01,qsurf02,...)"

        call abort
      endif ! of if (count.eq.nqmx)

      IF(nq.GE.1) THEN
         DO iq=1,nq
           txt=tnom(iq)

           ierr=NF_REDEF(nid)
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
             CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,1,idim2,nvarid)
#else
           ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,1,idim2,nvarid)
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to define ',trim(txt)
             CALL abort
           ENDIF
           ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
     &                        "tracer on surface")
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to define ',trim(txt),
     &               ' title attribute'
             CALL abort
           ENDIF
           ierr=NF_ENDDEF(nid)
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to swich out of define mode'
             CALL abort
           ENDIF
           
#ifdef NC_DOUBLE
            ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,qsurf(1,iq))
#else
            ierr=NF_PUT_VAR_REAL (nid,nvarid,qsurf(1,iq))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to write ',trim(txt)
             CALL abort
           ENDIF
         ENDDO ! of DO iq=1,nq
      ENDIF ! of IF(nq.GE.1)

c close file
      ierr = NF_CLOSE(nid)

      RETURN

      END