load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" 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" ;*************************************************** ; load ncl_preferences.ncl ;*************************************************** function which_script() local script begin script="cross_section" return(script) end if (isfilepresent("~/ncl_preferences.ncl")) then loadscript("~/ncl_preferences.ncl") else if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl")) then loadscript( "~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl") else print(" ") print("'ncl_preferences.ncl' does not exist in $home or $home/palm/current_version/trunk/SCRIPTS/NCL/") print(" ") exit end if end if begin ; *************************************************** ; set default parameter values and strings if not assigned in prompt or parameter list ; *************************************************** if (file_1 .EQ. "File in") then print(" ") print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'") print(" ") exit else file_in = file_1 end if if (.not. isfilepresent(file_in)) then print(" ") print("1st input file: '"+file_in+"' does not exist") print(" ") exit end if if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm")then print(" ") print("'format_out = "+format_out+"' is invalid and set to'x11'") print(" ") format_out="x11" end if if (sort .NE. "layer" .AND. sort .NE. "time") then print(" ") print("'sort'= "+sort+" is invalid and set to 'layer'") print(" ") sort = "layer" end if if (mode .NE. "Fill" .AND. mode .NE. "Line" .AND. mode .NE. "Both") then print(" ") print("'mode'= "+mode+" is invalid and set to 'Fill'") print(" ") mode = "Fill" end if if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND. fill_mode .NE. "CellFill") then print(" ") print("Your 'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'") print(" ") fill_mode = "AreaFill" end if if (shape .NE. 0 .AND. shape .NE. 1) then print(" ") print("'shape'= "+shape+" is invalid and set to 1") print(" ") shape = 1 end if if (xyc .NE. 0 .AND. xyc .NE. 1)then print(" ") print("'xyc'= "+xyc+" is invalid and set to 0") print(" ") xyc = 0 end if if (xzc .NE. 0 .AND. xzc .NE. 1)then print(" ") print("'xzc'= "+xzc+" is invalid and set to 0") print(" ") xyc = 0 end if if (yzc .NE. 0 .AND. yzc .NE. 1)then print(" ") print("'yzc'= "+yzc+" is invalid and set to 0") print(" ") xyc = 0 end if if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then print(" ") print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)") print(" ") exit end if if (xyc .EQ. 1 ) then if (xzc .EQ. 1 .OR. yzc .EQ. 1) then print(" ") print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)") print(" ") exit end if end if if (xzc .EQ. 1) then if (xyc .EQ. 1 .OR. yzc .EQ. 1) then print(" ") print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)") print(" ") exit end if end if if (yzc .EQ. 1) then if (xyc .EQ. 1 .OR. xzc .EQ. 1) then print(" ") print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)") print(" ") exit end if end if if (vector .NE. 0 .AND. vector .NE. 1) then print(" ") print("Please set 'vector' to 0 or 1") print(" ") exit end if ; *************************************************** ; open input file ; *************************************************** f = addfile( file_in, "r" ) vNam = getfilevarnames(f) print(" ") print("Variables in input file:") print("- "+ vNam) print(" ") dim = dimsizes(vNam) ; *************************************************** ; check for kind of input file ; *************************************************** check_3d = True do varn=0,dim-1 checkxy = isStrSubset( vNam(varn),"xy") checkxz = isStrSubset( vNam(varn),"xz") checkyz = isStrSubset( vNam(varn),"yz") if (checkxy .OR. checkxz .OR. checkyz) then check_3d = False break end if end do if (.not. check_3d) then if (xyc .EQ. 1)then if (.not. checkxy) then print(" ") print("Your input file doesn't have values for xy cross-sections") if (checkxz)then print("Select another input file or xz cross-sections") else print("Select another input file or yz cross-sections") end if print(" ") exit else print(" ") print("Your input file contains xy data") print(" ") end if end if if (xzc .EQ. 1)then if (.not. checkxz) then print(" ") print("Your input file doesn't have values for xz cross-sections") if (checkxy)then print("Select another input file or xy cross-sections") else print("Select another input file or yz cross-sections") end if print(" ") exit else print(" ") print("Your input file contains xz data") print(" ") end if end if if (yzc .EQ. 1)then if (.not. checkyz) then print(" ") print("Your input file doesn't have values for yz cross-sections") if (checkxy)then print("Select another input file or xy cross-sections") else print("Select another input file or xz cross-sections") end if print(" ") exit else print(" ") print("Your input file contains yz data") print(" ") end if end if else print(" ") print("Your input file: '"+file_in+"'") print("contains 3d or other data") print(" ") end if if (dim .EQ. 0) then print(" ") print("There is no data on file") print(" ") exit end if ; *************************************************** ; set up recourses ; *************************************************** cs_res = True cs_res@gsnYAxisIrregular2Linear = True cs_res@gsnDraw = False cs_res@gsnFrame = False cs_res@gsnMaximize = True cs_res@gsnPaperOrientation = "portrait" cs_res@gsnPaperWidth = 8.27 cs_res@gsnPaperHeight = 11.69 cs_res@gsnPaperMargin = 0.79 cs_res@tmXBLabelFontHeightF = .03 cs_res@tmYLLabelFontHeightF = .03 cs_res@tiXAxisFontHeightF = .03 cs_res@tiYAxisFontHeightF = .03 cs_res@tmXBMode ="Automatic" cs_res@tmYLMode ="Automatic" cs_res@lgTitleFontHeightF = .03 cs_res@lgLabelFontHeightF = .03 cs_res@txFontHeightF = .03 cs_res@cnLevelSelectionMode = "ManualLevels" cs_resP = True cs_resP@txString = f@title cs_resP@txFuncCode = "~" cs_resP@txFontHeightF = .02 if ( mode .eq. "Fill" ) then cs_res@cnFillOn = True cs_res@gsnSpreadColors = True cs_res@cnFillMode = fill_mode cs_res@lbOrientation = "Vertical" cs_res@cnLinesOn = False cs_res@cnLineLabelsOn = False end if if ( mode .eq. "Both" ) then cs_res@cnFillOn = True cs_res@gsnSpreadColors = True cs_res@cnFillMode = fill_mode cs_res@lbOrientation = "Vertical" cs_res@cnLinesOn = True cs_res@cnLineLabelsOn = True end if ; ********************************************* ; set up of start_time_step and end_time_step ; ********************************************* t_all = f->time nt = dimsizes(t_all) delta_t = t_all(nt-1)/nt ; **************************************************** ; start of time step and different types of mistakes that could be done ; **************************************************** if (start_time_step .EQ. -1.) then start_time_step=t_all(0)/3600 else if (start_time_step .GT. t_all(nt-1)/3600) print(" ") print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h") print(" ") print("Please select another 'start_time_step'") print(" ") exit end if if (start_time_step .LT. t_all(0)/3600) print(" ") print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h") print(" ") print("Please select another 'start_time_step'") print(" ") exit end if end if do i=0,nt-1 if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. start_time_step .LT. (t_all(i)+delta_t/2)/3600)then st=i break end if end do if (.not. isvar("st"))then print(" ") print("'start_time_step' = "+ start_time_step +"h is invalid") print(" ") print("Please select another 'start_time_step'") print(" ") exit end if ; **************************************************** ; end of time step and different types of mistakes that could be done ; **************************************************** if (end_time_step .EQ. -1.) then end_time_step = t_all(nt-1)/3600 else if (end_time_step .LT. t_all(0)/3600) print(" ") print("'end_time_step' = "+end_time_step+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h") print(" ") print("Please select another 'end_time_step'") print(" ") exit end if if (end_time_step .GT. t_all(nt-1)/3600) print(" ") print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h") print(" ") print("Please select another 'end_time_step'") print(" ") exit end if if (end_time_step .LT. start_time_step/3600) print(" ") print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h") print(" ") print("Please select another 'start_time_step' or 'end_time_step'") print(" ") exit end if end if do i=0,nt-1 if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. end_time_step .LT. (t_all(i)+delta_t/2)/3600)then et=i break end if end do if (.not. isvar("et"))then print(" ") print("'end_time_step' = "+ end_time_step +"h is invalid") print(" ") print("Please select another 'end_time_step'") print(" ") exit end if delete(start_time_step) start_time_step=round(st,3) delete(end_time_step) end_time_step=round(et,3) print(" ") print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step) print(" till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step) print(" ") no_time=(end_time_step-start_time_step)+1 ; **************************************************** ; set up variables for vector plot if required ; **************************************************** if (vector .EQ. 1) then if (vec1 .EQ. "component1") then print(" ") print("Please indicate Vector 1 ('vec1') for Vector-Plot or set 'vector' to 0") print(" ") exit end if if (vec2 .EQ. "component2") then print(" ") print("Please indicate Vector 2 ('vec2') for Vector-Plot or set 'vector' to 0") print(" ") exit end if end if check_vec1 = False check_vec2 = False check_vecp = False ; **************************************************** ; get data for plots ; **************************************************** if (xyc .EQ. 1) then do varn=0,dim-1 if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then x_d = f->$vNam(varn)$ xdim = dimsizes(x_d)-1 delta_x = x_d(1)-x_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then y_d = f->$vNam(varn)$ ydim = dimsizes(y_d)-1 delta_y = y_d(1)-y_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then z_d = f->$vNam(varn)$ zdim = dimsizes(z_d)-1 delta_z = 0 break else if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then z_d = f->$vNam(varn)$ zdim = dimsizes(z_d)-1 delta_z = -1.d break else zdim = 0 delta_z = -1.d end if end if end do end if if (xzc .EQ. 1) then do varn=0,dim-1 if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x") then x_d = f->$vNam(varn)$ xdim = dimsizes(x_d)-1 delta_x = x_d(1)-x_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then z_d = f->$vNam(varn)$ zdim = dimsizes(z_d)-1 delta_z = z_d(1)-z_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then y_d = f->$vNam(varn)$ ydim = dimsizes(y_d)-1 delta_y = y_d(1)-y_d(0) break else if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then y_d = f->$vNam(varn)$ ydim = dimsizes(y_d)-1 delta_y = -1.d break else ydim = 0 delta_y = -1.d end if end if end do end if if (yzc .EQ. 1) then do varn=0,dim-1 if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then y_d = f->$vNam(varn)$ ydim = dimsizes(y_d)-1 delta_y = y_d(1)-y_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then z_d = f->$vNam(varn)$ zdim = dimsizes(z_d)-1 delta_z = z_d(1)-z_d(0) break end if end do do varn=0,dim-1 if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x") x_d = f->$vNam(varn)$ xdim = dimsizes(x_d)-1 delta_x = x_d(1)-x_d(0) break else if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then x_d = f->$vNam(varn)$ xdim = dimsizes(x_d)-1 delta_x = -1.d break else xdim = 0 delta_x = -1.d end if end if end do end if ; **************************************************** ; set up ranges of x-, y- and z-coordinates ; **************************************************** if (xs .EQ. -1.d) then xs = x_d(0) if (delta_x .EQ. -1) then print(" ") print("You cannot choose a start value for x, there are preseted layers for x") print(" ") xstart=0 end if else if (delta_x .EQ. -1) then print(" ") print("You cannot choose a start value for x, there are preseted layers for x") print(" ") xstart=0 else if (xs .LT. 0-delta_x/2) then print(" ") print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m") print(" ") exit end if if (xyc .EQ. 1 .OR. xzc .EQ. 1) then if (xs .GE. xdim*delta_x) then print(" ") print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m") print(" ") exit end if else if (xs .GT. xdim*delta_x) then print(" ") print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m") print(" ") exit end if end if end if end if if (ys .EQ. -1.d) then ys = y_d(0) if (delta_y .EQ. -1) then print(" ") print("You cannot choose a start value for y, there are preseted layers for y") print(" ") ystart=0 end if else if (delta_y .EQ. -1) then print(" ") print("You cannot choose a start value for y, there are preseted layers for y") print(" ") ystart=0 else if (ys .LT. 0-delta_y/2) then print(" ") print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m") print(" ") exit end if if (xyc .EQ. 1 .OR. yzc .EQ. 1) then if (ys .GE. ydim*delta_y) then print(" ") print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m") print(" ") exit end if else if (ys .GT. ydim*delta_y) then print(" ") print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m") print(" ") exit end if end if end if end if if (zs .EQ. -1) then zs = 0 if (delta_z .EQ. -1) then print(" ") print("You cannot choose a start value for z, there are preseted layers for z") print(" ") end if else if (delta_z .EQ. -1) then print(" ") print("You cannot choose a start value for z, there are preseted layers for z") print(" ") else if (zs .LT. 0) then print(" ") print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0") print(" ") exit end if if (xzc .EQ. 1 .OR. yzc .EQ. 1) then if (zs .GE. zdim) then print(" ") print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim) print(" ") exit end if else if (zs .GT. zdim) then print(" ") print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim) print(" ") exit end if end if end if end if if (xe .EQ. -1) then xe = x_d(xdim) if (delta_x .EQ. -1) then print(" ") print("You cannot choose an end value for x, there are preseted layers for x") print(" ") xend=xdim end if else if (delta_x .EQ. -1) then print(" ") print("You cannot choose an end value for x, there are preseted layers for x") print(" ") xend=xdim else if (xe .GT. xdim*delta_x) then print(" ") print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m") print(" ") exit end if if (xyc .EQ. 1 .OR. xzc .EQ. 1) then if (xe .LE. 0-delta_x/2) print(" ") print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m") print(" ") exit end if if (xe .LE. xs) then print(" ") print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m") print(" ") exit end if if (xe .EQ. xs+1) then print(" ") print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m") print(" ") exit end if else if (xe .LT. 0-delta_x/2) print(" ") print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m") print(" ") exit end if if (xe .LT. xs) then print(" ") print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m") print(" ") exit end if end if end if end if if (ye .EQ. -1) then ye = y_d(ydim) if (delta_y .EQ. -1) then print(" ") print("You cannot choose an end value for y, there are preseted layers for y") print(" ") yend=ydim end if else if (delta_y .EQ. -1) then print(" ") print("You cannot choose an end value for y, there are preseted layers for y") print(" ") yend=ydim else if (ye .GT. ydim*delta_y) then print(" ") print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m") print(" ") exit end if if (xyc .EQ. 1 .OR. yzc .EQ. 1) then if (ye .LE. 0-delta_y/2) print(" ") print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m") print(" ") exit end if if (ye .LE. ys) then print(" ") print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m") print(" ") exit end if if (ye .EQ. ys+1) then print(" ") print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m") print(" ") exit end if else if (ye .LT. 0-delta_y/2) print(" ") print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m") print(" ") exit end if if (ye .LT. ys) then print(" ") print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m") print(" ") exit end if end if end if end if if (ze .EQ. -1) then ze = zdim if (delta_z .EQ. -1) then print(" ") print("You cannot choose an end value for z, there are preseted layers for z") print(" ") ze = zdim end if else if (delta_z .EQ. -1) then print(" ") print("You cannot choose an end value for z, there are preseted layers for z") print(" ") ze = zdim else if (ze .GT. zdim) then print(" ") print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim) print(" ") exit end if if (xzc .EQ. 1 .OR. yzc .EQ. 1) then if (ze .LE. 0) print(" ") print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0") print(" ") exit end if if (ze .LE. zs) then print(" ") print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs) print(" ") exit end if if (ze .EQ. zs+1) then print(" ") print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs) print(" ") exit end if else if (ze .LT. 0) print(" ") print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0") print(" ") exit end if if (ze .LT. zs) then print(" ") print("range end for z-coordinate = "+ze+" is lower than start range = "+zs) print(" ") exit end if end if end if end if if (delta_x .NE. -1) then do i=0,xdim if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then xstart=i break end if end do do i=0,xdim if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then xend=i break end if end do end if if (delta_y .NE. -1) then do i=0,ydim if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then ystart=i break end if end do do i=0,ydim if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then yend=i break end if end do end if if( shape .eq. 1 ) then if (xyc .EQ. 1)then cs_res@vpWidthF = (xe-xs)/(ye-ys) cs_res@vpHeightF = 1 end if if (xzc .EQ. 1)then cs_res@vpWidthF = (xe-xs)/(delta_x*(ze-zs)) cs_res@vpHeightF = 1 end if if (yzc .EQ. 1)then cs_res@vpWidthF = (ye-ys)/(delta_y*(ze-zs)) cs_res@vpHeightF = 1 end if end if delete(xs) delete(xe) delete(ys) delete(ye) xs=xstart xe=xend ys=ystart ye=yend if (xyc .EQ. 1 .OR. xzc .EQ.1)then if (xe .EQ. xs+1) then print(" ") print("range end for x-coordinate="+xe*delta_x+"m) must be at least two") print("more gridpoints(="+2*delta_x+"m) greater than start range="+xs*delta_x+"m)") print(" ") exit end if end if if (xyc .EQ. 1 .OR. yzc .EQ.1)then if (ye .EQ. ys+1) then print(" ") print("range end for y-coordinate="+ye*delta_y+"m) must be at least two") print("more gridpoints(="+2*delta_y+"m) greater than start range="+ys*delta_y+"m)") print(" ") exit end if end if if (xzc .EQ. 1 .OR. yzc .EQ.1)then if (ze .EQ. zs+1) then print(" ") print("range end for x-coordinate="+ze+") must be at least two") print("more gridpoints greater than start range="+zs+")") print(" ") exit end if end if if (xyc .EQ. 1) then no_layer = (ze-zs)+1 data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float) end if if (xzc .EQ. 1) then no_layer = (ye-ys)+1 data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float) end if if (yzc .EQ. 1) then no_layer = (xe-xs)+1 data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float) end if MinVal = new(dim,float) MaxVal = new(dim,float) ; **************************************************** ; define inner and outer loops depending on "sort" ; **************************************************** if ( xyc .eq. 1 ) then lays = zs laye = ze end if if ( xzc .eq. 1 ) then lays = ys laye = ye end if if ( yzc .eq. 1) then lays = xs laye = xe end if if ( sort .eq. "time" ) then lis = start_time_step lie = end_time_step los = 0 loe = laye-lays else lis = 0 lie = laye-lays los = start_time_step loe = end_time_step end if check = True v1=0 v2=0 no_var=0 n=0 do varn=dim-1,0,1 if (vector .EQ. 1) then check_vec1 = isStrSubset( vec1,","+vNam(varn)+",") check_vec2 = isStrSubset( vec2,","+vNam(varn)+",") end if if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then check = False else check = True end if if(check) then no_var=no_var+1 if (vector .EQ. 1) then if (","+vNam(varn)+"," .EQ. vec1) then v1=v1+1 end if if (","+vNam(varn)+"," .EQ. vec2) then v2=v2+1 end if end if if (xyc .EQ. 1) then data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) end if if ( xzc .eq. 1 ) then data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) end if if ( yzc .eq. 1) then data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) end if if (check_vec1) then vect1=data(varn,:,:,:,:) end if if (check_vec2) then vect2=data(varn,:,:,:,:) end if data!0 = "var" data!1 = "t" data!2 = "z" data!3 = "y" data!4 = "x" MinVal(varn) = min(data(varn,:,:,:,:)) MaxVal(varn) = max(data(varn,:,:,:,:)) end if end do var_input=new(no_var,string) numb=0 no_var=0 do varn=dim-1,0,1 if (vector .EQ. 1) then check_vec1 = isStrSubset( vec1,","+vNam(varn)+",") check_vec2 = isStrSubset( vec2,","+vNam(varn)+",") end if if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then check = False else check = True end if if(check) then var_input(numb)=vNam(varn) numb=numb+1 end if if (var .NE. "all") then check = isStrSubset( var,","+vNam(varn)+"," ) end if if (check) then no_var = no_var+1 end if end do if (no_var .EQ. 0) then print(" ") print("The variables var='"+var+"' do not exist on your input file;") print("be sure to have one comma berfore and after each variable") print(" ") exit end if if (vector .EQ. 1) then if (v1 .EQ. 0)then print(" ") print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:") print("- "+var_input) print("be sure to have one comma berfore and after the variable") print(" ") exit end if if (v2 .EQ. 0)then print(" ") print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:") print("- "+var_input) print("be sure to have one comma berfore and after the variable") print(" ") exit end if end if ; *************************************************** ; open workstation(s) ; *************************************************** wks_ps = gsn_open_wks(format_out,file_out) gsn_define_colormap(wks_ps,"rainbow+white") plot=new((/no_time*no_layer*no_var/),graphic) page = 0 n=0 print(" ") print("Plots sorted by '"+sort+"'") print(" ") ; *************************************************** ; create plots ; *************************************************** if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then do lo = los, loe do li = lis, lie level = "=" + data&z(lo) + "m" vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2 vecres@gsnLeftString = "t=" + data&t(li) +"s z"+level vecres@tiXAxisString = " " plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres) n=n+1 end do end do end if do varn=dim-1,0,1 if (vector .EQ. 1) then check_vecp = isStrSubset( plotvec,","+vNam(varn)+",") end if if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then check = False else check = True end if if (var .NE. "all") then check = isStrSubset( var,","+vNam(varn)+"," ) end if if(check) then space=(MaxVal(varn)-MinVal(varn))/24 cs_res@cnMinLevelValF = MinVal(varn) cs_res@cnMaxLevelValF = MaxVal(varn) cs_res@cnLevelSpacingF = space ; **************************************************** ; loops over time and layer ; **************************************************** do lo = los, loe do li = lis, lie ; **************************************************** ; xy cross section ; **************************************************** if ( xyc .eq. 1 ) then cs_res@tiXAxisString = "x [m]" cs_res@tiYAxisString = "y [m]" cs_res@gsnLeftString = "Plot of "+vNam(varn) if ( sort .eq. "time" ) then if ( data&z(lo) .eq. -1)then if (delta_z .EQ. -1) then level = "-average" else level = "=" + data&z(lo) + "m" end if else level = "=" + data&z(lo) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h z"+level plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res) if (vector .EQ. 1 .AND. check_vecp) then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " vecres@gsnLeftString = " " vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if (data&z(li) .eq. -1) then if (delta_z .EQ. -1) then level = "-average" else level = "=" + data&z(li) + "m" end if else level = "=" + data&z(li) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h z"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res) if (vector .EQ. 1 .AND. check_vecp) then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres) overlay(plot(n), plot_vec) end if end if end if ; **************************************************** ; xz cross section ; **************************************************** if ( xzc .eq. 1 ) then cs_res@tiXAxisString = "x [m]" cs_res@tiYAxisString = "z [m]" cs_res@gsnLeftString = "Plot of "+vNam(varn) if ( sort .eq. "time" ) then if ( data&y(lo) .eq. -1 ) then level = "-average" else level = "=" + data&y(lo) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "s y"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res) if (vector .EQ. 1 .AND. check_vecp) then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if ( data&y(li) .eq. -1 ) then level = "-average" else level = "=" + data&y(li) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "s y"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res) if (vector .EQ. 1 .AND. check_vecp) then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres) overlay(plot(n), plot_vec) end if end if end if ; **************************************************** ; yz cross section ; **************************************************** if ( yzc .eq. 1 ) then cs_res@tiXAxisString = "y [m]" cs_res@tiYAxisString = "z [m]" cs_res@gsnLeftString = "Plot of "+vNam(varn) if ( sort .eq. "time" ) then if ( data&x(lo) .eq. -1 ) then level = "-average" else level = "=" + data&x(lo) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "s x"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res) if (vector .EQ. 1 .AND. check_vecp) then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if ( data&x(li) .eq. -1 ) then level = "-average" else level = "=" + data&x(li) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "s x"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res) if (vector .EQ. 1 .AND. check_vecp)then vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag vecres@vcRefLengthF = 0.05 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres) overlay(plot(n), plot_vec) end if end if end if n=n+1 end do end do end if end do ; *************************************************** ; merge plots onto one page ; *************************************************** if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP) print(" ") print("Outputs to .eps or .epsi have only one frame") print(" ") else do np = 0,no_layer*no_time*(no_var+1)-1,no_lines*no_columns if ( np + no_lines*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_lines,no_columns/),cs_resP) else gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP) end if end do end if else if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP) print(" ") print("Outputs to .eps or .epsi have only one frame") print(" ") else do np = 0,no_layer*no_time*no_var-1,no_lines*no_columns if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_lines,no_columns/),cs_resP) else gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP) end if end do end if end if print(" ") print("Output to: " + file_out +"."+ format_out) print(" ") end