pro read_const,filename,spec,rec,lat,lon,plevs,arr,tplevs=tplevs,thlevs=thlevs,$ type=type,nohdr=nohdr,zm=zm,rec_init=rec_init,ddtype=ddtype,reduce=reduce,trpp=trpp ; New optional keyword 'trpp' to read in tropopause pressure (from idaily only) ; ; New optional keyword: tplevs. true 3D pressure levels. It is always read for ; non-monthly averaged data, but only returned when tplevs is set. ; However, when thlevs is used, tplevs must also be used to ; calculate true pressures and return accurate thetas for the UTLS. ; new option keyword: reduce. Set it equal to the number of days per month ; you want returned. Only applies to idaily and overpass data types. Set ; reduce equal to 1 to 15. (reduce=15 means get every other day). ; Basic reader for any GMI const output. Needs to read in lat, long, and plevs. ; You must know the file name/directory of the data and the species ; number. (Get that from @get_species_combo, for example.) Option to interpolate ; to theta (an input). ; Must input how many time steps (records) to read in. ; Default scaling is to ppb. ; ; **** Creating a new optional keyword 'nohdr', to be defined in cases where the ; file read in does not contain the usual met data (e.g., plevs). Set 'nohdr' equal ; to the number of vertical levels in the run. This will identify the model grid ; being read in. Choices are 28, 33, and 42 for the GMI strat model (as of 2005). ; ************ NEW GEOS5 Aura choice: 72 levels. **************************** ; ; Use keyword 'type' if the field's name is anything but 'const'. (e.g., 'kel' or ; 'const_freq2'. ; ; Optional rec_init allows you to read in data starting from somewhere other than ; the beginning of the file. Set rec_init equal to rec number to start from ; Keyword 'rec' should be the total number of records to be read in. if keyword_set(nohdr) eq 0 then nohdr=0 if nohdr eq 23 then begin ; assume GISS II' 4x5 model im=72 & jm=46 & km=nohdr lat=indgen(jm)*4-90 & lon=indgen(im)*5 plevs=[987.770, 959.742, 921.013, 864.958, 776.289, 649.401, 504.167, 373.711, 281.984, 219.814, $ 171.403, 133.50, 101.60, 71.20, 43.90, 24.7, 13.8999, 7.31502, 3.04504, 0.960541, 0.303009, $ 0.0880737, 0.0166016] endif if nohdr eq 28 then begin ; assume 4x5 model, e.g. wmo_fvg im=72 & jm=46 & km=nohdr lat=indgen(jm)*4-90 & lon=indgen(im)*5 plevs=[918.316, 772.170, 649.033, 543.019, 441.362, 348.182,$ 276.080, 224.975, 188.703, 161.825, 139.835, 118.835, 100.947, 87.8748, 77.2146, $ 67.1829, 56.6158, 44.3917, 31.6228, 21.5444, 14.6780, 10.0000, 6.81291, 4.64159, $ 3.16228, 2.15443, 1.33353, 0.656174] endif ; if nohdr eq 33 then begin ; assume 2x2.5 model, e.g. hindcast runs im=144 & jm=91 & km=nohdr lat=indgen(jm)*2-90 & lon=findgen(im)*2.5 plevs=[918.316, 772.170, 649.033, 543.019, 441.362, 348.182,$ 276.080, 224.975, 188.703, 161.825, 139.835, 118.835, 100.947, 87.8748, 77.2146, $ 67.1829, 56.6158, 44.3917, 31.6228, 21.5444, 14.6780, 10.0000, 6.81291, 4.64159, $ 3.16228, 2.15443, 1.33353, 0.681290,0.316223, 0.146774, 0.0681520, 0.0316162, 0.0146637] endif ; if nohdr eq 42 then begin ; assume 2x2.5 model, e.g. runs with 5-yr fvgcm winds. im=144 & jm=91 & km=nohdr lat=indgen(jm)*2-90 & lon=findgen(im)*2.5 plevs=[992.556, 970.555, 929.649, 867.161, 787.702, 696.7955, 600.524, $ 510.455, 433.895, 368.818, 313.501, 266.481, 226.513, 192.5395, 163.6615, $ 139.115, 118.2502, 100.5145, 85.43898, 72.55784, 61.49566, 52.01594, $ 43.90968, 36.99269, 31.08889, 26.04909, 21.76096, 18.12431, 15.05023, $ 12.46014, 9.765464, 6.938725, 4.7273, 3.220675, 2.19422, 1.38914, $ 0.7320795, 0.339801, 0.1577215, 0.073208, 0.03398, 0.015772] endif if nohdr eq 72 then begin ; assume 2x2.5 model with 72 vertical levels (GEOS5) im=144 & jm=91 & km=nohdr lat=indgen(jm)*2-90 & lon=findgen(im)*2.5 plevs=[992.524, 977.521, 962.514, 947.524, 932.503, 917.481, 902.476, 887.511, 872.518, 857.475, $ 842.458, 827.492, 810.028, 787.517, 762.483, 737.473, 712.470, 687.477, 656.253, 618.735, $ 581.234, 543.764, 506.258, 468.728, 431.254, 393.763, 356.237, 312.781, 266.481, 226.513, $ 192.5395, 163.6615, 139.115, 118.2502, 100.5145, 85.43898, 72.55784, 61.49566, 52.01594, $ 43.90968, 36.99269, 31.08889, 26.04909, 21.76096, 18.12431, 15.05023, 12.46014, 10.2849, $ 8.45639, 6.91832, 5.63180, 4.56169, 3.67650, 2.94832, 2.35259, 1.86788, 1.47565, 1.15998, $ 0.907287, 0.705957, 0.546293, 0.420424, 0.321783, 0.244938, 0.185422, 0.139599, 0.104524, $ 0.0776725, 0.0567925, 0.0401425, 0.0263500, 0.0150000] endif ; Get lat, lon, and plevs if nohdr eq 0 then read_nc_met,filename,lat,lon,plevs km=n_elements(plevs) im=n_elements(lon) & jm=n_elements(lat) arr=fltarr(im,jm,km,rec) if keyword_set(trpp) then trpp=fltarr(im,jm,rec) ; Set type eq 'const' if keyword was not set if keyword_set(type) eq 0 then type='const' ; Always read in the true pressures tplevs=fltarr(im,jm,km,rec) if type eq 'const' or type eq 'kel' or type eq 'metwater' then psftype='psf' if type eq 'const_freq1' or type eq 'kel_freq1' or type eq 'potentialVorticity_freq1' $ or type eq 'metwater_freq1' then psftype='psf_freq1' if type eq 'const_overpass' or type eq 'kel_overpass' then psftype='psf_overpass' ; If we aren't reading from the first record in the file, offset by 'rec_init ; For all data types except hourly ozone if keyword_set(ddtype) eq 0 then ddtype='none' ; ---------Calculate day indices if reduce is set ----------------- if keyword_set(reduce) and ddtype ne 'amonthly' then begin if reduce ge 5 then begin int=fix(27/(reduce-1)) rx=indgen(reduce)*int endif ; if reduce lt 5 then begin if reduce eq 1 then rx=[0] if reduce eq 2 then rx=[0,14] if reduce eq 3 then rx=[0,10,20] if reduce eq 4 then rx=[0,7,14,21] endif endif ; ------------------------------------------------------ if ddtype ne 'column' then begin if keyword_set(rec_init) then m1=rec_init-1 else m1=0 for m=m1,m1+rec-1 do begin data=read_nc_3d(filename,m+1,im,jm,km,type,spnum=spec) ; We will only scale constituents to ppb. scal=strpos(type,'const') if type eq 'metwater' then scal=1 if scal ge 0 then arr(*,*,*,m-m1)=data*1e9 else arr(*,*,*,m-m1)=data ; Do not calculate true plevs if this is a monthly avg ; Change this!!! Allow true plevs for monthly averages if ddtype ne 'rst' and ddtype ne 'none' then begin read_nc_psfdata,filename,m+1,im,jm,km,psftype,levs tplevs(*,*,*,m-m1)=levs endif ; read_nc_psf can read in any 2D variable, such as psf or trpp if keyword_set(trpp) then begin if ddtype eq 'idaily' then data=read_nc_psf(filename,m+1,im,jm,'tropopausePress_freq1') if strmid(ddtype,0,4) eq 'over' then $ data=read_nc_psf(filename,m+1,im,jm,'tropopausePress_overpass') trpp(*,*,m-m1)=data endif endfor ; If keyword_set(reduce) and this is not 'amonthly', then reduce number of days returned if keyword_set(reduce) and ddtype ne 'amonthly' then begin arr=arr(*,*,*,rx) tplevs=tplevs(*,*,*,rx) if keyword_set(trpp) then trpp=trpp(*,*,rx) endif ; endif if ddtype eq 'column' then begin arr=fltarr(im,jm,rec) if keyword_set(rec_init) then m1=rec_init-1 else m1=0 for m=m1,m1+rec-1 do begin data=read_nc_2d_spec(filename,m+1,im,jm,spec,type) ; data=read_nc_2d(filename,m+1,spec,im,jm,type) arr(*,*,m-m1)=data endfor ; Reduce number of days if keyword_set(reduce) if keyword_set(reduce) then arr=arr(*,*,rx) endif ; ********* Optional: Put data on theta levels ********************** ; ******** Special treatment for theta levels <320K ***************** if keyword_set(thlevs) ne 0 then begin tmp=fltarr(im,jm,km,rec) if ddtype ne 'amonthly' then thplevs=(1000/tplevs)^0.288 if ddtype eq 'amonthly' then thplevs=(1000/plevs)^0.288 if type eq 'const' then type2='kel' if type eq 'metwater' then type2='kel' if type eq 'const_freq1' then type2='kel_freq1' if type eq 'const_freq2' then type2='kel_freq2' if type eq 'const_overpass' then type2='kel_overpass' tk=strpos(type,'kel') if tk ge 0 then type2=type if type eq 'potentialVorticity_freq1' then type2='kel_freq1' ; This reads in temperature for m=m1,m1+rec-1 do begin data=read_nc_3d(filename,m+1,im,jm,km,type2) tmp(*,*,*,m-m1)=data endfor ; **** Calculate theta values for arr & reduce number of days if keyword_set(reduce) **** if keyword_set(reduce) then begin tmp=tmp(*,*,*,rx) rec=n_elements(rx) theta=fltarr(im,jm,km,rec) endif ; Pressure levels are now a 4D field just like the constituent when it's not 'amonthly' if ddtype ne 'amonthly' then theta=tmp*thplevs if ddtype eq 'amonthly' then begin theta=tmp for k=0,km-1 do theta(*,*,k,*)=tmp(*,*,k,*)*thplevs(k) endif ; Interpolate constituent to the input thlevs ; This method consistently works for thlevs ge 320K, which is always there. Thlevs lt 320K are of ; interest at high latitudes but not extant globally. Examine each profile individually and ; interpolate where sensible. Do not use plevs at 800mb and above - avoid the lower troposphere. if thlevs(0) ge 320 then begin nk=n_elements(thlevs) arrth=fltarr(im,jm,nk,rec) for i=0,im-1 do for j=0,jm-1 do for m=0,rec-1 do $ arrth(i,j,*,m)=interpol(arr(i,j,*,m),theta(i,j,*,m),thlevs) arr=arrth endif ; if thlevs(0) lt 320 then begin a=where(thlevs lt 320) & na=n_elements(a) b=where(thlevs ge 320) & nb=n_elements(b) nk=na+nb ; First do interpolation for theta 320K and above arrth=fltarr(im,jm,nk,rec) for i=0,im-1 do for j=0,jm-1 do for m=0,rec-1 do $ arrth(i,j,b(0):nk-1,m)=interpol(arr(i,j,*,m),theta(i,j,*,m),thlevs(b(0):nk-1) ) ; Second do interpolation for levels below 320K, one profile at a time. Choose plevs 800mb -100 mb. ; Check for monotonicity of theta levels in this pressure range and check for minimum theta value ; in the profile. This is the min value for which you can interpolate. p=where(plevs lt 800 and plevs gt 100) np=n_elements(p) mono=fltarr(np-1) for i=0,im-1 do for j=0,jm-1 do for m=0,rec-1 do begin prof=theta(i,j,p,m) ; Check to see if the lowest theta of the profile is below the highest thlev that is under 320K. if min(prof) gt thlevs(a(na-1)) then goto, nextprof for n=0,np-2 do mono(n)=prof((n+1))-prof(n) if min(mono gt 0) then begin x=where(thlevs(a) gt min(prof)) arrprof=arr(i,j,p,m) arrth(i,j,x,m)=interpol(arrprof,prof,thlevs(x)) endif ; Completed one profile. Loop over profiles in 4D theta array. nextprof: endfor arr=arrth endif ; ****** FINAL: endif for optional thlevs ***************************** endif ; For thlevs lt 320K, there will be missing values (0). Don't average over them if zm=1 if keyword_set(zm) then arr=avg_good(arr,0,0) end