; Example script to produce plots for a WRF real-data run, ; with the ARW coordinate dynamics option. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" undef("ARR2") function ARR2(a0:numeric,b0:numeric,tk:numeric) local xi begin xi = a0 * exp( -b0 /tk ) return(xi) end begin year = 2006 month = 7 do day = 18, 19 bl_thres1 = 70. bl_thres2 = 50. dirname = "/gpfs/proj2/dasg006/jeff/case_us_chem105_megan/" dirname = "/gpfs/proj2/dasg006/jeff/case_us_chem105-alt2/wrfout_all/" ; convert single digit months, days, and hours to double digits if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if ; create output file name type = "ps" fout = "Rates_BL"+bl_thres1+"_"+cmon+"-"+cday wks = gsn_open_wks(type,fout) ; Set some Basic Plot Information ; We would like to switch off the Initial Time and Footer info on the plot ; We also do not want titles on the color label bar res = True res@MainTitle = "" res@InitTime = False res@ValidTime = True res@Footer = False res@lbTitleOn = True res@lbTitleString = "ppbv" res@lbTitleFontHeightF = 0.015 pltres = True pltres1 = True pltres1@FramePlot = False mpres = True mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpGridLineColor = "Gray" mpres@mpLimbLineColor = "Gray" mpres@mpPerimLineColor = "Black" gsn_define_colormap(wks,"BlAqGrYeOrReVi200") opts = res opts@cnFillOn = True opts@gsnSpreadColorStart = 0 opts@cnLevelSelectionMode = "ExplicitLevels" opts@cnLevels = (/ 1.e3, 1.e4, 2.e4, 5.e4, 1.e5, 2.e5, 5.e5, 1.e6, 2.e6, 5.e6 /) opts_z = res opts_z@cnFillOn = False opts_z@gsnSpreadColorStart = 2 opts_z@cnLevelSelectionMode = "ExplicitLevels" opts_z@cnLevels = (/ 9700., 9750. /) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; The WRF ARW input files ifs = 0 ife = 0 ;7 ntms = ife-ifs+1 itime = new( (/ntms/), float) avg_utco_th2 = new( (/ntms/), float) avg_utco_th1 = new( (/ntms/), float) avg_utch4_th2 = new( (/ntms/), float) avg_utch4_th1 = new( (/ntms/), float) avg_uthc3_th2 = new( (/ntms/), float) avg_uthc3_th1 = new( (/ntms/), float) avg_utolt_th2 = new( (/ntms/), float) avg_utolt_th1 = new( (/ntms/), float) avg_utiso_th2 = new( (/ntms/), float) avg_utiso_th1 = new( (/ntms/), float) avg_uthcho_th2 = new( (/ntms/), float) avg_uthcho_th1 = new( (/ntms/), float) avg_utmacr_th2 = new( (/ntms/), float) avg_utmacr_th1 = new( (/ntms/), float) avg_utop1_th2 = new( (/ntms/), float) avg_utop1_th1 = new( (/ntms/), float) avg_utop2_th2 = new( (/ntms/), float) avg_utop2_th1 = new( (/ntms/), float) avg_uthket_th2 = new( (/ntms/), float) avg_uthket_th1 = new( (/ntms/), float) avg_utonit_th2 = new( (/ntms/), float) avg_utonit_th1 = new( (/ntms/), float) avg_utete_th2 = new( (/ntms/), float) avg_utete_th1 = new( (/ntms/), float) do ifl=ifs,ife ifn = floattointeger(floor(ifl/1)) ; divide by 1 for every 60 min hr = ifn*3 minutes = 0 ; ifl*30 - ifn*60 print(" ifn "+ifn+" minutes "+minutes+" hr "+hr) if(minutes .lt. 10) then minute = "0"+minutes else minute = minutes end if if(hr .lt. 10) then hour = "0"+hr else hour = hr end if if(hr .gt. 23) then day = day + 1 hr = hr-24 print(" ifn "+ifn+" minutes "+minutes+" hr "+hr+" day "+day) if(hr .lt. 10) then delete(hour) hour = "0"+hr else hour = hr end if end if if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if filename = "wrfout_d01_"+year+"-"+cmon+"-"+cday+"_"+hour+":"+minute+":00.nc" itime(ifl) = day + (hr*60. + minutes)/1440. fname1 = dirname+filename print(" reading from file "+fname1+" "+itime(ifl)) a = addfile(fname1, "r") times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) time = 0 res@TimeLabel = times(time) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; First get the variables we will need na = 6.022e23 ; molecules/mol mwa = 28.966 ; g/mol rd = 8.314e3/mwa ; J/kg/K xlat = a->XLAT(0,:,:) ; 2d arrays of latitude, longitude xlon = a->XLONG(0,:,:) pres = wrf_user_getvar(a,"pressure",time) ; hPa=mb tc = wrf_user_getvar(a,"tc",time) tk = tc + 273.15 rho = 100.*pres/(rd*tk) ; kg air/m3 mair = rho*1.e-6 *1.e3/mwa *na ; kgair/m3 m3/cm3 g/kg mol/g molecules/mol = molecules/cm3 z = wrf_user_getvar(a, "z",time) ; grid point height (m) print(" sfc avgs "+avg(pres(0,:,:))+" "+avg(tk(0,:,:))+" "+avg(rho(0,:,:))+" "+avg(mair(0,:,:))+" "+avg(z(0,:,:)) ) mdbz = wrf_user_getvar(a, "mdbz", time) ppm2cm3 = 1.e-6 /mwa*1.e3*rho*1.e-6*na ; molX/molair/ppm molair/gair gair/kgair kgair/m3 m3/cm3 molecules/molX = molec/cm3 oh = wrf_user_getvar(a,"ho",time) co = wrf_user_getvar(a,"co",time) ch4 = wrf_user_getvar(a,"ch4",time) hc3 = wrf_user_getvar(a,"hc3",time) olt = wrf_user_getvar(a,"olt",time) iso = wrf_user_getvar(a,"iso",time) hcho = wrf_user_getvar(a,"hcho",time) macr = wrf_user_getvar(a,"macr",time) op1 = wrf_user_getvar(a,"op1",time) op2 = wrf_user_getvar(a,"op2",time) hket = wrf_user_getvar(a,"hket",time) onit = wrf_user_getvar(a,"onit",time) ete = wrf_user_getvar(a,"ete",time) sr7 = wrf_user_getvar(a,"sr7",time) sr7 = sr7 * 100. ;----------------- ; convert to molec/cm3 oh = oh * ppm2cm3 co = co * ppm2cm3 ch4 = ch4 * ppm2cm3 hc3 = hc3 * ppm2cm3 olt = olt * ppm2cm3 iso = iso * ppm2cm3 hcho = hcho * ppm2cm3 macr = macr * ppm2cm3 op1 = op1 * ppm2cm3 op2 = op2 * ppm2cm3 hket = hket * ppm2cm3 onit = onit * ppm2cm3 ete = ete * ppm2cm3 dimo = dimsizes(oh) nz = dimo(0) ny = dimo(1) nx = dimo(2) ;---------------- ; oxidation with OH reactions that make a peroxy radical ;---------------- ; top 12 k58 = 1.5e-13*(1. + 2.439e-20*mair) ; CO+OH --> HO2 k61 = tk*tk *7.44e-18 * exp(-1361./tk) ; CH4+OH --> MO2 k63 = ARR2(5.26e-12,260.0,tk) ; HC3+OH --> 0.381 HO2 + 0.583 HC3P k66 = ARR2(1.96e-12,-438.0,tk) ; ETE+OH --> ETEP k67 = ARR2(5.72e-12,-500.0,tk) ; OLT+OH --> OLTP k70 = ARR2(2.54e-11,-410.0,tk) ; ISO+OH --> ISOP k76 = 1.00e-11 ; HCHO+OH--> HO2 k79 = 3.00e-12 ; HKET+OH--> HO2 + MGLY k82 = ARR2(1.86e-11,-175.0,tk) ; MACR+OH--> 0.49 HO2 k85 = ARR2(2.93e-12,-190.0,tk) ; OP1+OH --> 0.65 MO2 k86 = ARR2(3.40e-12,-190.0,tk) ; OP2+OH --> 0.44 HC3P k90 = ARR2(5.31e-12, 260.0,tk) ; ONIT+OH--> HC3P ; Loss Rates l58 = k58 * co*oh l61 = k61 * ch4*oh l63 = k63 * hc3*oh ; HC3+OH --> 0.381 HO2 + 0.583 HC3P l66 = k66 * ete*oh ; ETE+OH --> ETEP l67 = k67 * olt*oh l70 = k70 * iso*oh l76 = k76 * hcho*oh l79 = k79 * hket*oh l82 = k82 * macr*oh ; MACR+OH--> 0.49 HO2 l85 = k85 * op1*oh ; OP1+OH --> 0.65 MO2 l86 = k86 * op2*oh ; OP2+OH --> 0.44 HC3P l90 = k90 * onit*oh l58@description = "CO+OH rate" l61@description = "CH4+OH rate" l63@description = "HC3+OH rate" l66@description = "ETE+OH rate" l67@description = "OLT+OH rate" l70@description = "ISO+OH rate" l76@description = "HCHO+OH rate" l79@description = "HKET+OH rate" l82@description = "MACR+OH rate" l85@description = "OP1+OH rate" l86@description = "OP2+OH rate" l90@description = "ONIT+OH rate" l58@units = "molecules/cm3/s" l61@units = "molecules/cm3/s" l63@units = "molecules/cm3/s" l66@units = "molecules/cm3/s" l67@units = "molecules/cm3/s" l70@units = "molecules/cm3/s" l76@units = "molecules/cm3/s" l79@units = "molecules/cm3/s" l82@units = "molecules/cm3/s" l85@units = "molecules/cm3/s" l86@units = "molecules/cm3/s" l90@units = "molecules/cm3/s" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; get heights at 300 hPa level ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; z_plane = wrf_user_intrp3d( z, pres, "h", 300., 0., False) print(" max height at 300 mb "+max(z_plane) ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; get BL tracer at 10 km level ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; bl_plane = wrf_user_intrp3d( sr7, z, "h", 10000., 0., False) print(" max BL tracer at 10 km "+max(bl_plane) ) avg_z9750 = new( (/nz/), float) avg_z9750 = 0. i1_z300 = new( (/ny,nx/), integer) i2_z300 = new( (/ny,nx/), integer) i1_z300 = 0 ; -999 i2_z300 = 0 ; -999 ac = 0 bc = 0 if(any(bl_plane .ge. bl_thres2)) then do j=0,ny-1 do i=0,nx-1 if(bl_plane(j,i) .ge. bl_thres1) then i1_z300(j,i) = 1 ac = ac + 1 do k=0,nz-1 avg_z9750(k) = avg_z9750(k) + z(k,j,i) end do else if(bl_plane(j,i) .ge. bl_thres2 .and. bl_plane(j,i) .lt. bl_thres1) then i2_z300(j,i) = 1 bc = bc + 1 end if end if end do end do end if avg_z9750 = avg_z9750 / ac i1_z300@FillValue = -999 i2_z300@FillValue = -999 print(" number points with sr7>70, >50 "+ac+" "+bc +" avg z "+avg_z9750(30) ) z_9750 = z rco_9750 = l58 rco_9700 = l58 rch4_9750 = l61 rch4_9700 = l61 rhc3_9750 = l63 rhc3_9700 = l63 rolt_9750 = l67 rolt_9700 = l67 riso_9750 = l70 riso_9700 = l70 rhcho_9750 = l76 rhcho_9700 = l76 rmacr_9750 = l82 rmacr_9700 = l82 rop1_9750 = l85 rop1_9700 = l85 rop2_9750 = l86 rop2_9700 = l86 rhket_9750 = l79 rhket_9700 = l79 ronit_9750 = l90 ronit_9700 = l90 rete_9750 = l66 rete_9700 = l66 do k=0,nz-1 z_9750(k,:,:) = z(k,:,:)*i1_z300(:,:) rco_9750(k,:,:) = l58(k,:,:)*i1_z300(:,:) rco_9700(k,:,:) = l58(k,:,:)*i2_z300(:,:) rch4_9750(k,:,:) = l61(k,:,:)*i1_z300(:,:) rch4_9700(k,:,:) = l61(k,:,:)*i2_z300(:,:) rhc3_9750(k,:,:) = l63(k,:,:)*i1_z300(:,:) rhc3_9700(k,:,:) = l63(k,:,:)*i2_z300(:,:) rolt_9750(k,:,:) = l67(k,:,:)*i1_z300(:,:) rolt_9700(k,:,:) = l67(k,:,:)*i2_z300(:,:) riso_9750(k,:,:) = l70(k,:,:)*i1_z300(:,:) riso_9700(k,:,:) = l70(k,:,:)*i2_z300(:,:) rhcho_9750(k,:,:) = l76(k,:,:)*i1_z300(:,:) rhcho_9700(k,:,:) = l76(k,:,:)*i2_z300(:,:) rmacr_9750(k,:,:) = l82(k,:,:)*i1_z300(:,:) rmacr_9700(k,:,:) = l82(k,:,:)*i2_z300(:,:) rop1_9750(k,:,:) = l85(k,:,:)*i1_z300(:,:) rop1_9700(k,:,:) = l85(k,:,:)*i2_z300(:,:) rop2_9750(k,:,:) = l86(k,:,:)*i1_z300(:,:) rop2_9700(k,:,:) = l86(k,:,:)*i2_z300(:,:) rhket_9750(k,:,:) = l79(k,:,:)*i1_z300(:,:) rhket_9700(k,:,:) = l79(k,:,:)*i2_z300(:,:) ronit_9750(k,:,:) = l90(k,:,:)*i1_z300(:,:) ronit_9700(k,:,:) = l90(k,:,:)*i2_z300(:,:) rete_9750(k,:,:) = l66(k,:,:)*i1_z300(:,:) rete_9700(k,:,:) = l66(k,:,:)*i2_z300(:,:) end do if(ifl .eq. 0) then print(" assign 10 km planes ") co_plane = wrf_user_intrp3d( rco_9750, z, "h", 10000., 0., False) ch4_plane = wrf_user_intrp3d( rch4_9750, z, "h", 10000., 0., False) hc3_plane = wrf_user_intrp3d( rhc3_9750, z, "h", 10000., 0., False) olt_plane = wrf_user_intrp3d( rolt_9750, z, "h", 10000., 0., False) iso_plane = wrf_user_intrp3d( riso_9750, z, "h", 10000., 0., False) hcho_plane = wrf_user_intrp3d( rhcho_9750,z, "h", 10000., 0., False) macr_plane = wrf_user_intrp3d( rmacr_9750,z, "h", 10000., 0., False) op1_plane = wrf_user_intrp3d( rop1_9750, z, "h", 10000., 0., False) op2_plane = wrf_user_intrp3d( rop2_9750, z, "h", 10000., 0., False) hket_plane = wrf_user_intrp3d( rhket_9750,z, "h", 10000., 0., False) onit_plane = wrf_user_intrp3d( ronit_9750,z, "h", 10000., 0., False) ete_plane = wrf_user_intrp3d( rete_9750, z, "h", 10000., 0., False) print(" start plotting ") contour_z = wrf_contour(a,wks,z_plane,opts_z) contour_co = wrf_contour(a,wks,co_plane(:,:),opts) contour_ch4 = wrf_contour(a,wks,ch4_plane(:,:),opts) contour_hc3 = wrf_contour(a,wks,hc3_plane(:,:),opts) contour_olt = wrf_contour(a,wks,olt_plane(:,:),opts) contour_iso = wrf_contour(a,wks,iso_plane(:,:),opts) contour_hcho = wrf_contour(a,wks,hcho_plane(:,:),opts) contour_macr = wrf_contour(a,wks,macr_plane(:,:),opts) contour_op1 = wrf_contour(a,wks,op1_plane(:,:),opts) contour_op2 = wrf_contour(a,wks,op2_plane(:,:),opts) contour_hket = wrf_contour(a,wks,hket_plane(:,:),opts) contour_onit = wrf_contour(a,wks,onit_plane(:,:),opts) contour_ete = wrf_contour(a,wks,ete_plane(:,:),opts) plot = wrf_map_overlays(a,wks,(/contour_co,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_ch4,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_hc3,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_olt,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_iso,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_hcho,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_macr,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_op1,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_op2,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_hket,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_onit,contour_z/),pltres,mpres) plot = wrf_map_overlays(a,wks,(/contour_ete,contour_z/),pltres,mpres) print(" plotted high BL region ") end if k_10to12 = new( (/nz/), integer) k_10to12 = 0 kc = 0 do k=0,nz-1 if(avg_z9750(k) .gt. 9500. .and. avg_z9750(k) .lt. 12500.) then kc = kc + 1 k_10to12(k) = 1 print(" k, avg z "+k+" "+avg(z(k,:,:))+" "+avg_z9750(k) ) end if end do print(" number of model levels between 9500 and 12500 "+kc) rco_9750@_FillValue = 0. rco_9700@_FillValue = 0. rch4_9750@_FillValue = 0. rch4_9700@_FillValue = 0. rhc3_9750@_FillValue = 0. rhc3_9700@_FillValue = 0. rolt_9750@_FillValue = 0. rolt_9700@_FillValue = 0. riso_9750@_FillValue = 0. riso_9700@_FillValue = 0. rhcho_9750@_FillValue = 0. rhcho_9700@_FillValue = 0. rmacr_9750@_FillValue = 0. rmacr_9700@_FillValue = 0. rop1_9750@_FillValue = 0. rop1_9700@_FillValue = 0. rop2_9750@_FillValue = 0. rop2_9700@_FillValue = 0. rhket_9750@_FillValue = 0. rhket_9700@_FillValue = 0. ronit_9750@_FillValue = 0. ronit_9700@_FillValue = 0. rete_9750@_FillValue = 0. rete_9700@_FillValue = 0. n1 = num(.not.ismissing(rco_9750)) n2 = num(.not.ismissing(rco_9700)) print(" number of grid points in rco_9750 "+n1+" "+n2 ) ut_co_9750 = rco_9750 ut_co_9700 = rco_9700 ut_ch4_9750 = rch4_9750 ut_ch4_9700 = rch4_9700 ut_hc3_9750 = rhc3_9750 ut_hc3_9700 = rhc3_9700 ut_olt_9750 = rolt_9750 ut_olt_9700 = rolt_9700 ut_iso_9750 = riso_9750 ut_iso_9700 = riso_9700 ut_hcho_9750 = rhcho_9750 ut_hcho_9700 = rhcho_9700 ut_macr_9750 = rmacr_9750 ut_macr_9700 = rmacr_9700 ut_op1_9750 = rop1_9750 ut_op1_9700 = rop1_9700 ut_op2_9750 = rop2_9750 ut_op2_9700 = rop2_9700 ut_hket_9750 = rhket_9750 ut_hket_9700 = rhket_9700 ut_onit_9750 = ronit_9750 ut_onit_9700 = ronit_9700 ut_ete_9750 = rete_9750 ut_ete_9700 = rete_9700 do k=0,nz-1 ut_co_9750(k,:,:) = rco_9750(k,:,:)*k_10to12(k) ut_co_9700(k,:,:) = rco_9700(k,:,:)*k_10to12(k) ut_ch4_9750(k,:,:) = rch4_9750(k,:,:)*k_10to12(k) ut_ch4_9700(k,:,:) = rch4_9700(k,:,:)*k_10to12(k) ut_hc3_9750(k,:,:) = rhc3_9750(k,:,:)*k_10to12(k) ut_hc3_9700(k,:,:) = rhc3_9700(k,:,:)*k_10to12(k) ut_olt_9750(k,:,:) = rolt_9750(k,:,:)*k_10to12(k) ut_olt_9700(k,:,:) = rolt_9700(k,:,:)*k_10to12(k) ut_iso_9750(k,:,:) = riso_9750(k,:,:)*k_10to12(k) ut_iso_9700(k,:,:) = riso_9700(k,:,:)*k_10to12(k) ut_hcho_9750(k,:,:) = rhcho_9750(k,:,:)*k_10to12(k) ut_hcho_9700(k,:,:) = rhcho_9700(k,:,:)*k_10to12(k) ut_macr_9750(k,:,:) = rmacr_9750(k,:,:)*k_10to12(k) ut_macr_9700(k,:,:) = rmacr_9700(k,:,:)*k_10to12(k) ut_op1_9750(k,:,:) = rop1_9750(k,:,:)*k_10to12(k) ut_op1_9700(k,:,:) = rop1_9700(k,:,:)*k_10to12(k) ut_op2_9750(k,:,:) = rop2_9750(k,:,:)*k_10to12(k) ut_op2_9700(k,:,:) = rop2_9700(k,:,:)*k_10to12(k) ut_hket_9750(k,:,:) = rhket_9750(k,:,:)*k_10to12(k) ut_hket_9700(k,:,:) = rhket_9700(k,:,:)*k_10to12(k) ut_onit_9750(k,:,:) = ronit_9750(k,:,:)*k_10to12(k) ut_onit_9700(k,:,:) = ronit_9700(k,:,:)*k_10to12(k) ut_ete_9750(k,:,:) = rete_9750(k,:,:)*k_10to12(k) ut_ete_9700(k,:,:) = rete_9700(k,:,:)*k_10to12(k) end do ut_co_9750@_FillValue = 0. ut_co_9700@_FillValue = 0. ut_ch4_9750@_FillValue = 0. ut_ch4_9700@_FillValue = 0. ut_hc3_9750@_FillValue = 0. ut_hc3_9700@_FillValue = 0. ut_olt_9750@_FillValue = 0. ut_olt_9700@_FillValue = 0. ut_iso_9750@_FillValue = 0. ut_iso_9700@_FillValue = 0. ut_hcho_9750@_FillValue = 0. ut_hcho_9700@_FillValue = 0. ut_macr_9750@_FillValue = 0. ut_macr_9700@_FillValue = 0. ut_op1_9750@_FillValue = 0. ut_op1_9700@_FillValue = 0. ut_op2_9750@_FillValue = 0. ut_op2_9700@_FillValue = 0. ut_hket_9750@_FillValue = 0. ut_hket_9700@_FillValue = 0. ut_onit_9750@_FillValue = 0. ut_onit_9700@_FillValue = 0. ut_ete_9750@_FillValue = 0. ut_ete_9700@_FillValue = 0. print(" min/max ut_rate co_th2 "+min(ut_co_9750)+" "+max(ut_co_9750) ) avg_utco_th2(ifl) = avg(ut_co_9750) avg_utco_th1(ifl) = avg(ut_co_9700) avg_utch4_th2(ifl) = avg(ut_ch4_9750) avg_utch4_th1(ifl) = avg(ut_ch4_9700) avg_uthc3_th2(ifl) = avg(ut_hc3_9750) avg_uthc3_th1(ifl) = avg(ut_hc3_9700) avg_utolt_th2(ifl) = avg(ut_olt_9750) avg_utolt_th1(ifl) = avg(ut_olt_9700) avg_utiso_th2(ifl) = avg(ut_iso_9750) avg_utiso_th1(ifl) = avg(ut_iso_9700) avg_uthcho_th2(ifl) = avg(ut_hcho_9750) avg_uthcho_th1(ifl) = avg(ut_hcho_9700) avg_utmacr_th2(ifl) = avg(ut_macr_9750) avg_utmacr_th1(ifl) = avg(ut_macr_9700) avg_utop1_th2(ifl) = avg(ut_op1_9750) avg_utop1_th1(ifl) = avg(ut_op1_9700) avg_utop2_th2(ifl) = avg(ut_op2_9750) avg_utop2_th1(ifl) = avg(ut_op2_9700) avg_uthket_th2(ifl) = avg(ut_hket_9750) avg_uthket_th1(ifl) = avg(ut_hket_9700) avg_utonit_th2(ifl) = avg(ut_onit_9750) avg_utonit_th1(ifl) = avg(ut_onit_9700) avg_utete_th2(ifl) = avg(ut_ete_9750) avg_utete_th1(ifl) = avg(ut_ete_9700) n1 = num(.not.ismissing(ut_co_9750)) n2 = num(.not.ismissing(ut_co_9700)) print(" average UT O3 for z>th2, th1 grid points in 10-12 km "+avg_utco_th2(ifl)+" "+avg_utco_th1(ifl) ) print(" number of grid points averaged "+n1+" "+n2 ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Clean up delete(hr) delete(minutes) delete(hour) delete(minute) end do ; ifl loop data_1 = new( (/ntms,17/), float) data_2 = new( (/ntms,17/), float) data_1(:,0) = itime data_1(:,1) = avg_utco_th2 data_1(:,2) = avg_utch4_th2 data_1(:,3) = avg_uthc3_th2 data_1(:,4) = avg_utolt_th2 data_1(:,5) = avg_utiso_th2 data_1(:,6) = avg_uthcho_th2 data_1(:,7) = avg_utmacr_th2 data_1(:,8) = avg_utop1_th2 data_1(:,9) = avg_utop2_th2 data_1(:,10) = avg_uthket_th2 data_1(:,11) = avg_utonit_th2 data_1(:,12) = avg_utete_th2 data_2(:,0) = itime data_2(:,1) = avg_utco_th1 data_2(:,2) = avg_utch4_th1 data_2(:,3) = avg_uthc3_th1 data_2(:,4) = avg_utolt_th1 data_2(:,5) = avg_utiso_th1 data_2(:,6) = avg_uthcho_th1 data_2(:,7) = avg_utmacr_th1 data_2(:,8) = avg_utop1_th1 data_2(:,9) = avg_utop2_th1 data_2(:,10) = avg_uthket_th1 data_2(:,11) = avg_utonit_th1 data_2(:,12) = avg_utete_th1 opt1 = True opt1@fout = "./rates_BL"+bl_thres1+"_"+cmon+"-"+cday+".dat" write_matrix(data_1, "13e13.6", opt1) opt1@fout = "./rates_BL"+bl_thres2+"_"+cmon+"-"+cday+".dat" write_matrix(data_2, "13e13.6", opt1) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; day do loop delete(day) delete(cday) end