source: palm/trunk/SCRIPTS/nc2vdf.ncl @ 1494

;--------------------------------------------------------------------------------!
; This file is part of PALM.
;
; PALM is free software: you can redistribute it and/or modify it under the terms
; of the GNU General Public License as published by the Free Software Foundation,
; either version 3 of the License, or (at your option) any later version.
;
; PALM is distributed in the hope that it will be useful, but WITHOUT ANY
; WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
; A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License along with
; PALM. If not, see <http://www.gnu.org/licenses/>.
;
; Copyright 1997-2014  Leibniz Universitaet Hannover
;------------------------------------------------------------------------------!
;
; Current revisions:
; -----------------
; 
; 
; Former revisions:
; -----------------
; $Id: nc2vdf.ncl 1328 2014-03-21 11:00:33Z maronga $
;
; 1310 2014-03-14 08:01:56Z raasch
; adapted for use on desktop machines
;
; 1067 2012-11-26 10:13:14Z maronga
; Code put under GNU GPL
;
; 1062 2012-11-21 15:57:18Z maronga
; full batch mode support, loading configuration from an external file
; fallback to interactive mode if no configuration present
; batch mode does not support loading more than one variable per file if there
; is more than one file specified (but there is a workaround)
; you cannot add time steps to a variable from an additional file in batch mode
;
; 1046 2012-11-09 14:38:45Z maronga
; Initial revision
;
; Description:
; ------------
; This NCL script coverts PALM NetCDF data to vdf, VAPOR's own data format. In
; order to run this script, NCL version 5.2.0 or higher is required.
; Default setting will be loaded from .nc2vdf.config.
; The script does only work on hicegate2, bicegate2,
; uv, on the imuk server. It can be adapted for use on any desktop machine.
;------------------------------------------------------------------------------!

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

alt_varnames = False
man = True
path = ""

if(isfilepresent(inputfile)) then
   print("*** Loading " + inputfile)
   loadscript(inputfile)
   man = False
else
    print("+++ Input file does not exist.")
    man_resp = systemfunc("read -p 'Enter manual mode? (y/n) ' RESP; echo $RESP")
    if (man_resp .eq. "y") then
       man = True
    else
       print("+++ Exiting...")
       exit
    end if
end if


begin

if (.not. man) then
;   print loaded data from input file
    print("*** loading " + dimsizes(files) + ":")
    do i=0,dimsizes(files)-1
       print("    " + i + ") " + files(i))
    end do

    print("*** will create dataset " + outputfile)
end if

;******************************************************
;give path to the directory where NetCDF files are stored:
    print("*** Checking hostname...")
    sys        = "unknown"
    query      = systemfunc("echo $HOSTNAME")
    check_muk  = isStrSubset(query,"muk.uni-hannover.de")
    check_hice = isStrSubset(query,"hicegate2")
    check_bice = isStrSubset(query,"bicegate2")
    check_huv  = isStrSubset(query,"huv")
    check_buv  = isStrSubset(query,"buv")




    if ( check_muk .eq. True ) then
       print("*** nc2vdf will execute at IMUK")
       sys  = "IMUK"
    else
       if ( check_hice .eq. True .or. check_bice .eq. True ) then
          print("*** nc2vdf will execute on hicegate2/bicegate2")
          sys = "HLRN"
       else
          if ( check_huv .eq. True .or. check_buv .eq. True ) then
             print("*** nc2vdf will execute on UV")
             sys = "HLRN"
             else
                if ( check_workstation .eq. True) then
                print("*** nc2vdf will execute on a workstation PC/Mac")
                sys = "workstation"
                path = install_path
          end if
          end if
       end if
    end if
    if ( sys .eq. "unknown" ) then
       print("+++ unknown system. Exiting...")
       exit
    end if


workpath = systemfunc("echo $PWD") + "/"


;******************************************************
;******************************************************
if man then
;specify name of output file:
outputfile      = "output.vdf" ;<<<<------------- !!!
end if
;******************************************************
;******************************************************

x_min = new(1,double)
x_max = new(1,double)
x_min = 99.0
x_max = -99.0
y_min = new(1,double)
y_max = new(1,double)
y_min = 99.0
y_max = -99.0
z_min = new(1,double)
z_max = new(1,double)
z_min = 99.0
z_max = -99.0

vars3d      = ""



if man then
    ;**** Get names of variables appearing in the VDF file, then get their dimensions with min/max x,y,z
    print(" ")
    print("Choose the *.nc files containing 3D data by typing their indices (seperated by ','). Start with the first files of the timeseries.")
    print(" ")
    files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
    print(" "+files_avail)


    files_in_ary  = toint(str_split(systemfunc("read file; echo $file"), ","))
    nofiles       = dimsizes(files_in_ary)
    file_in       = new(nofiles,string)

    do i=0,nofiles-1
        file_in(i)    = files_avail(files_in_ary(i))
    end do
else
    file_in = files
    nofiles = dimsizes(file_in)
end if





if (nofiles .ge. 2) then
nvars   = 1
if man then
varString  = new((/nofiles,nvars/),string)
end if; man
dimNames   = new((/nofiles,nvars,4/),string)
end if


do i=0,nofiles-1
    print("Loading file...")
    f             = addfile(workpath+file_in(i),"r")

    varNames      = getfilevarnames(f)


    if man then
    print(" ")
    print("The following variables were found in "+file_in(i)+":")
    print(" ")
    print(" "+varNames)
    print(" ")
    print("Type the index of the VARIABLES (not dimensions) you want to visualize in Vapor (seperated by ',')")

    vars       = toint(str_split(systemfunc("read vars; echo $vars"), ","))


    nvars = dimsizes(vars)
    end if; man

    if (nofiles .eq. 1) then
    if man then
    varString  = new((/nofiles,nvars/),string)
    end if; man
    dimNames   = new((/nofiles,nvars,4/),string)
    end if



    do n=0,nvars-1
       if man then
       varString(i,n) = varNames(vars(n))
       end if; man

       if alt_varnames then
            varString(i,n) = varNames(vars(i,n))
       end if

       if (vars3d .ne. "") then
       vars3d       = vars3d + ":" + varString(i,n)
       else
       vars3d       = varString(i,n)
       end if


       dims          = getfilevardimsizes(f, varString(i,n))
       dimNames(i,n,:) = getfilevardims(f, varString(i,n))



       time     = f->$dimNames(i,n,0)$
       t_max    = max(time)
       z        = f->$dimNames(i,n,1)$
       y        = f->$dimNames(i,n,2)$
       x            = f->$dimNames(i,n,3)$

       if (z_min .ge. min(z)) then
          z_min = min(z)
       end if

       if (z_max .le. max(z)) then
          z_max = max(z)
       end if

       if (y_min .ge. min(y)) then
          y_min = min(y)
       end if

       if (y_max .le. max(y)) then
          y_max = max(y)
       end if

       if (x_min .ge. min(x)) then
          x_min = min(x)
       end if

       if (x_max .le. max(x)) then
          x_max = max(x)
       end if

       delete(z)
       delete(y)
       delete(x)
    end do

end do

;print(varString)

if man then
delete(files_avail)
end if

if man then
    print(" ")
    periodic_resp              = systemfunc("read -p 'Periodic lateral boundaries? (y/n) ' RESP; echo $RESP")
    if (periodic_resp .eq. "y") then
        periodic = True
    else
        periodic = False
    end if
end if

if (periodic) then
   px           = 1
   py           = 1
else
   px           = 0
   py           = 0
end if
print(" ")
print(" ")

if man then
    grid_resp              = systemfunc("read -p 'Is grid stretching used? (y/n) YES will cut of heights above dz_stretch_level? ' RESP; echo $RESP")
    if (grid_resp .eq. "y") then
        grid_stretch = True
    else
        grid_stretch = False
    end if
end if

if (grid_stretch) then
    if man then
       stretch_level    = stringtointeger(systemfunc("read -p 'Please type height of dz_stretch_level (in meters) ' RESP; echo $RESP"))
    end if; man
    z_max = stretch_level ;**** limit z_max
end if

print(" ")

;*******************************************************************************




if (any(varNames .EQ. "zwwi")) then
   var2d     = " -vars2dxy HGT:topo"
   topo      = True
   print(" ")
   print("Found topography - will add 'zwwi' (w-grid) topography to *.vdf file.")
   print(" ")
else
   var2d     = " "
   topo      = False
end if

if man then
print(" ")

print(dims(0)+" time steps found in "+file_in)
ts       = systemfunc("read -p 'Include all time steps in *.vdf file? (y/n) ' RESP; echo $RESP")

if (ts .eq. "y") then
   t_start      = 0
   t_end        = dimsizes(time)-1
else
   t_start      = stringtointeger(systemfunc("read -p 'Type index of first time step: ' ts; echo $ts"))
   t_end        = stringtointeger(systemfunc("read -p 'Type index of last time step: ' ts; echo $ts"))
end if
end if

if (.not. man) then
    if (.not. spec_ts) then
           t_start      = 0
           t_end        = dimsizes(time)-1
    end if
end if

if man then
    print(" ")
    askreflevel              = systemfunc("read -p 'Specify a maximum refinement level? (y/n) No will set the default value.  ' RESP; echo $RESP")
    if (askreflevel .eq. "y") then
       reflevel    = systemfunc("read -p 'Please type maximum refinement level ' RESP; echo $RESP")
       reflevel_str = "-level "+ reflevel
    else
       reflevel_str = ""
    end if
    print(" ")
else
    if spec_reflevel then
        reflevel_str = "-level "+ reflevel
    else
       reflevel_str = ""
    end if
end if



;*******************************************************************************
;*******************************************************************************
;*******************************************************************************
print(" ")
print("Creating *.vdf file")
print(" ")
;vars3d    = str_join(varString,":")

a = systemfunc(path+"vdfcreate -dimension "+dims(3)+"x"+dims(2)+"x"+dims(1)+" -numts " + (t_end - t_start + 1) + var2d + " "+reflevel_str+" -vars3d "+vars3d+" -extents "+x_min+":"+y_min+":"+z_min+":"+x_max+":"+y_max+":"+z_max+" -periodic "+px+":"+py+":0 "+workpath+outputfile)
delete(a)
print("Done.")


;*******************************************************************************
print("Populating *.vdf file")

do i=0,nofiles-1
    do n=0,nvars-1
       dimnamesstring = str_join(dimNames(i,n,::-1),":")
       print(" ")
       print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
       print(" ")
       do t=t_start,t_end
          print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
          a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varString(i,n)+" -ncdfvar "+varString(i,n)+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(i))
          delete(a)
       end do
    end do
end do

if topo then
   print(" ")
   print("Adding topography (1)...")
   print(" ")
   do t=t_start,t_end
      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
      a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
      delete(a)
   end do
   print(" ")
   print("Adding topography (2)...")
   print(" ")
   do t=t_start,t_end
      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
      a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname topo -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
      delete(a)
   end do
end if


;*******************************************************************************
if man then
;**** Read more files and add them to VDF/VDC
addmore = "y"
do while (addmore .eq. "y")

   addmore       = systemfunc("read -p 'Add time steps or 2D cross sections from additional *.nc files? (2d/time/n) ' RESP; echo $RESP")

;**** Add time steps from additional 3D files
    if (addmore .eq. "time") then

    delete(f)

    do i=0,nofiles-1

        print(" ")
        print("Choose the *.nc file that contains 3D data for "+varString(i,0)+" by typing itÂŽs index.")
        print(" ")
        files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
        print(" "+files_avail)
        print(" ")
        file_in(i)    = files_avail(toint(systemfunc("read file; echo $file")))

    end do

    f        = addfile(workpath+file_in(0),"r")

    time     = f->time

    a = systemfunc(path+"vdfedit -addts "+dimsizes(time)+" "+workpath+outputfile)
    delete(a)

    do i=0,nofiles-1
        do n=0,nvars-1
           dimnamesstring = str_join(dimNames(i,n,::-1),":")
           print(" ")
           print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
           print(" ")
           do t=t_end,dimsizes(time)-1
              print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
              a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varString(i,n)+" -ncdfvar "+varString(i,n)+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(i))
              delete(a)
           end do
        end do
    end do

    if topo then
       print(" ")
       print("Adding topography (1)...")
       print(" ")
       do t=t_end,dimsizes(time)-1
          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
          a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
          delete(a)
       end do
       print(" ")
       print("Adding topography (2)...")
       print(" ")
       do t=t_end,dimsizes(time)-1
          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
          a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
          delete(a)
       end do
    end if



    t_end       = dimsizes(time)-1



    end if


;**** Add 2D cross sections
   if (addmore .eq. "2d") then
      delete(f)
      delete(file_in)

      print(" ")
      print("Choose additional file to add 2D cross sections. Type index.")
      print(" ")
      files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
      print(" "+files_avail)

      file_in       = files_avail(toint(systemfunc("read file; echo $file")))
      delete(files_avail)

      f             = addfile(workpath+file_in,"r")

      varNames2d      = getfilevarnames(f)

      crosstype       = systemfunc("read -p 'Please specify type of cross section (xy/xz/yz) ' RESP; echo $RESP")
      print(" ")
      print("The following variables were found in the file:")
      print(" ")
      print(" "+varNames2d)
      print(" ")
      print("Type the index of the variables you want to add to the *.vdf file (seperated by ',')")

      vars2d       = toint(str_split(systemfunc("read vars; echo $vars"), ","))



      nvars2d = dimsizes(vars2d)

      varString2d  = new(dimsizes(vars2d),string)
      dimNames2d   = new((/dimsizes(vars2d),4/),string)

      do n=0,nvars2d-1
         varString2d(n) = varNames2d(vars2d(n))

         dims2d     = getfilevardimsizes(f,varString2d(n))
         dimNames2d(n,:) = getfilevardims(f,varString2d(n))


      end do

      vars2dstr    = str_join(varString2d,":")

      a = systemfunc(path+"vdfedit -addvars2d"+crosstype+" "+vars2dstr+" "+workpath+outputfile)
      delete(a)

      do n=0,nvars2d-1
         dimnamesstring = str_join(dimNames2d(n,::-1),":")

         do t=0,dims2d(0)-1
            a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varNames2d(vars2d(n))+" -ncdfvar "+varNames2d(vars2d(n))+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in)
            delete(a)
         end do
      end do
      delete(varNames2d)
   end if

   addmore       = systemfunc("read -p 'Add additional *.nc files? (y/n) ' RESP; echo $RESP")
end do
end if; man
;*******************************************************************************
print("Done.")
;*************
end