; Builds the global joint histogram for cloud optical thickness vs effective radius: ; for ice and liquid water. ; The MOD08_D3 dataset is used. The COSP modis function is copied and used here. ; The data is calcualted to monthly means for each year. load "myNclFunctions.ncl" idir = "/nobackup/rossby20/sm_marst/amipobs/modis/MOD08_D3/" years = tostring(ispan(2000,2017,1)) nyears = dimsizes(years) doys = ispan(1,366,1) ndoys = dimsizes(doys) ; cosp modis setup: intervals tauBounds = (/0.3, 1.3, 3.6, 9.4, 23., 60., 10000. /) ntau = 6 nreffliq = 6 ; Number of bins for tau/ReffLiq joint-histogram nreffice = 6 ; Number of bins for tau/ReffICE joint-histogram reffLIQ_Bounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/) reffICE_Bounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/) undef("GetVar") function GetVar(f:file,vname:string) local var begin print("Extracting variable = "+vname) var = short2flt(f->$vname$) var&YDim_mod08 = f->YDim var&XDim_mod08 = f->XDim var!0 = "lat" var!1 = "lon" copy_VarAtts(latatts,var&lat) copy_VarAtts(lonatts,var&lon) var = where(var.eq.-9999.0.or.var.lt.0.3,default_fillvalue("float"),var) var@_FillValue = default_fillvalue("float") if(vname.eq."Cloud_Effective_Radius_Ice_Mean" .or. \ vname.eq."Cloud_Effective_Radius_Liquid_Mean") then var = where(.not.ismissing(var),var*1.0e-6,var) ; convert to meters end if ;printVarInfo(var,"var") return(var) end undef("hist2d") ; histogram for joint modis products function hist2d(cot:numeric,reffliq:numeric,reffice:numeric) local dims,jointHist,tau,ref,csumliq,csumice,ice_hist,liq_hist begin rad = 4.0*atan(1.0)/180.0 clat = cos(cot&lat*rad) dims = dimsizes(cot) jointHist_liq = new((/dims(0),dims(1),ntau,nreffliq/),"float",0.0) jointHist_liq = 0.0 jointHist_ice = new((/dims(0),dims(1),ntau,nreffice/),"float",0.0) jointHist_ice = 0.0 liq_hist = new((/ntau,nreffliq/),"float",0.0) liq_hist = 0.0 ice_hist = new((/ntau,nreffice/),"float",0.0) ice_hist = 0.0 ;print("Calculating histogram...") print("Total valid points (liq) = "+num(.not.ismissing(cot).and..not.ismissing(reffliq))) csumliq = 0.0 do tau = 0,ntau-1 do ref = 0,nreffliq-1 ;print("Tau bins: "+tauBounds(tau)+" "+tauBounds(tau+1)) ;print("Reffliq: "+reffLIQ_Bounds(ref)+" "+reffLIQ_Bounds(ref+1)) jointHist_liq(:,:,tau,ref) = where(cot.ge.tauBounds(tau) .and. \ cot.lt.tauBounds(tau+1) .and. reffliq.ge.reffLIQ_Bounds(ref) .and.\ reffliq.lt.reffLIQ_Bounds(ref+1),1.0,0.0) csumliq = csumliq + sum(jointHist_liq(:,:,tau,ref)) ;print("npoints found (liq): "+ sum(jointHist_liq(:,:,tau,ref))) ; area weigted mean liq_hist(tau,ref) = wgt_areaave_Wrap(jointHist_liq(:,:,tau,ref),clat,1.0,1) end do end do print("Total valid points (ice) = "+num(.not.ismissing(cot).and..not.ismissing(reffice))) csumice = 0.0 do tau = 0,ntau-1 do ref = 0,nreffice-1 ;print("Tau bins: "+tauBounds(tau)+" "+tauBounds(tau+1)) ;print("Reffice: "+reffICE_Bounds(ref)+" "+reffICE_Bounds(ref+1)) jointHist_ice(:,:,tau,ref) = where(cot.ge.tauBounds(tau) .and. \ cot.lt.tauBounds(tau+1) .and. reffice.ge.reffICE_Bounds(ref) .and.\ reffice.lt.reffICE_Bounds(ref+1),1.0,0.0) csumice = csumice + num(jointHist_ice(:,:,tau,ref).eq.1.0) ;print("Csum = "+csumice) ;print("npoints found: "+ num(jointHist_ice(:,:,tau,ref).eq.1.0)) ice_hist(tau,ref) = wgt_areaave_Wrap(jointHist_ice(:,:,tau,ref),clat,1.0,1) end do end do ;print("Cumlative sum (liq) = "+csumliq) ;print("Cumlative sum (ice) = "+csumice) return([/liq_hist,ice_hist/]) end ;------ ; Main- ;------ begin First = True valid = new((/nyears,366/),"integer",0) theflnm = "modis_cot_vs_reffice_liq_jhist.nc" setfileoption("nc","format","netcdf4") setfileoption("nc","headerReserveSpace",64000) setfileoption("nc","preFill",False) if (isfilepresent(theflnm)) then system("rm " + theflnm) end if f = addfile(theflnm,"c") dim_names = (/"year","doy","tau","reffliq","reffice"/) dim_sizes = (/nyears,366,ntau+1,nreffliq,nreffice/) dim_unlimited = (/False,False,False,False,False/) filedimdef(f,dim_names,dim_sizes,dim_unlimited) atts = True atts@description = "Modis Jhist COT vs Reff Ice/Liq" atts@timestamp = systemfunc("date") print("Defining file attributes") fileattdef(f,atts) print("Defining file variables") filevardef(f,"year","float","year") filevardef(f,"doy","float","doy") filevardef(f,"tau","float","tau") filevardef(f,"reffice","float","reffice") filevardef(f,"reffliq","float","reffliq") filevardef(f,"cot_vs_reffice","float",(/"year","doy","tau","reffice"/)) filevardef(f,"cot_vs_reffliq","float",(/"year","doy","tau","reffliq"/)) f->tau = (/0.15,0.80,2.45,6.5,16.2,41.5,100.0/) f->reffliq = (/8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/) f->reffice = (/1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/) f->year = years f->doy = doys valid = 0 do y = 0,nyears - 1 do doy = 0,ndoys - 1 dd = sprinti("%0.3i",doys(doy)) cmd := idir+years(y)+"/"+dd ;print(cmd) chk = systemfunc("test -d "+cmd+"; echo $?") if(chk .eq. 0) then cmd := cmd+"/MOD08_D3*.hdf" ifile := systemfunc ("ls "+cmd) print(ifile) else ;print("Directory missing for => "+cmd) valid(y,doy) = 1 ; there is no data continue end if a := addfile(ifile,"r") cotice = GetVar(a,"Cloud_Optical_Thickness_Ice_Mean") ;printVarInfo(cotice,"cotice") ;Qplot(cotice) cotliq = GetVar(a,"Cloud_Optical_Thickness_Liquid_Mean") reffice = GetVar(a,"Cloud_Effective_Radius_Ice_Mean") reffliq = GetVar(a,"Cloud_Effective_Radius_Liquid_Mean") if(First) then dim = dimsizes(cotliq) cot_vs_reffice_jhist_awgt = new((/nyears,366,ntau+1,nreffice/),"float") cot_vs_reffice_jhist_awgt@_FillValue = default_fillvalue("float") cot_vs_reffice_jhist_awgt = default_fillvalue("float") cot_vs_reffliq_jhist_awgt = new((/nyears,366,ntau+1,nreffliq/),"float") cot_vs_reffliq_jhist_awgt@_FillValue = default_fillvalue("float") cot_vs_reffliq_jhist_awgt = default_fillvalue("float") First = False end if jhist = hist2d(cotliq,reffliq,reffice) cot_vs_reffliq_jhist_awgt(y,doy,1:7,:) = jhist[0] cot_vs_reffice_jhist_awgt(y,doy,1:7,:) = jhist[1] end do end do delete(f) end