;--------------------------------------------------------------------------------! ; 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 . ; ; Copyright 1997-2012 Leibniz University Hannover ;------------------------------------------------------------------------------! ; ; Current revisions: ; ----------------- ; adapted for use on desktop machines ; ; Former revisions: ; ----------------- ; $Id: nc2vdf.ncl 1110 2013-03-07 11:42:45Z raasch $ ; ; 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