Changeset 2432 for trunk/LMDZ.MARS/util/lslin.F90

Timestamp:

Nov 14, 2020, 3:27:44 PM (4 years ago)

Author:

emillour

Message:

Mars GCM:

• Make a single README file (get rid of REAME.exec) and update "compile" script to (hopefully) better illustrate how to compile utilities
• Update utilities zrecast and lslin: zrecast.F90: Force the vertical interpolation of any variable starting by "rho" to be in log space (as it should for density, molecular concentration in mol.cm-3 etc..). Set missing value to 1.E20 lslin.F90: fix various problems in the output.

FF + EM

File:

• trunk/LMDZ.MARS/util/lslin.F90 (modified) (26 diffs)

trunk/LMDZ.MARS/util/lslin.F90

@@ -27 +27 @@
 character (len=50) :: units_alt
 ! var(): units of altitude which can be 'km' or 'Pa' (after zrecast)
+character (len=50) :: altlong_name,altunits,altpositive
+! altlong_name(): [netcdf] altitude "long_name" attribute
+! altunits(): [netcdf] altitude "units" attribute
+! altpositive(): [netcdf] altitude "positive" attribute
+
 character (len=100) :: outfile
 ! outfile(): output file name
@@ -35 +40 @@
 character (len=1) :: average
 ! average: the character 'y' to average within the Ls timestep
-integer :: nid,varid,ierr,miss
+integer :: nid,varid,ierr,miss,validr
 ! nid: [netcdf] ID # of input file
 ! varid: [netcdf] ID # of a variable
@@ -55 +60 @@
 ! tab_cntrl(length): array, stored in the input file,
 !                which contains various control parameters.
-real, dimension(:), allocatable:: lat,lon,alt,time,lsgcm,timels
+real, dimension(:), allocatable:: lat,lon,alt,time,lsgcm,timels,ctl
 ! lat(): latitude coordinates (read from input file)
 ! lon(): longitude coordinates (read from input file)
@@ -62 +67 @@
 ! lsgcm(): time coordinate (in unevenly spaced "Ls")
 ! timels(): new time coordinates (evenly spaced "Ls"; written to output file)
-integer :: latlen,lonlen,altlen,timelen,Nls,Sls
+! ctl(): array, stores controle array
+integer :: latlen,lonlen,altlen,timelen,Nls,Sls,ctllen
 ! latlen: # of elements of lat() array
 ! lonlen: # of elements of lon() array
@@ -72 +78 @@
 ! nbvarfile: total # of variables (in netcdf file)
 ! ndim: [netcdf] # of dimensions (3 or 4) of a variable
-integer :: latdim,londim,altdim,timedim
+integer :: latdim,londim,altdim,timedim,ctldim
 ! latdim: [netcdf] "latitude" dim ID
 ! londim: [netcdf] "longitude" dim ID
 ! altdim: [netcdf] "altdim" dim ID
 ! timedim: [netcdf] "timedim" dim ID
-integer :: latvar,lonvar,altvar,timevar
+! ctldim: [netcdf] "controle" dim ID
+integer :: latvar,lonvar,altvar,timevar,ctlvar
 ! latvar: [netcdf] ID of "latitude" variable
 ! lonvar: [netcdf] ID of "longitude" variable
@@ -240 +247 @@
    else
       ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
-      units_alt="                                                    "
+      units_alt=""
       ierr=NF_GET_ATT_TEXT(nid,varid,"units",units_alt)
    endif
    write(*,*) "altlen: ",altlen
+! Get altitude attributes to handle files with any altitude type
+   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
+   ierr=nf_get_att_text(nid,altvar,'units',altunits)
+   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
+
 
 ! Allocate lat(), lon() and alt()
@@ -251 +263 @@
 
 ! Read lat(),lon() and alt() from input file
-
-
-
-
-
       ierr = NF_GET_VAR_REAL(nid,latvar,lat)
       ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
@@ -264 +271 @@
 
 
+   ierr=NF_INQ_DIMID(nid,"index",ctldim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) ' Dimension <index> is missing in file '//trim(infile)
+      ctllen=0
+      !stop ""  
+   endif
+   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'Field <controle> is missing in file '//trim(infile)
+      ctllen=0
+      !stop ""
+   else
+      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
+   endif
+
+   allocate(ctl(ctllen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate ctl(ctllen)"
+     stop
+   endif
+
+   if (ctllen .ne. 0) then
+     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+     if (ierr.ne.0) then
+       write(*,*) "Error: failed to load controle"
+       write(*,*) NF_STRERROR(ierr)
+       stop
+     endif
+   endif ! of if (ctllen .ne. 0)
+
+
+
 !==============================================================================
 ! 2.2. Create (and initialize latitude, longitude and altitude in) output file
 !==============================================================================
 
-call initiate(outfile,lat,lon,alt,nout,&
-           latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
+  ! Initialize output file's lat,lon,alt and time dimensions
+   call initiate(outfile,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+         altlong_name,altunits,altpositive,&
+         nout,latdimout,londimout,altdimout,timedimout,timevarout)
+
+  ! Initialize output file's aps,bps and phisinit variables
+   call init2(nid,lonlen,latlen,altlen,altdim,&
+              nout,londimout,latdimout,altdimout)
 
 !==============================================================================
@@ -328 +373 @@
    endif
 
-   day_ini = tab_cntrl(4)                                               
+   day_ini = nint(tab_cntrl(4))                                               
    day_ini = modulo(day_ini,669)                                                    
    write(*,*) 'day_ini', day_ini
@@ -392 +437 @@
 
 ! Build timels()
-timels(1) = 0.01*nint(100.*ls0) ! Ehouarn: what the !!!???!!!
+timels(1) = 0.01*nint(100.*ls0) 
 do k=2,Nls
    timels(k) = timels(k-1) + deltals
@@ -482 +527 @@
 ! 3.3 Write this variable's definition and attributes to the output file
 !==============================================================================
-   units="                                                    "
-   title="                                                    "
-   ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
-   ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
-
-   call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
+   if (ndim.ge.3) then
+    units=""
+    title=""
+    ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
+    ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
+
+    call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
+   end if
 
 !==============================================================================
 ! 3.4 Read variable
-!==============================================================================
-
-
-
-
-   ierr = NF_GET_VAR_REAL(nid,varid,var3d)
-   miss = NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
-   miss = NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
+!==== ==========================================================================
+
+   if (ndim.ge.3) then
+      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
+      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
+      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
+   end if
 
 
@@ -582 +628 @@
 
 
-
-   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
-
-
-   if (ierr.ne.NF_NOERR) then
+   if (ndim.ge.3) then
+
+     ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
+
+
+     if (ierr.ne.NF_NOERR) then
       write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
       stop ""
-   endif
+     endif
 
 ! In case there is a "missing value" attribute and "valid range"
-   if (miss.eq.NF_NOERR) then
-      call missing_value(nout,varidout,valid_range,missing)  
-   end if
-
-   deallocate(var3d)
-   deallocate(var3dls)
-   deallocate(var3dxy)
+     if (miss.eq.NF_NOERR) then
+      call missing_value(nout,varidout,missing)  
+     end if
+
+     deallocate(var3d)
+     deallocate(var3dls)
+     deallocate(var3dxy)
+   endif !  if (ndim.ge.3) 
 
 enddo ! of do j=1,nbvar loop
@@ -641 +689 @@
 contains
 
-!************************************
-subroutine initiate (filename,lat,lon,alt,&
-                     nout,latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
+!******************************************************************************
+
+subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+                     altlong_name,altunits,altpositive,&
+                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
 !==============================================================================
 ! Purpose:
@@ -661 +711 @@
 character (len=*), intent(in):: filename
 ! filename(): the file's name
-real, dimension(:), intent(in):: lat
+integer,intent(in) ::latlen,lonlen,altlen,ctllen
+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, intent(in):: ctl(ctllen)
+! ctl(): controle
+character (len=*), intent(in) :: altlong_name
+! altlong_name(): [netcdf] altitude "long_name" attribute
+character (len=*), intent(in) :: altunits
+! altunits(): [netcdf] altitude "units" attribute
+character (len=*), intent(in) :: altpositive
+! altpositive(): [netcdf] altitude "positive" attribute
 integer, intent(out):: nout
 ! nout: [netcdf] file ID
@@ -679 +738 @@
 integer, intent(out):: timevarout
 ! timevarout: [netcdf] "Time" (considered as a variable) ID
-character (len=50) :: units_alt
-! var(): units of altitude which can be m or Pa (after zrecast)
-
 
 !==============================================================================
@@ -688 +744 @@
 !integer :: latdim,londim,altdim,timedim
 integer :: nvarid,ierr
+integer :: ctldimout
 ! nvarid: [netcdf] ID of a variable
 ! ierr: [netcdf]  return error code (from called subroutines)
@@ -698 +755 @@
 ! 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
@@ -706 +764 @@
 !==============================================================================
 
-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)
+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
 
 !==============================================================================
@@ -718 +809 @@
 !==============================================================================
 
-call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
+call def_var(nout,"Time","Ls (Solar Longitude)","degrees",1,&
              (/timedimout/),timevarout,ierr)
+! Switch to netcdf define mode
+ierr = NF_REDEF (nout)
+ierr = NF_ENDDEF(nout)
 
 !==============================================================================
@@ -729 +823 @@
 
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <latitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -738 +837 @@
 
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
-
-!==============================================================================
-! 4. Write "altitude" (data and attributes)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <longitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 5. Write "altitude" (data and attributes)
 !==============================================================================
 
@@ -748 +852 @@
 ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
 
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"altitude")
-!********************* temlporaire ****************
-!ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,units_alt)
-!if((units_alt(1:2).eq.'Pa').or.(units_alt(1:2).eq.'pa')) then
-!   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"down")
-!else
-!   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
-!end if
-!********************* temlporaire ****************
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,'Pa')
-   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',4,'down')
-!********************* temlporaire ****************
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
 
 ! End netcdf define mode
@@ -765 +860 @@
 
 ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <altitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 6. Write "controle" (data and attributes)
+!==============================================================================
+
+if (size(ctl).ne.0) then
+   ! Switch to netcdf define mode
+   ierr = NF_REDEF (nout)
+
+   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
+
+   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
+
+   ! End netcdf define mode
+   ierr = NF_ENDDEF(nout)
+
+   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
+   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,altdimid, &
+                 outfid,londimout,latdimout,altdimout)
+!==============================================================================
+! Purpose:
+! Copy aps(), bps() and phisinit() from input file to output file
+!==============================================================================
+! Remarks:
+! The NetCDF files must be open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: infid  ! NetCDF output file ID
+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 altitude
+integer, intent(in) :: altdimid ! ID of altitude dimension
+integer, intent(in) :: outfid ! NetCDF output file ID
+integer, intent(in) :: londimout ! longitude dimension ID
+integer, intent(in) :: latdimout ! latitude dimension ID
+integer, intent(in) :: altdimout ! altitude dimension ID
+!==============================================================================
+! Local variables:
+!==============================================================================
+real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
+real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
+integer :: apsid,bpsid,phisinitid
+integer :: ierr
+integer :: tmpdimid ! temporary dimension ID
+integer :: tmpvarid ! temporary variable ID
+integer :: tmplen ! temporary variable length
+logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
+
+
+!==============================================================================
+! 1. Read data from input file
+!==============================================================================
+
+! hybrid coordinate aps
+  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
+  write(*,*) "oops Failed to get aps ID. OK"
+  aps_ok=.false.
+else
+  ! 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)
+  ierr=NF_INQ_DIMLEN(infid,tmpdimid,tmplen)
+  if (tmplen.ne.altlen) then
+    write(*,*) "tmplen,altlen", tmpdimid, altdimid
+    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"
+     aps_ok=.false.
+    endif
+    aps_ok=.true.
+  endif
+endif
+
+! hybrid coordinate bps
+  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
+  write(*,*) "oops: Failed to get bps ID. OK"
+  bps_ok=.false.
+else
+  ! Check that bps() is the right size
+  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
+  ierr=NF_INQ_DIMLEN(infid,tmpdimid,tmplen)
+  if (tmplen.ne.altlen) 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"
+      bps_ok=.false.
+    endif
+    bps_ok=.true.
+  endif
+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)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Failed to get phisinit ID. OK"
+  phisinit = 0.
+  phis = .false.
+else
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading phisinit"
+  endif
+  phis = .true.
+endif
+
+
+!==============================================================================
+! 2. Write
+!==============================================================================
+
+!==============================================================================
+! 2.2. Hybrid coordinates aps() and bps()
+!==============================================================================
+
+IF (aps_ok) then 
+  ! 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
+  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
+  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(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
+  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write bps"
+  endif
+ENDIF ! of IF (bps_ok)
+
+!==============================================================================
+! 2.2. phisinit()
+!==============================================================================
+
+IF (phis) 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
+  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write phisinit"
+  endif
+
+ENDIF ! of IF (phis)
+
+
+! Cleanup
+deallocate(aps)
+deallocate(bps)
+deallocate(phisinit)
+
+end subroutine init2
+
+!******************************************************************************
 subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
 !==============================================================================
@@ -816 +1129 @@
 
 end subroutine def_var
-!************************************
+
+!******************************************************************************
+subroutine  missing_value(nout,nvarid,missing)
+!==============================================================================
+! Purpose:
+! Write "valid_range" and "missing_value" attributes (of a given
+! variable) to a netcdf file
+!==============================================================================
+! Remarks:
+! NetCDF file must be open
+! Variable (nvarid) ID must be set
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc"  ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: nout
+! nout: [netcdf] file ID #
+integer, intent(in) :: nvarid
+! varid: [netcdf] variable ID #
+!real, dimension(2), intent(in) :: valid_range
+! valid_range(2): [netcdf] "valid_range" attribute (min and max)
+real, intent(in) :: missing
+! missing: [netcdf] "missing_value" attribute
+
+!==============================================================================
+! Local variables:
+!==============================================================================
+integer :: ierr
+! ierr: [netcdf] subroutine returned error code
+!      INTEGER nout,nvarid,ierr
+!      REAL missing
+!      REAL valid_range(2)
+
+! Switch to netcdf dataset define mode
+ierr = NF_REDEF (nout)
+if (ierr.ne.NF_NOERR) then
+   print*,'missing_value: '
+   print*, NF_STRERROR(ierr)
+endif
+
+
+!********* valid range not used in Lslin ****************
+! Write valid_range() attribute
+!ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
+!
+!if (ierr.ne.NF_NOERR) then
+!   print*,'missing_value: valid_range attribution failed'
+!   print*, NF_STRERROR(ierr)
+!   !write(*,*) 'NF_NOERR', NF_NOERR
+!   !CALL abort
+!   stop ""
+!endif
+!*********************************************************
+
+! Write "missing_value" attribute
+ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
+if (ierr.NE.NF_NOERR) then
+   print*, 'missing_value: missing value attribution failed'
+   print*, NF_STRERROR(ierr)
+!    WRITE(*,*) 'NF_NOERR', NF_NOERR
+!    CALL abort
+   stop ""
+endif
+
+! End netcdf dataset define mode
+ierr = NF_ENDDEF(nout)
+
+end subroutine  missing_value
+
+!******************************************************************************
 subroutine sol2ls(sol,Ls)
 !==============================================================================
@@ -921 +1308 @@
 integer nd
 !  internal
-integer i,j
+integer i
 real y_undefined
                                                                                 
@@ -955 +1342 @@
 end subroutine interpolf
      
-end
-
-
-
-!******************************************************************************
-!******************************************************************************
-subroutine  missing_value(nout,nvarid,valid_range,missing)
-!==============================================================================
-! Purpose:
-! Write "valid_range" and "missing_value" attributes (of a given
-! variable) to a netcdf file
-!==============================================================================
-! Remarks:
-! NetCDF file must be open
-! Variable (nvarid) ID must be set
-!==============================================================================
-
-implicit none
-
-include "netcdf.inc"  ! NetCDF definitions
-
-!==============================================================================
-! Arguments:
-!==============================================================================
-integer, intent(in) :: nout
-! nout: [netcdf] file ID #
-integer, intent(in) :: nvarid
-! varid: [netcdf] variable ID #
-real, dimension(2), intent(in) :: valid_range
-! valid_range(2): [netcdf] "valid_range" attribute (min and max)
-real, intent(in) :: missing
-! missing: [netcdf] "missing_value" attribute
-
-!==============================================================================
-! Local variables:
-!==============================================================================
-integer :: ierr
-! ierr: [netcdf] subroutine returned error code
-!      INTEGER nout,nvarid,ierr
-!      REAL missing
-!      REAL valid_range(2)
-
-! Switch to netcdf dataset define mode
-ierr = NF_REDEF (nout)
-if (ierr.ne.NF_NOERR) then
-   print*,'missing_value: '
-   print*, NF_STRERROR(ierr)
-endif
-
-! Write valid_range() attribute
-#ifdef NC_DOUBLE
-ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
-#else
-ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
-#endif
-
-if (ierr.ne.NF_NOERR) then
-   print*,'missing_value: valid_range attribution failed'
-   print*, NF_STRERROR(ierr)
-   !write(*,*) 'NF_NOERR', NF_NOERR
-   !CALL abort
-   stop ""
-endif
-
-! Write "missing_value" attribute
-#ifdef NC_DOUBLE
-ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
-#else
-ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
-#endif
-
-if (ierr.NE.NF_NOERR) then
-   print*, 'missing_value: missing value attribution failed'
-   print*, NF_STRERROR(ierr)
-!    WRITE(*,*) 'NF_NOERR', NF_NOERR
-!    CALL abort
-   stop ""
-endif
-
-! End netcdf dataset define mode
-ierr = NF_ENDDEF(nout)
-
-end subroutine  missing_value
-!******************************************************************************
+end program lslin