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.config or .ncl.config.default ;*************************************************** function which_script() local script begin script="cross_section" return(script) end if (isfilepresent("$PALM_BIN/../../.ncl.config")) then loadscript("$PALM_BIN/../../.ncl.config") else if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then loadscript( "$PALM_BIN/NCL/.ncl.config.default") else palm_bin_path = getenv("PALM_BIN") print(" ") print("Neither the personal configuration file '.ncl.config' exists in") print("~/palm/current_version") print("nor the default configuration file '.ncl.config.default' exists in") print(palm_bin_path + "/NCL") print(" ") exit end if end if begin ;*************************************************** ; Retrieving the NCL version used ;*************************************************** ncl_version_ch = systemfunc("ncl -V") ncl_version = stringtofloat(ncl_version_ch) ; *************************************************** ; Retrieving the double quote character ; *************************************************** dq=str_get_dq() ; *************************************************** ; set default parameter values and strings if not assigned in prompt or parameter list ; *************************************************** if (file_1 .EQ. "File in") then print(" ") print("Declare input file 'file_1=' in '.ncl.config' or prompt") print(" ") exit else file_in = file_1 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" .AND. format_out .NE. "png")then print(" ") print("'format_out = "+format_out+"' is invalid and set to'x11'") print(" ") format_out="x11" end if if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then print(" ") print("Output of png files not available") print("png output is avaiable with NCL version 5.2.0 and higher ") print("NCL version used: " + ncl_version_ch) print(" ") exit 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("'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("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("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("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("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("Set 'vector' to 0 or 1") print(" ") exit end if if (norm_x .EQ. 0) then print(" ") print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0") print(" ") norm_x = 1.0 end if if (norm_y .EQ. 0) then print(" ") print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0") print(" ") norm_y = 1.0 end if if (norm_z .EQ. 0) then print(" ") print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0") print(" ") norm_z = 1.0 end if ; *************************************************** ; open input file ; *************************************************** file_in_1 = False if (isStrSubset(file_in, ".nc"))then start_f = -2 end_f = -2 file_in_1 = True end if if (start_f .EQ. -1)then print(" ") print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)") print(" ") exit end if if (end_f .EQ. -1)then print(" ") print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)") print(" ") exit end if files=new(end_f-start_f+1,string) if (file_in_1)then if (isfilepresent(file_in))then files(0)=file_in else print(" ") print("1st input file: '"+file_in+"' does not exist") print(" ") exit end if else if (start_f .EQ. 0)then if (isfilepresent(file_in+".nc"))then files(0)=file_in+".nc" do i=1,end_f if (isfilepresent(file_in+"."+i+".nc"))then files(i)=file_in+"."+i+".nc" else print(" ") print("Input file: '"+file_in+"."+i+".nc' does not exist") print(" ") exit end if end do else print(" ") print("Input file: '"+file_in+".nc' does not exist") print(" ") exit end if else do i=start_f,end_f if (isfilepresent(file_in+"."+i+".nc"))then files(i-start_f)=file_in+"."+i+".nc" else print(" ") print("Input file: '"+file_in+"."+i+".nc' does not exist") print(" ") exit end if end do end if end if f=addfiles(files,"r") f_att=addfile(files(0),"r") ListSetType(f,"cat") vNam = getfilevarnames(f_att) 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 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 vecres = True cs_res@gsnYAxisIrregular2Linear = True cs_res@gsnDraw = False cs_res@gsnFrame = False cs_res@gsnMaximize = True cs_res@tmXBLabelFontHeightF = font_size cs_res@tmYLLabelFontHeightF = font_size cs_res@tiXAxisFontHeightF = font_size cs_res@tiYAxisFontHeightF = font_size cs_res@tmXBMode ="Automatic" cs_res@tmYLMode ="Automatic" cs_res@cnLevelSelectionMode = "ManualLevels" cs_res@lbLabelFontHeightF = font_size_legend cs_res@lbLabelStride = legend_label_stride cs_resP = True cs_resP@txString = f_att@title cs_resP@txFuncCode = "~" cs_resP@txFontHeightF = 0.015 if ( mode .eq. "Fill" ) then cs_res@cnFillOn = True cs_res@gsnSpreadColors = True cs_res@cnFillMode = fill_mode 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@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("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("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("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("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("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("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("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("Indicate Vector 1 ('vec1') for Vector-Plot or set 'vector' to 0") print(" ") exit end if if (vec2 .EQ. "component2") then print(" ") print("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_att->$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_att->$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_att->$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_att->$vNam(varn)$ zdim = dimsizes(z_d)-1 delta_z = -1.d break else zdim = 0 delta_z = -1.d end if end if if (vNam(varn) .eq. "zu1_xy") then zu1 = f_att->$vNam(varn)$ 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_att->$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_att->$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_att->$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_att->$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_att->$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_att->$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_att->$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_att->$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(" ") zs = 0 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 vecres@vpWidthF = (xe-xs)/(ye-ys) vecres@vpHeightF = 1 if (xe-xs .GT. ye-ys)then if (no_rows .LE. no_columns)then cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" else cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" end if else if (no_rows .GE. no_columns)then cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" else cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" end if end if end if if (xzc .EQ. 1)then cs_res@vpWidthF = (xe-xs)/(delta_x*(ze-zs)) cs_res@vpHeightF = 1 vecres@vpWidthF = (xe-xs)/(delta_x*(ze-zs)) vecres@vpHeightF = 1 if (xe-xs .GT. (delta_x*(ze-zs)))then if (no_rows .LE. no_columns)then cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" else cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" end if else if (no_rows .GE. no_columns)then cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" else cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" end if end if end if if (yzc .EQ. 1)then cs_res@vpWidthF = (ye-ys)/(delta_y*(ze-zs)) cs_res@vpHeightF = 1 vecres@vpWidthF = (ye-ys)/(delta_y*(ze-zs)) vecres@vpHeightF = 1 if (ye-ys .GT. (delta_y*(ze-zs)))then if (no_rows .LE. no_columns)then cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" else cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" end if else if (no_rows .GE. no_columns)then cs_res@gsnPaperOrientation = "portrait" vecres@gsnPaperOrientation = "portrait" cs_res@lbOrientation = "Vertical" else cs_res@gsnPaperOrientation = "landscape" vecres@gsnPaperOrientation = "landscape" cs_res@lbOrientation = "Horizontal" end if end if end if end if delete(xs) delete(xe) delete(ys) delete(ye) xs=xstart xe=xend ys=ystart ye=yend if (xyc .EQ. 1) then d= (ye-ys+1)/(major_ticks_y-1) e=(xe-xs+1)/(major_ticks_x-1) array_yl =new(major_ticks_y,integer) array_yl_labels=new(major_ticks_y,double) array_minor_yl =new((major_ticks_y-1)*4,double) array_xb =new(major_ticks_x,integer) array_xb_labels=new(major_ticks_x,double) array_minor_xb =new((major_ticks_x-1)*4,double) array_yl(0)=ys array_xb(0)=xs array_yl_labels(0)=y_d(ys) array_xb_labels(0)=x_d(xs) do ar=1,major_ticks_y-1 array_yl(ar)=d*(ar-1)+d-1 array_yl_labels(ar) = y_d(array_yl(ar)) if (ar .GT. 0) do min_ar=0,3 array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+(y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1) end do end if end do do br=1,major_ticks_x-1 array_xb(br)=e*(br-1)+e-1 array_xb_labels(br) = x_d(array_xb(br)) if (br .GT. 0) do min_br=0,3 array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1) end do end if end do end if if (xzc .EQ. 1) then d=(ze-zs+1)/(major_ticks_y-1) e=(xe-xs+1)/(major_ticks_x-1) array_yl =new(major_ticks_y,integer) array_yl_labels=new(major_ticks_y,double) array_minor_yl =new((major_ticks_y-1)*4,double) array_xb =new(major_ticks_x,integer) array_xb_labels=new(major_ticks_x,double) array_minor_xb =new((major_ticks_x-1)*4,double) array_yl(0)=zs array_xb(0)=xs array_yl_labels(0)=z_d(zs) array_xb_labels(0)=x_d(xs) do ar=1,major_ticks_y-1 array_yl(ar)=d*(ar-1)+d-1 array_yl_labels(ar) = z_d(array_yl(ar)) if (ar .GT. 0) do min_ar=0,3 array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1) end do end if end do do br=1,major_ticks_x-1 array_xb(br)=e*(br-1)+e-1 array_xb_labels(br) = x_d(array_xb(br)) if (br .GT. 0) do min_br=0,3 array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1) end do end if end do end if if (yzc .EQ. 1) then d=(ze-zs+1)/(major_ticks_y-1) e=(ye-ys+1)/(major_ticks_x-1) array_yl =new(major_ticks_y,integer) array_yl_labels=new(major_ticks_y,double) array_minor_yl =new((major_ticks_y-1)*4,double) array_xb =new(major_ticks_x,integer) array_xb_labels=new(major_ticks_x,double) array_minor_xb =new((major_ticks_x-1)*4,double) array_yl(0)=zs array_xb(0)=ys array_yl_labels(0)=z_d(zs) array_xb_labels(0)=y_d(ys) do ar=1,major_ticks_y-1 array_yl(ar)=d*(ar-1)+d-1 array_yl_labels(ar) = z_d(array_yl(ar)) if (ar .GT. 0) do min_ar=0,3 array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1) end do end if end do do br=1,major_ticks_x-1 array_xb(br)=e*(br-1)+e-1 array_xb_labels(br) = y_d(array_xb(br)) if (br .GT. 0) do min_br=0,3 array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+(y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1) end do end if end do end if if (axes_explicit .EQ. 1)then cs_res@tmYLMode = "Explicit" cs_res@tmXBMode = "Explicit" if (xyc .EQ. 1)then cs_res@tmYLValues = y_d(array_yl) cs_res@tmXBValues = x_d(array_xb) cs_res@tmYLLabels = array_yl_labels/norm_y cs_res@tmXBLabels = array_xb_labels/norm_x if (norm_x .NE. 1.)then cs_res@tiXAxisString = "x ["+unit_x+"]" else cs_res@tiXAxisString = "x [m]" end if if (norm_y .NE. 1.)then cs_res@tiYAxisString = "y ["+unit_y+"]" else cs_res@tiYAxisString = "y [m]" end if end if if (xzc .EQ. 1)then cs_res@tmYLValues = z_d(array_yl) cs_res@tmXBValues = x_d(array_xb) cs_res@tmYLLabels = array_yl_labels/norm_z cs_res@tmXBLabels = array_xb_labels/norm_x if (norm_x .NE. 1.)then cs_res@tiXAxisString = "x ["+unit_x+"]" else cs_res@tiXAxisString = "x [m]" end if if (norm_z .NE. 1.)then cs_res@tiYAxisString = "z ["+unit_z+"]" else cs_res@tiYAxisString = "z [m]" end if end if if (yzc .EQ. 1)then cs_res@tmYLValues = z_d(array_yl) cs_res@tmXBValues = y_d(array_xb) cs_res@tmYLLabels = array_yl_labels/norm_z cs_res@tmXBLabels = array_xb_labels/norm_y if (norm_y .NE. 1.)then cs_res@tiXAxisString = "y ["+unit_y+"]" else cs_res@tiXAxisString = "y [m]" end if if (norm_z .NE. 1.)then cs_res@tiYAxisString = "z ["+unit_z+"]" else cs_res@tiYAxisString = "z [m]" end if end if cs_res@tmYLMinorValues = array_minor_yl cs_res@tmXBMinorValues = array_minor_xb else if (axes_explicit .EQ. 0)then cs_res@tmYLMinorPerMajor = 4 cs_res@tmXBMinorPerMajor = 4 else print(" ") print("'axes_explicit' is invalid and set to 0") print(" ") axes_explicit = 0 cs_res@tmYLMinorPerMajor = 4 cs_res@tmXBMinorPerMajor = 4 end if if (xyc .EQ. 1)then cs_res@tiXAxisString = "x [m]" cs_res@tiYAxisString = "y [m]" end if if (xzc .EQ. 1)then cs_res@tiXAxisString = "x [m]" cs_res@tiYAxisString = "z [m]" end if if (yzc .EQ. 1)then cs_res@tiXAxisString = "y [m]" cs_res@tiYAxisString = "z [m]" end if end if 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) unit = new(dim,string) ; **************************************************** ; 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 = lays loe = laye else lis = lays lie = laye los = start_time_step loe = end_time_step end if check = True v1=0 v2=0 no_var=0 n=0 no_zu1=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 (var .NE. "all") then check = isStrSubset( var,","+vNam(varn)+"," ) 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 temp = f[:]->$vNam(varn)$ data_att = f_att->$vNam(varn)$ if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \ .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \ .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \ .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy" \ .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \ .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\ .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy" \ .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then ;these variables depend on zu1_xy and that's why they have only one z-layer data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe) no_zu1=no_zu1+1 else data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) end if delete(temp) end if if ( xzc .eq. 1 ) then temp = f[:]->$vNam(varn)$ data_att = f_att->$vNam(varn)$ data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) delete(temp) end if if ( yzc .eq. 1) then temp = f[:]->$vNam(varn)$ data_att = f_att->$vNam(varn)$ data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) delete(temp) 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,start_time_step:end_time_step,0:(ze-zs),0:(ye-ys),0:(xe-xs))) MaxVal(varn) = max(data(varn,start_time_step:end_time_step,0:(ze-zs),0:(ye-ys),0:(xe-xs))) unit(varn) = data_att@units delete(data_att) 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 before and after each variable") print(" ") exit end if 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 (var .NE. "all") then check = isStrSubset( var,","+vNam(varn)+"," ) end if if(check) then var_input(numb)=vNam(varn) numb=numb+1 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 before 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 before 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 before and after the variable") print(" ") exit end if end if ; *************************************************** ; open workstation(s) ; *************************************************** if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then format_out@wkPaperSize = "A4" end if if ( format_out .EQ. "png" ) then format_out@wkWidth = 1000 format_out@wkHeight = 1000 end if wks_ps = gsn_open_wks(format_out,file_out) gsn_define_colormap(wks_ps,"rainbow+white") if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then plot=new((/no_time*no_layer*no_var - no_time*(no_layer-1)*no_zu1 \ + no_time*no_layer/),graphic) else plot=new((/no_time*no_layer*no_var - no_time*(no_layer-1)*no_zu1/),graphic) end if dim_plot=dimsizes(plot) 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 if (xyc .EQ. 1)then if (sort .EQ. "time")then if(z_d(lo) .eq. -1.d) then level = "z-average" else level = "z=" + z_d(lo) + "m" end if else if(z_d(li) .eq. -1.d) then level = "z-average" else level = "z=" + z_d(li) + "m" end if end if end if if (xzc .EQ. 1)then if (sort .EQ. "time")then if(y_d(lo) .eq. -1.d) then level = "y-average" else level = "y=" + y_d(lo) + "m" end if else if(y_d(li) .eq. -1.d) then level = "y-average" else level = "y=" + y_d(li) + "m" end if end if end if if (yzc .EQ. 1)then if (sort .EQ. "time")then if(x_d(lo) .eq. -1.d) then level = "x-average" else level = "x=" + x_d(lo) + "m" end if else if(x_d(li) .eq. -1.d) then level = "x-average" else level = "x=" + x_d(li) + "m" end if end if end if 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@tmXBLabelFontHeightF = font_size vecres@tmYLLabelFontHeightF = font_size vecres@tiXAxisFontHeightF = font_size vecres@tiYAxisFontHeightF = font_size if (sort .EQ. "time")then vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h "+level else vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h "+level end if if (xyc .EQ. 1) then vecres@tiXAxisString = "x [m]" vecres@tiYAxisString = "y [m]" if (sort .EQ. "time")then plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres) else plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),vect2(lo,li-lis,:,:),vecres) end if end if if (xzc .EQ. 1) then vecres@tiXAxisString = "x [m]" vecres@tiYAxisString = "z [m]" if (sort .EQ. "time")then plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),vect2(li,:,lo-los,:),vecres) else plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),vect2(lo,:,li-lis,:),vecres) end if end if if (yzc .EQ. 1) then vecres@tiXAxisString = "y [m]" vecres@tiYAxisString = "z [m]" if (sort .EQ. "time")then plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),vect2(li,:,:,lo-los),vecres) else plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),vect2(lo,:,:,li-lis),vecres) end if end if 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@gsnLeftString = "Plot of "+vNam(varn) cs_res@gsnRightString = unit(varn) if ( sort .eq. "time" ) then if ( z_d(lo) .eq. -1)then if (delta_z .EQ. -1) then level = "-average" else level = "=" + z_d(lo) + "m" end if else level = "=" + z_d(lo) + "m" end if if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \ .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \ .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \ .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy" \ .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \ .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\ .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy" \ .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then loe = 0 level = "=" + zu1(0) + "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-los,:,:),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-los,:,:),vect2(li,lo-los,:,:),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if (z_d(li) .eq. -1) then if (delta_z .EQ. -1) then level = "-average" else level = "=" + z_d(li) + "m" end if else level = "=" + z_d(li) + "m" end if if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \ .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \ .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \ .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy" \ .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \ .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\ .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy" \ .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then lie = 0 level = "=" + zu1(0) + "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-lis,:,:),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-lis,:,:),vect2(lo,li-lis,:,:),vecres) overlay(plot(n), plot_vec) end if end if end if ; **************************************************** ; xz cross section ; **************************************************** if ( xzc .eq. 1 ) then cs_res@gsnLeftString = "Plot of "+vNam(varn) if ( sort .eq. "time" ) then if ( y_d(lo) .eq. -1 ) then level = "-average" else level = "=" + y_d(lo) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo-los,:),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-los,:),vect2(li,:,lo-los,:),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if ( y_d(li) .eq. -1 ) then level = "-average" else level = "=" + y_d(li) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li-lis,:),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-lis,:),vect2(lo,:,li-lis,:),vecres) overlay(plot(n), plot_vec) end if end if end if ; **************************************************** ; yz cross section ; **************************************************** if ( yzc .eq. 1 ) then cs_res@gsnLeftString = "Plot of "+vNam(varn) if ( sort .eq. "time" ) then if ( x_d(lo) .eq. -1 ) then level = "-average" else level = "=" + x_d(lo) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo-los),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-los),vect2(li,:,:,lo-los),vecres) overlay(plot(n), plot_vec) end if end if if ( sort .eq. "layer" ) then if ( x_d(li) .eq. -1 ) then level = "-average" else level = "=" + x_d(li) + "m" end if cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li-lis),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-lis),vect2(lo,:,:,li-lis),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 ; *************************************************** no_frames = 0 if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) 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,dim_plot-1,no_rows*no_columns if ( np + no_rows*no_columns .gt. dim_plot-1) then gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP) no_frames = no_frames + 1 else gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP) no_frames = no_frames + 1 end if end do end if else if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. dim_plot .gt. no_rows*no_columns) then gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP) print(" ") print("Outputs to .eps or .epsi have only one frame") print(" ") else do np = 0,dim_plot-1,no_rows*no_columns if ( np + no_rows*no_columns .gt. dim_plot-1) then gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP) no_frames = no_frames + 1 else gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP) no_frames = no_frames + 1 end if end do end if end if if (format_out .EQ. "png" ) then png_output = new((/no_frames/), string) j = 0 do i=0, no_frames-1 j = i + 1 if (j .LE. 9) then png_output(i) = file_out+".00000"+j+".png" end if if (j .GT. 9 .AND. j .LE. 99) then png_output(i) = file_out+".0000"+j+".png" end if if (j .GT. 99 .AND. j .LE. 999) then png_output(i) = file_out+".000"+j+".png" end if if (j .GT. 999) then png_output(i) = file_out+".00"+j+".png" end if ;using imagemagick's convert for reducing the white ;space around the plot cmd = "convert -geometry 1000x1000 -density 300 -trim " + \ png_output(i) + " " + png_output(i) system(cmd) end do print(" ") print("Output to: "+ png_output) print(" ") else print(" ") print("Output to: " + file_out +"."+ format_out) print(" ") end if end