;*********************************************** ; write_radar_qvp.ncl ; writes QVP diagrams based on radar at all available levels. ; ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" undef("set_res1_list") function set_res1_list() begin res = True res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnLevels = ispan(-9,48,3) ; set the contour levels cmap := read_colormap_file("NCV_rainbow2") printVarSummary(cmap) cmap := cmap(50:, :) res@cnFillPalette = cmap ; set color map ; res@pmLabelBarOrthogonalPosF = -0.1 res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Distance (km)" ; opts_fft@lbLabelBarOn = False ; res@gsnDraw = False ; don't draw ; res@gsnFrame = False ; don't advance frame res@cnLinesOn = False ; this gives a warning (once) that it's not a valid setting, but it is. ; I think it happens because I reuse the opts for different variables. res@cnLineLabelsOn = False ; opts_fft@lbLabelBarOn = False res@cnInfoLabelOn = False return(res) end undef("set_res2_list") function set_res2_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap res@cnLevels = (/-2.0,-1.0,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ 2.5,3.0/) res@cnFillColors = (/0,90,85,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Distance (km)" return(res) end undef("set_res3_list") function set_res3_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnFillPalette = "MPL_gist_ncar" ; set color map ; printVarSummary(cmap) ; cmap := cmap(75:, :) ; dbres@cnFillPalette := cmap ; set color map res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; dfres@cnLevels = ispan(-2,6,1) ; set the contour levels res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res@cnLevels = (/ 0.70,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.80,0.81,\ 0.82,0.83,0.84,0.85,0.86,0.87,0.88,\ 0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\ 0.97,0.98,0.99,1.0/) res@cnFillColors = (/ 13,17,21,25,29,33,37,41,\ 45,49,53,57,61,65,69,73,75,77,\ 81 ,83 ,85,\ 87,89,91,93,95,126,122,117,115,113,111,110/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Distance (km)" return(res) end undef("set_res4_list") function set_res4_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnFillPalette = "NCV_rainbow2" res@cnLevels = (/-1.0,-0.5,-0.25,0,0.01,0.02,0.06,0.1,0.2,0.3,0.4,0.5,0.6,\ 0.7,0.8,0.9,1.0/) res@cnFillColors = (/0,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Distance (km)" return(res) end external RADAR_WRITE_QVP "./radar_write_qvp.so" begin ;main program kdp_find = True ; kdp_find = False filenames = ((/"KOKX_V06_20130209_000138.nc","KOKX_V06_20130209_003053.nc", \ "KOKX_V06_20130209_010213.nc", "KOKX_V06_20130209_013125.nc", \ "KOKX_V06_20130209_020046.nc", "KOKX_V06_20130209_023001.nc", \ "KOKX_V06_20130209_025915.nc","KOKX_V06_20130209_032833.nc",\ "KOKX_V06_20130209_035753.nc"/)) n_files = dimsizes(filenames) filename = filenames(0) f = addfile(filename,"r") var := f->Reflectivity(:,:,:) dimsvar = dimsizes(var) ; n_r = dimsvar(2) n_e = dimsvar(0) n_r = 3*dimsvar(2)/6 ; n_z = 500 n_z = 101 dbz_z = new((/n_e+2,n_z/),float) zdr_z = new((/n_e+2,n_z/),float) cc_z = new((/n_e+2,n_z/),float) kdp_z = new((/n_e+2,n_z/),float) dbz_1d = new((/n_z/),float) zdr_1d = new((/n_z/),float) cc_1d = new((/n_z/),float) kdp_1d = new((/n_z/),float) z_z = new((/n_files,n_z/),float) dbz_zr = new((/n_e+2,n_r/),float) zdr_zr = new((/n_e+2,n_r/),float) cc_zr = new((/n_e+2,n_r/),float) kdp_zr = new((/n_e+2,n_r/),float) z_zr= new((/n_e+2,n_r/),float) dbz_z = 0 zdr_z = 0 cc_z = 0 kdp_z = 0 z_z = 0. dbz_zr = -999 zdr_zr = -999 cc_zr = -999 kdp_zr = -999 z_zr = 0. dbz_z@_FillValue = -999 zdr_z@_FillValue = -999 cc_z@_FillValue = -999 kdp_z@_FillValue = -999 z_10m = new(n_z,float) z_10m(0) = 0 do k=1,n_z-1 ; z_10m(k) = z_10m(k-1)+10 z_10m(k) = z_10m(k-1)+50 end do yr = str_get_cols(filenames, 9, 12) year = stringtoint(yr) print("year = " + year) mnth = str_get_cols(filenames, 13, 14) month = stringtoint(mnth) print("month = " + month) dy = str_get_cols(filenames, 15, 16) day = stringtoint(dy) print("day = " + day) hr = str_get_cols(filenames, 18, 19) hour = stringtoint(hr) print("hour = " + hour) mn = str_get_cols(filenames, 20, 21) minute = stringtoint(mn) print("minute = " + minute) do i_file = 0,n_files-1,1 ; do i_file = 1,1,1 ; do i_file = 0,8,1 ; do i_file = 3,3,1 ; do i_file = 8,8,4 print("i_file = " + i_file) filename = filenames(i_file) f = addfile(filename,"r") print("i_file = " + i_file) print("filename = " + filename) do i_elev = 0,n_e-1 ; do i_elev = 0,0,1 ; do i_elev= 0,n_elev(0)-1,1 print("i_elev = " + i_elev) elev := f->elevationR ; printVarSummary(elev) var := f->Reflectivity(i_elev,:,:) ; printVarSummary(var) ; print(X) dimsvar2 = dimsizes(var) ; print("dimsvar = " + dimsvar2) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data angless := f->azimuthR angles := f->azimuthR(i_elev,:) dims_ang = dimsizes(angles) ; print("dims_ang = " + dims_ang) j_s = 0 j_e = dims_ang(0)-1 dsizes_dbz = dimsizes(dbz) i_s = 0 i_e = 3*dsizes_dbz(1)/6-1 print("i_e = " + i_e) dbz_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(dbz_sub) var := f->DifferentialReflectivity(i_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data zdr_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(zdr_sub) ;; Correlation Coefficient var := f->CorrelationCoefficient(i_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz := (dbz * scale_factor)+add_offset ;; now scale the data cc_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(cc_sub) ; Differential Phase var := f->DifferentialPhase(i_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data wk_smooth121(dbz) wk_smooth121(dbz) dfp_o := dbz dfp_sub := dbz(j_s:j_e,i_s:i_e) dfp_sub@_FillValue = default_fillvalue(typeof(dfp_sub)) replace_ieeenan (dfp_sub,dfp_sub@_FillValue,0) ; calculate KDP ; filtered dp dfp_rave := dfp_o dfp_rave = 0. dims_o = dimsizes(dfp_o) do j = 0,dims_o(0)-1 dfp_rave(j,:) = runave(dfp_o(j,:),18,0) end do if (kdp_find)then kdp := dfp_o kdp = 0. ; define variable to hold points for determing slope of "line." ; print("calculating kdp") dfp_x = new(9,float) dfp_y = new(9,float) do j = 0,dims_o(0)-1 do i = 5,dims_o(1)-5 n_dfp = 0 check = ismissing(dfp_rave(j,i)) if (check.eq.False)then ; check = ismissing(cc_o(j,i)) ; if (check.eq.False)then ; if (cc_o(j,i).ge.0.9)then do ii = -4 + i,4+i dfp_x(n_dfp) = ii dfp_y(n_dfp) = dfp_rave(j,ii) n_dfp = n_dfp + 1 end do rc = regline(dfp_x,dfp_y) kdp(j,i) = rc ; end if ; end if end if end do end do kdp_sub := kdp(j_s:j_e,i_s:i_e) else kdp_sub = dfp_sub end if ; printVarSummary(kdp_sub) ;---To see a list of all the variables on the file ; angles=90-angles angles_sub := angles(j_s:j_e) DEGTORAD = 0.017453292519943 xcenter = 0.0 ycenter = 0.0 ; Create 2D coordinate arrays. ; x:=f->distanceR x_r:=f->distanceR printVarSummary(x_r) ; print("x_r = " + x_r) ; print(X) ; x_sub = x(i_s:i_e) ; angles2d := conform(dbz_sub,angles_sub,0) ; x2d := conform(dbz_sub,x_sub,1)/1000 ; convert to km ; xarr := xcenter + x2d * cos(DEGTORAD * angles2d) ; yarr := ycenter + x2d * sin(DEGTORAD * angles2d) ; one_degree = 0.0174533 ; radians ; z_zr2 := dbz_sub x_xr2 := dbz_sub x_w := x_r/1000. do i = 0,i_e radius = x_r(i) ; this is the hypotenuse; ; see:https://www.mathsisfun.com/sine-cosine-tangent.html do i_ang = 0,dimsvar(1)-1 f_deg = one_degree*elev(i_elev,i_ang) sin_f = sin(f_deg) ; opposite over hypotenuse loc_hgt = radius*sin_f z_zr2(i_ang,i) = loc_hgt cos_f = cos(f_deg) ; adjacent over hypotenuse loc_dist = radius*cos_f x_xr2(i_ang,i) = loc_dist ; print("tan_f = " + tan_f) end do end do xr1 = dim_avg_n_Wrap(x_xr2,0) do i = 0,i_e radius = xr1(i)/1000. if (radius.le.50-10)then x_w(i) = 1. end if if (radius.gt.50-1.and.radius.le.50)then x_w(i) = 1./(radius - (50-1)) end if if (radius.gt.50)then x_w(i) = 1./(radius - (50-1))^2 end if end do x_w_sub = x_w(i_s:i_e) ; print("x_w = " + x_w) ; multiple by weighting function do i = 0,i_e dbz_sub(i_elev,i) = dbz_sub(i_elev,i)*x_w_sub(i) zdr_sub(i_elev,i) = zdr_sub(i_elev,i)*x_w_sub(i) cc_sub(i_elev,i) = cc_sub(i_elev,i)*x_w_sub(i) kdp_sub(i_elev,i) = kdp_sub(i_elev,i)*x_w_sub(i) end do ; average over azimuth (angles) ; radius = x_r(i_e) ; z_zr(i_elev,i_e) = dim_sum_n_Wrap( dbz_zr1 = dim_avg_n_Wrap(dbz_sub,0) zdr_zr1 = dim_avg_n_Wrap(zdr_sub,0) cc_zr1 = dim_avg_n_Wrap(cc_sub,0) kdp_zr1 = dim_avg_n_Wrap(kdp_sub,0) zr1 = dim_avg_n_Wrap(z_zr2,0) ; printVarSummary(dbz_zr1) ; printVarSummary(dbz_zr) dbz_zr(i_elev,:) = dbz_zr1 zdr_zr(i_elev,:) = zdr_zr1 cc_zr(i_elev,:) = cc_zr1 kdp_zr(i_elev,:) = kdp_zr1 ; printVarSummary(zr1) z_zr(i_elev,:) = zr1 ; printMinMax(dbz_zr1,False) ; printMinMax(zdr_zr1,False) ; printMinMax(cc_zr1,False) ; printMinMax(kdp_zr1,False) ; printMinMax(zr1,False) ; printVarSummary(dbz_zr1) ; printMinMax(dbz_zr,False) ; print("kdp_zr = " + kdp_zr(0,:)) ; print("z_zr2 = " + z_zr2(i_elev,:)) ; linear interpolation xi=z_zr(i_elev,:) fi=dbz_zr(i_elev,:) xo=z_10m dbz_z(i_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) ; print("f1 = " + dbz_z(i_elev,:)) fi=zdr_zr(i_elev,:) zdr_z(i_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) ; print("f2 = " + zdr_z(i_elev,:)) fi=cc_zr(i_elev,:) cc_z(i_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) ; print("f3 = " + cc_z(i_elev,:)) fi=kdp_zr(i_elev,:) kdp_z(i_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) ; print("f4 = " + kdp_z(i_elev,:)) ; do i = 0,i_e end do ; Start "HI" resolution part of radar observations ii_elev = -1 do i_elev = 0,3,3 ii_elev = ii_elev + 1 iii_elev = n_e+ii_elev elev := f->elevationR_HI ; printVarSummary(elev) var := f->Reflectivity_HI(i_elev,:,:) ; printVarSummary(var) ; print(X) dimsvar2 = dimsizes(var) ; print("dimsvar = " + dimsvar2) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data angless := f->azimuthR_HI angles := f->azimuthR_HI(ii_elev,:) dims_ang = dimsizes(angles) ; print("dims_ang = " + dims_ang) j_s = 0 j_e = dims_ang(0)-1 dsizes_dbz = dimsizes(dbz) ; i_s = 0 ; i_e = 3*dsizes_dbz(1)/6-1 dbz_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(dbz_sub) var := f->DifferentialReflectivity_HI(ii_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data zdr_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(zdr_sub) ;; Correlation Coefficient var := f->CorrelationCoefficient_HI(ii_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz := (dbz * scale_factor)+add_offset ;; now scale the data cc_sub := dbz(j_s:j_e,i_s:i_e) ; printVarSummary(cc_sub) ; Differential Phase var := f->DifferentialPhase_HI(ii_elev,:,:) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data wk_smooth121(dbz) wk_smooth121(dbz) dfp_o := dbz dfp_sub := dbz(j_s:j_e,i_s:i_e) dfp_sub@_FillValue = default_fillvalue(typeof(dfp_sub)) replace_ieeenan (dfp_sub,dfp_sub@_FillValue,0) ; calculate KDP ; filtered dp dfp_rave := dfp_o dfp_rave = 0. dims_o = dimsizes(dfp_o) do j = 0,dims_o(0)-1 dfp_rave(j,:) = runave(dfp_o(j,:),18,0) end do x_r:=f->distanceR_HI ; printVarSummary(x_r) if (kdp_find)then kdp := dfp_o kdp = 0. ; define variable to hold points for determing slope of "line." ; print("calculating kdp") dfp_x = new(9,float) dfp_y = new(9,float) do j = 0,dims_o(0)-1 do i = 5,dims_o(1)-5 n_dfp = 0 check = ismissing(dfp_rave(j,i)) if (check.eq.False)then ; check = ismissing(cc_o(j,i)) ; if (check.eq.False)then ; if (cc_o(j,i).ge.0.9)then do ii = -4 + i,4+i dfp_x(n_dfp) = ii dfp_y(n_dfp) = dfp_rave(j,ii) n_dfp = n_dfp + 1 end do rc = regline(dfp_x,dfp_y) kdp(j,i) = rc ; end if ; end if end if end do end do kdp_sub := kdp(j_s:j_e,i_s:i_e) else kdp_sub := dfp_sub end if ; printVarSummary(kdp_sub) ;---To see a list of all the variables on the file ; one_degree = 0.0174533 ; radians ; z_zr2 := dbz_sub x_xr2 := dbz_sub do i = 0,i_e radius = x_r(i) ; this is the hypotenuse; ; see:https://www.mathsisfun.com/sine-cosine-tangent.html do i_ang = 0,dimsvar2(0)-1 f_deg = one_degree*elev(i_elev,i_ang) sin_f = sin(f_deg) ; opposite over hypotenuse loc_hgt = radius*sin_f z_zr2(i_ang,i) = loc_hgt cos_f = cos(f_deg) ; adjacent over hypotenuse loc_dist = radius*cos_f x_xr2(i_ang,i) = loc_dist ; print("tan_f = " + tan_f) end do end do ; Calculate weighting function do i = 0,i_e radius = xr1(i)/1000. if (radius.le.50-10)then x_w(i) = 1. end if if (radius.gt.50-1.and.radius.le.50)then x_w(i) = 1./(radius - (50-1)) end if if (radius.gt.50)then x_w(i) = 1./(radius - (50-1))^2 end if end do ; print("x_w = " + x_w) x_w_sub = x_w(i_s:i_e) ; radius = x_r(i_e) ; print("second use of dbz") ; printVarSummary(dbz_sub) ; multiply by weighting function do i = 0,i_e dbz_sub(i_elev,i) = dbz_sub(ii_elev,i)*x_w_sub(i) zdr_sub(ii_elev,i) = zdr_sub(ii_elev,i)*x_w_sub(i) cc_sub(ii_elev,i) = cc_sub(ii_elev,i)*x_w_sub(i) kdp_sub(ii_elev,i) = kdp_sub(ii_elev,i)*x_w_sub(i) end do ; average around the azimuth (all angles) dbz_zr1 = dim_avg_n_Wrap(dbz_sub,0) zdr_zr1 = dim_avg_n_Wrap(zdr_sub,0) cc_zr1 = dim_avg_n_Wrap(cc_sub,0) kdp_zr1 = dim_avg_n_Wrap(kdp_sub,0) ; zr1 = dim_sum_n_Wrap(z_zr2,0) zr1 = dim_avg_n_Wrap(z_zr2,0) xr1 = dim_avg_n_Wrap(x_xr2,0) dbz_zr(iii_elev,:) = dbz_zr1 zdr_zr(iii_elev,:) = zdr_zr1 cc_zr(iii_elev,:) = cc_zr1 kdp_zr(iii_elev,:) = kdp_zr1 ; printVarSummary(zr1) z_zr(iii_elev,:) = zr1 ; linear interpolation xi=z_zr(iii_elev,:) fi=dbz_zr(iii_elev,:) xo=z_10m dbz_z(iii_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) fi=zdr_zr(iii_elev,:) zdr_z(iii_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) fi=cc_zr(iii_elev,:) cc_z(iii_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) fi=kdp_zr(iii_elev,:) kdp_z(iii_elev,:) = linint1_Wrap (xi,fi, False, xo, 0) ; print("f8 = " + kdp_z(iii_elev,:)) ; do i = 0,i_e end do ; average over elevations ; I am not sure why the 0 position is not defined. Perhaps there is no data. dbz_1d(0) = 0 dbz_1d = dim_avg_n_Wrap(dbz_z,0) ; print("dbz_1d = " + dbz_1d) zdr_1d = dim_avg_n_Wrap(zdr_z,0) zdr_1d(0) = 0 ; print("zdr_1d = " + zdr_1d) cc_1d = dim_avg_n_Wrap(cc_z,0) cc_1d(0) = 0 ; print("cc_1d = " + cc_1d) kdp_1d = dim_avg_n_Wrap(kdp_z,0) kdp_1d(0) = 0 print("kdp_1d(0) = " + kdp_1d(0)) ; write out data RADAR_WRITE_QVP::radar_write_qpv(year(i_file),month(i_file),day(i_file),hour(i_file),minute(i_file),\ dbz_1d,zdr_1d,cc_1d,kdp_1d,z_10m,n_z) end do ; filename end ; main program