Index: trunk/LMDZ.MARS/util/localtime.F90
===================================================================
--- trunk/LMDZ.MARS/util/localtime.F90	(revision 2202)
+++ trunk/LMDZ.MARS/util/localtime.F90	(revision 2204)
@@ -115,10 +115,11 @@
 
 !==============================================================================
-! 1.3. Open the first input file
+! 1.3. Open the input file
 !==============================================================================
 
 ierr = NF_OPEN(file,NF_NOWRITE,nid)
 if (ierr.NE.NF_NOERR) then
-   write(*,*) 'ERROR: Pb opening file '//file
+   write(*,*) 'ERROR: Pb opening file '//trim(file)
+   write(*,*) NF_STRERROR(ierr)
    stop ""
 endif
@@ -126,7 +127,12 @@
 ierr=NF_INQ_NVARS(nid,nbvarfile)
 ! nbvarfile now set to be the (total) number of variables in file
-
-!==============================================================================
-! 1.4. Ask for (output) "Time" axis type
+if (ierr.NE.NF_NOERR) then
+   write(*,*) 'ERROR: Pb with NF_INQ_NVARS'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 1.4. List of variables that should not be processed
 !==============================================================================
 
@@ -148,5 +154,5 @@
 
 !==============================================================================
-! 1.5. Get (and check) list of variables to concatenate
+! 1.5. Get (and check) list of variables to process
 !==============================================================================
 write(*,*)
@@ -165,7 +171,13 @@
 enddo
 
-! Nnotconcat: # of variables that won't be concatenated
+! Nnotconcat: # of variables that won't be processed
 ! nbvarfile: total # of variables in file
-allocate(var(nbvarfile-Nnotconcat))
+allocate(var(nbvarfile-Nnotconcat),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
+  write(*,*) "  nbvarfile=",nbvarfile
+  write(*,*) "  Nnotconcat=",Nnotconcat
+  stop
+endif
 
 
@@ -187,5 +199,5 @@
             ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
          enddo
-! Variables names from the file are catched
+! Variables names from the file are stored in var()
    nbvar=nbvarfile-Nnotconcat
    do i=1,nbvar
@@ -202,5 +214,5 @@
 !==============================================================================
  filename=file(1:len_trim(file)-3)//"_LT.nc"
- write(*,*) filename
+ write(*,*) trim(filename)
 
 !==============================================================================
@@ -208,14 +220,12 @@
 !==============================================================================
 
-   if (i/=1) then
       write(*,*) 
       write(*,*) "opening "//trim(file)//"..."
       ierr = NF_OPEN(file,NF_NOWRITE,nid)
       if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Pb opening file '//file
+         write(*,*) 'ERROR: Pb opening file '//trim(file)
          write(*,*) NF_STRERROR(ierr)
          stop ""
       endif
-   endif
 
 !==============================================================================
@@ -224,7 +234,11 @@
 
    ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
+      stop ""  
+   endif
    ierr=NF_INQ_VARID(nid,"latitude",latvar)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <latitude> is missing in file'//file
+      write(*,*) 'ERROR: Field <latitude> is missing in file '//trim(file)
       stop ""  
    endif
@@ -233,7 +247,11 @@
 
    ierr=NF_INQ_DIMID(nid,"longitude",londim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
+      stop ""  
+   endif
    ierr=NF_INQ_VARID(nid,"longitude",lonvar)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <longitude> is missing in file'//file
+      write(*,*) 'ERROR: Field <longitude> is missing in file'//trim(file)
       stop "" 
    endif
@@ -242,7 +260,11 @@
 
    ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
+      stop ""  
+   endif
    ierr=NF_INQ_VARID(nid,"altitude",altvar)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <altitude> is missing in file'//file
+      write(*,*) 'ERROR: Field <altitude> is missing in file'//trim(file)
       stop ""
    endif
@@ -251,7 +273,11 @@
 
    ierr=NF_INQ_DIMID(nid,"index",ctldim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <index> is missing in file '//trim(file)
+      stop ""  
+   endif
    ierr=NF_INQ_VARID(nid,"controle",ctlvar)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'Field <controle> is missing in file'//file
+      write(*,*) 'Field <controle> is missing in file'//trim(file)
       ctllen=0
       !stop ""
@@ -267,19 +293,71 @@
 
 ! First call; initialize/allocate
-      allocate(lat(latlen))
-      allocate(lon(lonlen))
-      allocate(alt(altlen))
-      allocate(ctl(ctllen))
+      allocate(lat(latlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lat(latlen)"
+        stop
+      endif
+      allocate(lon(lonlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lon(lonlen)"
+        stop
+      endif
+      allocate(alt(altlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate alt(altlen)"
+        stop
+      endif
+      allocate(ctl(ctllen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate ctl(ctllen)"
+        stop
+      endif
+
 #ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
+#else
+      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load latitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+
+#ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
+#else
+      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load longitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+
+#ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
-      if (ctllen .ne. -1) ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
-#else
-      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
-      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+#else
       ierr = NF_GET_VAR_REAL(nid,altvar,alt)
-      if (ctllen .ne. -1) ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
-#endif
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load altitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+
+      if (ctllen .ne. 0) then
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
+#else
+        ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+#endif
+        if (ierr.ne.0) then
+          write(*,*) "Error: failed to load controle"
+          write(*,*) NF_STRERROR(ierr)
+          stop
+        endif
+      endif ! of if (ctllen .ne. 0)
+
 !==============================================================================
 ! 2.4. Handle "Time" dimension from input file
@@ -291,10 +369,10 @@
    ierr=NF_INQ_DIMID(nid,"Time",timedim)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <Time> is missing in file'//file
+      write(*,*) 'ERROR: Dimension <Time> is missing in file'//trim(file)
       stop ""
    endif
    ierr=NF_INQ_VARID(nid,"Time",timevar)
    if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <Time> is missing in file'//file
+      write(*,*) 'ERROR: Field <Time> is missing in file'//trim(file)
       stop ""
    endif
@@ -313,5 +391,9 @@
 
    ! allocate time() array and fill it with values from input file
-   allocate(time(timelen))
+   allocate(time(timelen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate time(timelen)"
+     stop
+   endif
 
 #ifdef NC_DOUBLE
@@ -322,4 +404,5 @@
    if (ierr.NE.NF_NOERR) then
      write(*,*) "Error , failed to load Time"
+     write(*,*) NF_STRERROR(ierr)
      stop
    endif
@@ -331,5 +414,9 @@
       write(*,*) 'number of local time to be used ?'
       read(*,*) nhour
-      allocate(lt_hour(nhour))
+      allocate(lt_hour(nhour),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_hour(nhour)"
+        stop
+      endif
       write(*,*) 'list of selected local time (0<t<24)'
       do it=1,nhour
@@ -341,5 +428,9 @@
          stats=.true.
          timelen_lt=nhour
-         allocate(lt_out(timelen_lt))
+         allocate(lt_out(timelen_lt),stat=ierr)
+         if (ierr.ne.0) then
+           write(*,*) "Error: failed to allocate lt_hour(nhour)"
+           stop
+         endif
          do it=1,nhour
            lt_out(it)=lt_hour(it)/24.
@@ -355,5 +446,9 @@
         nsol = int(time(timelen)) - int(time(1))
         timelen_lt=nhour*nsol
-        allocate(lt_out(timelen_lt))
+        allocate(lt_out(timelen_lt),stat=ierr)
+        if (ierr.ne.0) then
+           write(*,*) "Error: failed to allocate lt_hour(nhour)"
+           stop
+         endif
         i=0
         do k=1,nsol
@@ -371,9 +466,19 @@
          timelen_tot=timelen+1
       endif      
-      allocate(lt_gcm(timelen_tot))
-      allocate(var_gcm(timelen_tot))
-
-      allocate(lt_outc(timelen_lt))
-         
+      allocate(lt_gcm(timelen_tot),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
+        stop
+      endif
+      allocate(var_gcm(timelen_tot),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
+        stop
+      endif
+      allocate(lt_outc(timelen_lt),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_outc(timelen_lt)"
+        stop
+      endif
 
 !==============================================================================
@@ -383,8 +488,9 @@
 
    ! Initialize output file's lat,lon,alt and time dimensions
-      call initiate (filename,lat,lon,alt,ctl,nout,&
-           latdimout,londimout,altdimout,timedimout,timevarout)
+     call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+           nout,latdimout,londimout,altdimout,timedimout,timevarout)
+
    ! Initialize output file's aps,bps and phisinit variables
-     call init2(nid,lonlen,latlen,altlen,&
+     call init2(nid,lonlen,latlen,altlen,altdim,&
                 nout,londimout,latdimout,altdimout)
 
@@ -454,5 +560,5 @@
          corner(4)=1
 
-	 ! length (along dimensions) of block of data to be written
+         ! length (along dimensions) of block of data to be written
          edges(1)=lonlen 
          edges(2)=latlen 
@@ -474,5 +580,5 @@
          corner(4)=1
 
-	 ! length (along dimensions) of block of data to be written
+         ! length (along dimensions) of block of data to be written
          edges(1)=lonlen
          edges(2)=latlen
@@ -481,6 +587,6 @@
       endif
 
-         units="                                                    "
-         title="                                                    "
+         units=" "
+         title=" "
          ierr=nf_get_att_text(nid,varid,"title",title)
          ierr=nf_get_att_text(nid,varid,"units",units)
@@ -597,5 +703,5 @@
 
 !******************************************************************************
-Subroutine initiate (filename,lat,lon,alt,ctl,&
+Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
                      nout,latdimout,londimout,altdimout,timedimout,timevarout)
 !==============================================================================
@@ -616,12 +722,13 @@
 character (len=*), intent(in):: filename
 ! filename(): the file's name
-real, dimension(:), intent(in):: lat
+real, intent(in):: lat(latlen)
 ! lat(): latitude
-real, dimension(:), intent(in):: lon
+real, intent(in):: lon(lonlen)
 ! lon(): longitude
-real, dimension(:), intent(in):: alt
+real, intent(in):: alt(altlen)
 ! alt(): altitude
-real, dimension(:), intent(in):: ctl
+real, intent(in):: ctl(ctllen)
 ! ctl(): controle
+integer,intent(in) ::latlen,lonlen,altlen,ctllen
 integer, intent(out):: nout
 ! nout: [netcdf] file ID
@@ -653,5 +760,6 @@
 ! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
 if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'ERROR: Impossible to create the file.'
+   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
+   write(*,*) NF_STRERROR(ierr)
    stop ""
 endif
@@ -661,12 +769,44 @@
 !==============================================================================
 
-ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
-ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
-ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
-if (size(ctl).ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
+ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+if (size(ctl).ne.0) then
+  ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
+  if (ierr.NE.NF_NOERR) then
+    WRITE(*,*)'initiate: error failed to define dimension <index>'
+    write(*,*) NF_STRERROR(ierr)
+    stop ""
+  endif
+endif
 ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <Time>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 ! End netcdf define mode
 ierr = NF_ENDDEF(nout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to switch out of define mode'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -689,4 +829,9 @@
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
 #endif
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <latitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -702,4 +847,9 @@
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
 #endif
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <longitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -728,4 +878,9 @@
 ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
 #endif 
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <altitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -753,9 +908,14 @@
    ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
 #endif
+   if (ierr.NE.NF_NOERR) then
+      WRITE(*,*)'initiate: error failed writing variable <controle>'
+      write(*,*) NF_STRERROR(ierr)
+      stop ""
+   endif
 endif
 
 end Subroutine initiate
 !******************************************************************************
-subroutine init2(infid,lonlen,latlen,altlen, &
+subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
                  outfid,londimout,latdimout,altdimout)
 !==============================================================================
@@ -777,5 +937,6 @@
 integer, intent(in) :: lonlen ! # of grid points along longitude
 integer, intent(in) :: latlen ! # of grid points along latitude
-integer, intent(in) :: altlen ! # of grid points along latitude
+integer, intent(in) :: altlen ! # of grid points along altitude
+integer, intent(in) :: altdimid ! ID of altitude dimension
 integer, intent(in) :: outfid ! NetCDF output file ID
 integer, intent(in) :: londimout ! longitude dimension ID
@@ -789,4 +950,5 @@
 integer :: apsid,bpsid,phisinitid
 integer :: ierr
+integer :: tmpdimid ! temporary dimension ID
 integer :: tmpvarid ! temporary variable ID
 logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
@@ -798,5 +960,10 @@
 
 ! hybrid coordinate aps
-  allocate(aps(altlen))
+  allocate(aps(altlen),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate aps(altlen)"
+    stop
+  endif
+
 ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
 if (ierr.ne.NF_NOERR) then
@@ -804,13 +971,26 @@
   aps_ok=.false.
 else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
-  if (ierr.ne.NF_NOERR) then
-   stop "error: Failed reading aps"
+  ! Check that aps() is the right size (it most likely won't be if
+  ! zrecast has be used to generate the input file)
+  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
+  if (tmpdimid.ne.altdimid) then
+    write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
+    aps_ok=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
+    if (ierr.ne.NF_NOERR) then
+     stop "init2 error: Failed reading aps"
+    endif
+    aps_ok=.true.
   endif
-  aps_ok=.true.
 endif
 
 ! hybrid coordinate bps
-  allocate(bps(altlen))
+  allocate(bps(altlen),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate bps(altlen)"
+    stop
+  endif
+
 ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
 if (ierr.ne.NF_NOERR) then
@@ -818,14 +998,25 @@
   bps_ok=.false.
 else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading bps"
+  ! Check that bps() is the right size
+  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
+  if (tmpdimid.ne.altdimid) then
+    write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
+    bps_ok=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
+    if (ierr.ne.NF_NOERR) then
+      stop "init2 Error: Failed reading bps"
+    endif
+    bps_ok=.true.
   endif
-  bps_ok=.true.
 endif
 
 ! ground geopotential phisinit
+allocate(phisinit(lonlen,latlen),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
+  stop
+endif
 ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
-allocate(phisinit(lonlen,latlen))
 if (ierr.ne.NF_NOERR) then
   write(*,*) "Failed to get phisinit ID. OK"
@@ -835,5 +1026,5 @@
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
   if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading phisinit"
+    stop "init2 Error: Failed reading phisinit"
   endif
   phis = .true.
@@ -850,40 +1041,40 @@
 
 IF (aps_ok) then 
-! define aps
-call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
-             (/altdimout/),apsid,ierr)
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to def_var aps"
-endif
-
-! write aps
-#ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
-#else
-ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
-#endif
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to write aps"
-endif
-ENDIF 
+  ! define aps
+  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
+               (/altdimout/),apsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to def_var aps"
+  endif
+
+  ! write aps
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write aps"
+  endif
+ENDIF ! of IF (aps_ok)
 
 IF (bps_ok) then 
-! define bps
-call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
-             (/altdimout/),bpsid,ierr)
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to def_var bps"
-endif
-
-! write bps
-#ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
-#else
-ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
-#endif
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to write bps"
-endif
-ENDIF 
+  ! define bps
+  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
+               (/altdimout/),bpsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to def_var bps"
+  endif
+
+  ! write bps
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write bps"
+  endif
+ENDIF ! of IF (bps_ok)
 
 !==============================================================================
@@ -891,25 +1082,24 @@
 !==============================================================================
 
-
 IF (phis) THEN
 
-!define phisinit
- call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
-            (/londimout,latdimout/),phisinitid,ierr)
- if (ierr.ne.NF_NOERR) then
+  !define phisinit
+  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
+              (/londimout,latdimout/),phisinitid,ierr)
+  if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to def_var phisinit"
   endif
 
-! write phisinit
-#ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
-#else
-ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
-#endif
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to write phisinit"
-endif
-
-END IF
+  ! write phisinit
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write phisinit"
+  endif
+
+ENDIF ! of IF (phis)
 
 
@@ -1051,7 +1241,4 @@
 
 end subroutine  missing_value
-!******************************************************************************
-
-end program localtime
 
 !*****************************************************************************
@@ -1096,2 +1283,6 @@
 end subroutine interpolf
 
+!******************************************************************************
+
+end program localtime
+
