#!/mnt/sportr1/raid1/cbblanke/miniforge3/envs/new_pycrtm/bin/python #was #!/usr/bin/env python3 #Adding: #At each location, write lat, lon, TIWV, CLW, SWS #How to get TIWV/PW? #PW(layer) = q * dp /g # ampr_retrieval.py 4/9/26-? # Clay Blankenship with ideas by Brent Roberts # Adapt binned_retrieval.py to work with AMPR. ############################################################ # binned_retrieval.py started on 12/6/2022 by Clay Blankenship # based on <1dvar-graphic_logsws_logcloud.py>. # That code did the retrieval including log(SWS), log(CLW), log(IWP) # This version will use different background, EOFs, and covariances for SST/Qcon bins. # The mean(background), EOFs, and covariances were calculated by # and are read from save5/evecs5_2003_2017.npz ################################################# #8/17/2026. Gets mostly good retrievals with 30% (guess) bad. #next step is to make a map of good/bad locations # from MerraProf import * import os, h5py, sys import numpy as np import datetime as dt from matplotlib import pyplot as plt import matplotlib.dates as mdates from pyCRTM import pyCRTM, profilesCreate,WATER_CLOUD,ICE_CLOUD from old_read_tbs import old_read_tbs import netCDF4 as nc from bilin import pseudoindex import pdb import time #for timing from global_land_mask import globe import cartopy.crs as ccrs import cartopy.feature as cfeature from pathlib import Path import pickle #from contextlib import redirect_stdout def merra_geti(lon): #find nearest neighbor index for merra #This works on a scalar or a numpy array #lon interval is 0.625 degrees, starting at -180 #ii=np.floor((lon+180.3125)/0.625).astype(int) #CBB fixed issue of i not looping back to 0 when over 179.6875 ii=np.floor(np.mod((lon+180.3125),360)/0.625).astype(int) if ii<0 or ii>=576: pdb.set_trace() return(ii) def merra_getj(lat): #find nearest neighbor index for merra. Grid is -90 to 90 by 0.5 degrees return np.floor(2*(lat+90.25)).astype(int) def getifloat(lon): #find nearest neighbor index for merra #This works on a scalar or a numpy array #lon interval is 0.625 degrees, starting at -180 #ii=np.floor((lon+180.3125)/0.625).astype(int) #returns integer #CBB fixed issue of i not looping back to 0 when over 179.6875 #ii=np.mod((lon+180.3125),360)/0.625 ii=np.mod((lon+180.),360)/0.625 if ii<0 or ii>=576: pdb.set_trace() return(ii) def getjfloat(lat): #find nearest neighbor index for merra. Grid is -90 to 90 by 0.5 degrees return 2*(lat+90.) def findindex(target,inarr): #Give the index of the element of inarr that is closest to target mydist=np.abs(np.array(inarr)-target) #mydist=np.abs(inarr-target) myindex=np.argmin(mydist) return(myindex) def read_bkstats(infile): #Read background profiles and covariance matrices for each of the [SST,qcon] bins #Additional info on clouds is in save4 mydata=np.load(infile) #Problem, missing this file. Need to regenerate it. #avgprof is calculated in units of T [K] and log Q [g/kg] #Binned quantities are by qcon bin (0-9) and then SST bin (0-8) #isst=int((SST-271.15)/4) #bin center SST[0 to 8]=[273,277....301,305] or [0,4....28,32] #jqcon=int(((qcon+2)/0.4)) #bin center qcon[0 to 9]=[-1.8, -1.4, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.4, 1.8] bin_avgprof=mydata['avgprof'] #(10,9,58) levels, T then w. The avg prof of the historical data. bin_count=mydata['count'] #(10,9) bin_evecs=mydata['evecs'] #(10,9,5,58) bin_evals=mydata['evals'] #(10,9,5) bin_sws=mydata['sws'] #(10,9) myplevs=mydata['plevs'] #(30): 10 to 1013 bin_bkcov=mydata['bkcov'] #(10,9,8,8) avgclwprof=mydata['clwprof'] #(10,9,29) avgiceprof=mydata['iceprof'] #(10,9,29) mydata.close() return(bin_avgprof, bin_count, bin_evecs, bin_evals, bin_sws, myplevs, bin_bkcov, avgclwprof, avgiceprof) def SetProfile(NLVL,pprof,tprof,qprof,clwprof,iceprof,oprof, mytsfc,SWS,yyyy,mm,dd,eia, ncloud=0): #profiles=SetProfile_merra( NLVL, pprof, tprof, np.exp(logqprof), clwprof, iceprof, oprof, tsfc[jj,ii], SWS, yyyy,mm,dd) nprof=1 #ncloud=0 #CBB debugging for n in range(0,nprof): #CBB test 4/13/26 trying to fix phase error #profiles = profilesCreate( nprof, NLVL-1, nClouds=2, nAerosols=0) #needed for CRTM, makes empty profile #profiles = profilesCreate( nprof, NLVL-1, nClouds=ncloud, nAerosols=0) #needed for CRTM, makes empty profile #print('call profilesCreate with ncloud=',ncloud) profiles = profilesCreate( nprof, NLVL-1, nClouds=ncloud, nAerosols=1) #needed for CRTM, makes empty profile profiles.Angles[n,0] = eia #CBB sensor zenith angle goes here profiles.Angles[n,1] = 999.9 #Sensor azimuth, not needed I guess? profiles.Angles[n,2] = 100.0 # 100 degrees, zenith below horizon. (CBB source zenith, i.e. sun) profiles.Angles[n,3] = 0.0 # zero solar azimuth #profiles.Angles[n,4] = h5['scanAngle'][()] profiles.Angles[n,4] = eia #CBB sensor scan angle goes here profiles.DateTimes[n,0] = int(yyyy) #year profiles.DateTimes[n,1] = int(mm) #month. CBB verify profiles.DateTimes[n,2] = int(dd) #day. CBB verify #profiles.O3[n,:] = np.sqrt(oprof[0:-1]*oprof[1:]) profiles.O3[n,:] = oprof[0:29] #CBB before 4/14/26, it was like this: # #profiles.aerosols[n,:,0,0] = np.asarray(h5['aerosolConcentration']) # profiles.aerosols[n,:,0,0] = 0.0 # #profiles.aerosols[n,:,0,1] = np.asarray(h5['aerosolEffectiveRadius']) # profiles.aerosols[n,:,0,1] = 3.0 # #profiles.aerosolType[n] = INVALID_AEROSOL # #profiles.aerosolType[n] = h5['aerosolType'][()] # profiles.aerosolType[n] = 1 # #print('aerosolType',profiles.aerosolType[n]) #profiles.aerosols[n,:,0,0] = np.asarray(profiles.aerosolConcentration) #profiles.aerosols[n,:,0,1] = np.asarray(profiles.aerosolEffectiveRadius) profiles.aerosols[n,:,0,0] = 0.0 profiles.aerosols[n,:,0,1] = 3.0 #profiles.aerosolType[n] = profiles.aerosolType profiles.aerosolType[n] = 1 #print('aerosolType',profiles.aerosolType[n]) #print('profiles.aerosols is ',profiles.aerosols) #print('profiles.clouds is ',profiles.clouds) profiles.climatology[n] = 6 #print('climatology',profiles.climatology) s=profiles.Salinity #print(type(profiles),type(profiles.Salinity),type(s)) profiles.Salinity[n] = 33.0 # just use salinity out of S2m for the moment. #print('salinity',profiles.Salinity) profiles.LAI[n] = 0.0 #CBB These are missing? #profiles.SurfType=[] #profiles.SurfGeom=[] # land, soil, veg, water, snow, ice profiles.surfaceTypes[n,0] = 1 #h5['landType'][()] (doesn't matter) profiles.surfaceTypes[n,1] = 1 #h5['soilType'][()] (doesn't matter) profiles.surfaceTypes[n,2] = 1 #h5['vegType'][()] (doesn't matter) profiles.surfaceTypes[n,3] = 1 #(1 for seawater, the only type.) h5['waterType'][()] profiles.surfaceTypes[n,4] = 1 #h5['snowType'][()] (doesn't matter) profiles.surfaceTypes[n,5] = 1 #h5['iceType'][()] (doesn't matter) profiles.surfaceFractions[n,:] = [0,1,0,0] # I think these are land/water/snow/ice?????? # CBB Will these vary? #profiles.windSpeed10m[n] = np.sqrt(u10m[jj,ii]*v10m[jj,ii]) #m/s #CBB update this profiles.windSpeed10m[n] = SWS profiles.windDirection10m[n] = 45.0 #deg E from N #profiles.surfaceTemperatures[n,:] = tsfc[jj, ii] profiles.surfaceTemperatures[n,:] = mytsfc #pprof is the layer pressures? profiles.Pi[n,:] = pprof profiles.P[n,:] = 0.5*(pprof[0:-1]+pprof[1:]) #should be "hypsometric avg"? #Put things that vary with profile here. #These will vary #profiles.T[n,:] = 0.5*(tprof[0:-1]+tprof[1:]) #profiles.Q[n,:] = np.sqrt(qprof[0:-1]*qprof[1:]) profiles.T[n,:] = tprof #profiles.Q[n,:] = np.exp(logqprof) #profiles.Q[n,:] = qprof*1000 #convert kg/kg to g/kg profiles.Q[n,:] = qprof #already in g/kg #CBB We have to set multiple cloud species (water,ice) #Input (MERRA) units for cloud: kg/kg (need air density) #output (CRTM) units for cloud: kg/m2 (int over layer, need layer thickness!) #This variable is (nprof,nlev) and does not iterate over cloud species if ncloud>0: profiles.cloudFraction[n,:] = np.ones(29) #CBB change if no cloud? #print('cloudType',profiles.cloudType[n]) delpa=100.0*(pprof[1:]-pprof[0:-1]) #print('setting clwprof*1000:',clwprof*1000) #profiles.cloudType[n,0]=np.asarray([WATER_CLOUD]) #Why this syntax? profiles.cloudType[n,0]=WATER_CLOUD profiles.clouds[n,:,0,0] = clwprof*delpa/9.81 profiles.clouds[n,:,0,1] = 10.0 #effective radius in microns #print('cloudFraction',profiles.cloudFraction[n,:]) #print('setting iceprof*1000:',iceprof*1000) #profiles.cloudType[n,1]=np.asarray([ICE_CLOUD]) #Why? if ncloud>1: profiles.cloudType[n,1]=ICE_CLOUD profiles.clouds[n,:,1,0] = iceprof*delpa/9.81 profiles.clouds[n,:,1,1] = 10.0 #print(clwprof*1000) # arg=input("Press enter, or type 9 to stop.") # if arg=='9': # pdb.set_trace() # return(profiles) def myformat(arr): fmt='%6.1f' outstring='' for i in range(0,len(arr)): x=fmt%arr[i] outstring=outstring+x return(outstring) def fdmodel(sensor_id,chan_set,lat,lon,pprof,tprof,logqprof, \ clwprof,iceprof,SST,SWS,eia): #CBB is profiles needed? If so, then why are tprof, logqprof... needed? np.set_printoptions(precision=2,suppress=True) #starts as True,True #If cloud goes negative, set that one to False #and set background value and variance appropriately. status=-1 #if (sensor_id == 'gmi_gpm'): # chan_set=[0,1,2,3,4,5,6,7,8] mycolors=['red','orange','gold','green','blue','lightblue','purple','lavender','pink','pink'] nchan=len(chan_set) #print('chan_set is ',chan_set) #print('nchan',nchan) qprof=np.exp(logqprof) print('clwprof:',clwprof) profiles=SetProfile( NLVL, pprof, tprof, qprof, clwprof, iceprof, oprof, SST, SWS, yyyy,mm,dd, eia) #CBB Check that Q is in right units: #print('profiles.q',profiles.Q[0,-1]) if profiles.Q[0,-1]>50 or profiles.Q[0,-1]<.1: print('Error in Q units or range') print(profiles.Q[0:]) pdb.set_trace() #TB calc for Merra crtmOb = pyCRTM() #Seg fault if omitted crtmOb.sensor_id = sensor_id crtmOb.profiles = profiles crtmOb.nThreads = 1 crtmOb.loadInst() #load the instrument, set number of channels crtmOb.channelSubset = chan_set+1 #CBB why? #print('Calling crtm runDirect') #print('channelSubset',crtmOb.channelSubset) crtmOb.runDirect() #CBB forward model for Merra #cals pycrtm.wrap_forward in pyCRTM.py->pycrtm.f90 #calctb = crtmOb.Bt[0,chan_set] calctb = crtmOb.Bt[0,:] if (verbosity>4): print('TBs are ',calctb) return(calctb) #def ret1dvar(sensor_id,chan_set,lat,lon,tbvec,xarr,tprof,logqprof,profiles,clwbk,icebk, \ # ff,merraprof,merratprof,merralogqprof,cloudstatus): #CBB 5/26/26 taking out cloudstatus for now (was last argument, a 2-length boolean) def ret1dvar(sensor_id,chan_set,lat,lon,tbvec,xarr,bkcov,evecs,evals,pprof,tprof,logqprof, \ clwbk,icebk,SST,SWS,ff,merraprof,merratprof,merralogqprof,eia,ncloud,tbmerra): doplot=False #CBB is profiles needed? If so, then why are tprof, logqprof... needed? np.set_printoptions(precision=2,suppress=True) #print('cloudstatus:',cloudstatus) print('tbvec',tbvec) #cloudstatus(first element: liquid, second: ice) is a boolean array of length 2. #starts as True,True #If cloud goes negative, set that one to False #and set background value and variance appropriately. yobs = tbvec #yobs = tbvec.data mdeltay=yobs-tbmerra print('Merra calc TBs are ',tbmerra) print('Delta-TB',mdeltay) mtbpen=np.sqrt(np.sum(mdeltay*mdeltay/tbpen_denom)/np.size(mdeltay)) #rmse #Needs to be normalized by TB error of each channel #xpen is already normalized by default? status=-1 #exsize=5 #for the spark plots exsize=6+ncloud #for the spark plots #CBB tring 6/2/22 #if (sensor_id == 'gmi_gpm'): # chan_set=[0,1,2,3,4,5,6,7,8] mycolors=['red','orange','gold','green','blue','lightblue','purple','lavender','pink','pink'] nchan=len(chan_set) exes=range(nchan) #Used to plot xarr exes2=range(0,exsize) #Set constant matrices rinv=np.diag(np.ones(nchan))*1.0 #inverse cov matrix of obs (TBs) bkcovinv=np.linalg.inv(bkcov) bsmallinv=np.eye(xarr.size) #diagonal matrix, size of #EOFs #CBB changing this 2/4/22 bsmallinv[5,5]=4.0 #equivalent to sigma=0.5 on the exponent if ncloud>0: bsmallinv[6,6]=4.0 #equivalent to sigma=0.5 on the exponent #Notes on this. value=1e6 (sigma-.001) ->haywire #val=0.25 (sigma=2) ->lots of leeway if ncloud>1: bsmallinv[7,7]=4.0 #Ice qprof=np.exp(logqprof) bkprof=np.append(tprof,logqprof) clwprof=clwbk iceprof=icebk print('clwprof:',clwprof) #profiles will consist of tprof, qprof, etc. initialized from background #merraprof (passed in already assembled) is the Merra version. print('call SetProfile for first guess, ncloud=',ncloud) profiles=SetProfile( NLVL, pprof, tprof, qprof, clwprof, iceprof, oprof, SST, SWS, yyyy,mm,dd, eia,ncloud=ncloud) print('cloud after SetProfile') print(profiles.clouds[0,:,0,0]*1000) #CBB Check that Q is in right units: #print('profiles.q',profiles.Q[0,-1]) if profiles.Q[0,-1]>50 or profiles.Q[0,-1]<.1: print('Error in Q units or range') print(profiles.Q[0:]) pdb.set_trace() #CBBFIX #if (cloudstatus[0]==False): # xarr[5]=0.0 # bsmallinv[5,5]=1.0e-6 #Is this right, set to a small value? #if (cloudstatus[1]==False): # xarr[6]=0.0 # bsmallinv[6,6]=1.0e-6 #Is this right, set to a small value? print('First guess profile:') #CBB Need to calculate xarr here if bkprof!=initial condition print('FG xarr is ',xarr) #print('pprof is ',pprof) print('tprof is ',tprof[[0,1,2,-3,-2,-1]]) print('ln qprof(g/kg) is ',logqprof[[0,1,2,-3,-2,-1]]) file2.write(f'initial xarr: {xarr} \n') print(f'initial xarr:',xarr,'\n') #CBB This is how to do it using the basis vectors origxarr=xarr.copy() oldxarr=xarr.copy() SWS=profiles.windSpeed10m[0] MSWS=merraprof.windSpeed10m[0] #The profile vector is bkprof plus xarr*EOF #So our first guess is the background if xarr=0 #CBB Initialize plot with original profile #CBB pprof has 30 nicely numbered levels (10,20...950,sfc) #the tprof and logqprof have 29 values on the in-between layers #Set up the axis from 0,29 and plot T,q on half-levels if doplot: fig,ax=plt.subplots(1,2,figsize=[12,6]) axwind=ax[0].inset_axes([0.05,0.1,0.3,0.25]) axwind.set_xlim(0,25.0) axwind.set_ylim(0,1.0) axwind.set_yticks([]) axwind.set_title('Surface wind (dots)') #axwind.set_xscale('log') #axwind.set_xlim(0.001,1.0) #axwind.set_xticks([0.001,0.01,0.1,1.0]) axwind.set_xticks([0,5,10,15,20,25]) axwind.plot(MSWS,0.55,marker='o',color='black') axcloud=ax[1].inset_axes([0.05,0.1,0.4,0.5]) axcloud.set_xlim(0.001,1.0) axcloud.set_ylim(0,1.0) axcloud.set_yticks([]) #axc=ax[1].secondary_xaxis('top') #axc=ax[1].twinx() axcloud.set_title('Cloud (dots)') axcloud.set_xscale('log') #axc.tick_params(axis='x',labelcolor='black') axcloud.set_xticks([0.001,0.01,0.1,1.0]) axcloud.plot(merraclwsum,0.24,marker='o',color='black') #for ice axcloud.plot(merraicesum,0.74,marker='o',color='black') #for clw yticks=np.array([0,2,4,6,8,10,12,14,18,22,25,27,29]) #indices of levels to put tick on yv=pprof[yticks] ystring=yticks.astype(str) for i in range(0,yticks.size): ystring[i]='{:,.0f}'.format(yv[i]) y=np.linspace(0.5,nlay-.5,nlay) #Start with Merra ax[0].plot(merratprof,y,linewidth=2,color='black') ax[0].text(0.58,0.95,'TBpen Xpen',transform=ax[0].transAxes) latlonst='{:,.0f}'.format(lat)+', '+'{:,.0f}'.format(lon) ax[0].text(0.1,0.1,latlonst,transform=ax[0].transAxes) ax[0].plot(merra_tsfc[jj,ii],nlay-1,marker='o',color='black') #ax[1].plot(logqprof,y,linewidth=2) ax[1].plot(np.exp(merralogqprof),y,linewidth=2,color='black') ax[1].set_xscale('log') ax[1].set_xlim(1.e-3,10) ax[0].set(xlabel='Temperature',ylabel='Pressure') ax[1].set(xlabel='Humidity') for n in range(0,2): ax[n].set_ylim(nlay-1,0.0) #ax[n].set_yticks(yv) ax[n].set_yticks(yticks) ax[n].set_yticklabels(ystring) plt.ion() plt.show() #We are going to retrieve xarr, then calculate prof from xarr #print('begin iteration with prof=',prof) #print('begin iteration with prof=',prof,file=ff) print() iter_tbpens=np.zeros(12) iter_xpens=np.zeros(12) tbpen=10000. di2=5000 #convergence criterion, initialize to high value xpen=0.0 lastpen=15000 #compare current totpen to lastpen stepfac=1.0 #TB calc for Merra(?) or should be FG? #### print('initialize CRTM') #### crtmOb = pyCRTM() #Seg fault if omitted #### #CBB Aerosol profile removed to avoid "phase errors" #### # Depends on presence or absence of aerosols and clouds #### #crtmOb.AerosolCoeff_File='AerosolCoeff.bin' #CBB TEST #### crtmOb.sensor_id = sensor_id #### crtmOb.profiles = merraprof #### crtmOb.nThreads = 1 #### print('crtmOb.loadInst()................................') #### crtmOb.loadInst() #load the instrument, set number of channels #### print('Total channels loaded:',crtmOb.nChanTotal) #### crtmOb.channelSubset = chan_set+1 #### print('call runDirect................................') #### crtmOb.runDirect() #CBB This is what runs the Forward Model #### #merratb = crtmOb.Bt[0,chan_set] #### tbcalc = crtmOb.Bt[0,:] #was merratb (WRONG name), this is from first guess #### print('tbvec/yobs:',yobs) #### print('tbcalc0 :',tbcalc) #tbvec/yobs: [148.78 95.83 95.34 149.68 # 165.89 113.05 112.71 166.88 # 188.83 137.79 137.53 189.78 # 226.97 181.38 181.58 227.87] #tbcalc : [149.67 95.29 166.49 115.01 187.59 137.73 228.33 186.15] #### deltay=yobs-tbcalc #### print('TBs are ',tbcalc) #### print('Delta-TB',deltay) #### #Is this needed yet? #### #tbpen=np.sqrt(np.sum(deltay*deltay)/np.size(deltay)) #rmse ## if y==ny-10: ## pdb.set_trace() #Decompose Merra profile into eofs here for comparison mprof=np.append(merratprof,merralogqprof) xmerra=np.matmul(evecs,mprof-bkprof)/evals #CBB Augment the xmerra vector (eof coeffs) with the cloud coeffs xmerra=np.append(xmerra,np.log(profiles.windSpeed10m[0])) if ncloud>0: xmerra=np.append(xmerra,xarr[6]) #CBB or call this xmerra_new? if ncloud>1: xmerra=np.append(xmerra,xarr[7]) #CBB this is ice, may be off mxpen=np.sqrt(np.sum(xmerra*xmerra)/np.size(xmerra)) #rmse file2.write(f'xmerra: {xmerra} \n') if doplot: ax[1].text(0.50,0.95,'Merra TBpen/Xpen',transform=ax[1].transAxes) ax[1].text(0.58,0.90,'{:5.1f}'.format(mtbpen),transform=ax[1].transAxes,color='black') ax[1].text(0.73,0.90,'{:5.1f}'.format(mxpen),transform=ax[1].transAxes,color='black') mystring='SST, Qcon: '+"{:.2f}".format(mysst)+" {:.2f}".format(myqcon) ax[1].text(0.55,0.85,mystring,transform=ax[1].transAxes) axtiny=ax[1].inset_axes([0.38,0.90,0.15,0.05]) axtiny.bar(exes,deltay) axtiny.set_ylim(-10,10) axtiny.axis('off') axtiny2=ax[1].inset_axes([0.83,0.90,0.15,0.05]) axtiny2.bar(exes2,xmerra) axtiny2.set_ylim(-3,3) axtiny2.axis('off') #A record of cloud and SWS clwlog=np.zeros(10) icelog=np.zeros(10) swslog=np.zeros(10) #Now the retrieval chisq=9999. for it in range(0,10): #Retrieval iterations. Change to a while loop? print () #print ('n=',it,':') crtmOb = pyCRTM() #Seg fault if omitted crtmOb.sensor_id = sensor_id crtmOb.profiles = profiles #This creates a "pointer" (OK?) #print('cloud at start of retrieval') #Verified same as "after SetProfile" #print(profiles.clouds[0,:,0,0]*1000) crtmOb.nThreads = 1 crtmOb.loadInst() #load the instrument, set number of channels #print('call pyCRTM') #CBB May have to define crtmOb.CloudCoeff_File beofore crtmOb.loadInst() *** #CBB can define crtmOb.channelSubset here # *** CBB temporary fd. TB calc to verify retrieval #profiles.T[n,:] = prof[0:29] #profiles.Q[n,:] = np.exp(prof[29:58]) #CBB Initial TB calculation before loop. #CBB Also used for testing to set the fake TBobs by adding a small perturbation. #sys.stdout=open('stdout.txt','w') #with open('stdout.txt','w') as ff: # with ff as sys.stdout: #with redirect_stdout(ff): crtmOb.channelSubset = chan_set+1 #starting with one vs. zero crtmOb.runDirect() #CBB This is what runs the Forward Model #sys.stdout.close() #Tbcalc = np.copy(crtmOb.Bt) #Tbcalc = crtmOb.Bt[0,chan_set] Tbcalc = crtmOb.Bt[0,:] if it==0: tbinit=Tbcalc if(np.isnan(np.sum(Tbcalc))): #check for any nan status=-3 #blew up print('tbcalc nan!') print() tiwv=-1 clwsum=-1 icesum=-1 break #exit loop #print('TBs: ',Tbcalc) #print('TBs are ',Tbcalc,file=ff) #deltay=np.squeeze(yobs-Tbcalc[0,:]) deltay=np.squeeze(yobs-Tbcalc) #print('deltaTb',deltay) prev_tbpen=tbpen #CBB changing to include normalization by nedt (guess 1.5k) #tbpen=np.sqrt(np.sum(deltay*deltay)/np.size(deltay)) #rmse tbpen=np.sqrt(np.sum(deltay*deltay/tbpen_denom)/np.size(deltay)) #rmse iter_tbpens[it]=tbpen if(it==0): origtbpen=tbpen xdiff=xarr-origxarr #CBB changing on 9/11/26 #xpen=np.sum(xdiff*xdiff)/np.size(xdiff) #Is this right time to calculate? xpen=np.sqrt(np.sum(xdiff*xdiff)/np.size(xdiff)) #Is this right time to calculate? iter_xpens[it]=xpen #print('xpen [xarr]',xpen,xarr) #print(myformat(deltay)+'\n',file=ff) #print('TB penalty:','%8.1f'%tbpen) #print('deltaTb',deltay,'TBpen:','%6.1f'%tbpen) #print('xpen,tbpen ', '%5.1f'%xpen,'%6.1f'%tbpen) totpen=xpen+tbpen if totpen>lastpen: stepfac=stepfac*0.5 lastpen=totpen if doplot: ax[0].text(0.58,0.90-0.05*it,'{:5.1f}'.format(tbpen),transform=ax[0].transAxes,color=mycolors[it]) ax[0].text(0.73,0.90-0.05*it,'{:5.1f}'.format(xpen),transform=ax[0].transAxes,color=mycolors[it]) ax[0].plot(tprof,y,color=mycolors[it]) ax[1].plot(np.exp(logqprof),y,color=mycolors[it]) if ncloud>0: clwsum=np.sum(profiles.clouds[0,:,0,0]) #should be mm if doplot: axcloud.plot(clwsum,0.2,marker='o',color=mycolors[it]) clwlog[it]=clwsum print('clw','%6.2f'%clwsum) if ncloud>1: icesum=np.sum(profiles.clouds[0,:,1,0]) if doplot: axcloud.plot(icesum,0.7,marker='o',color=mycolors[it]) icelog[it]=icesum print('ice','%6.2f'%icesum) swslog[it]=SWS #Not sure if this needs to be recalculated but print it for diagnostics delpa=100.0*(profiles.Pi[0,1:]-profiles.Pi[0,0:-1]) #29 levels, in Pa tiwv=0.001*np.dot(delpa,profiles.Q[0,:])/9.81 #in mm or kg/m2 #CBB fixed 9/14/26 print(f'TIWV: {tiwv:.1f} SWS: {SWS:.1f} ') if doplot: axwind.plot(np.exp(xarr[5]),0.5,marker='o',color=mycolors[it]) #CBB Plot a tiny bar plot of the delta-t axtiny=ax[0].inset_axes([0.38,0.90-0.05*it,0.15,0.05]) axtiny.bar(exes,deltay) axtiny.set_ylim(-10,10) axtiny.axis('off') #CBB Plot a tiny bar plot of the Xarr (EOF coeffiecients) axtiny2=ax[0].inset_axes([0.83,0.90-0.05*it,0.15,0.05]) axtiny2.bar(exes2,xarr) axtiny2.set_ylim(-3,3) axtiny2.axis('off') #Convergence test. #A bit confusing in Rodgers but di2 is calculated by taking the RHS of 5.31 (prod3+prod4) #and multiplying by prod5=deltax=xarr-oldxarr (5.29,30,31) #This quantity should be <0): di2=np.matmul(prod5,prod4+prod3) if (verbosity>4): print(' di2',di2) #if (di22*origtbpen): # status=-2 #increasing # print('No retrieval, TBpen increasing') # break #CBB: Now do the Jacobian calculation forwardEmissivity = crtmOb.surfEmisRefl[0,:] #I don't think we need this. crtmOb.surfEmisRefl = [] #with open('stdout.txt','w') as ff: #with redirect_stdout(ff): #CBB might need next line crtmOb.output_cloud_K=True #CBB attempting to get Jacobian for cloud crtmOb.runK() #produces crtmOb.TK and .QK #crtmOb.CloudConcentrationK has dimensions (nchan=8, 1, nlev=29, 1) prof=np.append(tprof,logqprof) #move? oldprof=np.copy(prof) #CBB: Solve for new profile based on delta-Y and Kmat ############################### #Construct kTB from crtmOb.TK and .QK #TK and QK are (1,22,29) #22 channels, 42 levels #kTb = np.concatenate((crtmOb.TK[0,chan_set,:],crtmOb.QK[0,chan_set,:]),axis=1) #not used #QK=crtmOb.QK[0,chan_set,:] #This was calculated for g/kg QK=crtmOb.QK[0,:,:] #This was calculated for g/kg qprof=np.exp(logqprof) #### CBB need to add extra factor of Q to account for Q being in log space #CBB: When we have the cloudK, we will need to set Kx[1-5] this way and then set the next two columns #from K(clw) and K(ice) #Kx = np.concatenate((crtmOb.TK[0,chan_set,:],1000*QK*qprof),axis=1) #CBB for the 5 EOFs #Kx = np.concatenate((crtmOb.TK[0,chan_set,:],QK*qprof),axis=1) #CBB for the 5 EOFs Kx = np.concatenate((crtmOb.TK[0,:,:],QK*qprof),axis=1) #CBB for the 5 EOFs Kwind=crtmOb.WindSpeedK #if cloudopt=='Yes': #if cloudstatus[0]: if ncloud>0: #Kclw=crtmOb.cloudK[0,chan_set,:,0] #profile,channel,species #Kclw=crtmOb.cloudK[0,:,:,0] #profile,channel,species Kclw=crtmOb.CloudConcentrationK[:,0,:,0] #dims are nchan, ?, nlev, ? #print('Kclw',Kclw[0,:]) else: Kclw=np.zeros([nchan,29]) #if cloudstatus[1]: if ncloud>1: Kclw=crtmOb.CloudConcentrationK[:,0,:,1] #dims are nchan, ?, nlev, species(?) else: Kice=np.zeros([nchan,29]) #CBB Dimensions are not right. Kmat should end up (nchan? x ncoeffs) #where coeffs = size of state vector (5 EOFs plus 2 cloud profiles) #think of the cloud profiles as 2 extra eofs kmat5=np.matmul(Kx,np.transpose(evecs)) #Jacobian wrt. eigenvectors [nchan,58][58,5]=[9,5] # Kmat test was here (moved to comments at end) #CBB Verified that Kx and kmat are accurate by forward differencing #kmat tells the delta-TB for a given change in the principal component vector #We need to augment kmat5 with the changes for SWS and the 2 cloud vectors kmatsws=Kwind*SWS kmat=np.append(kmat5,np.transpose(kmatsws),axis=1) #kmatclw is dTb/dxarr[5] #CBB the following conversion verified by forward differencing kmatclw=np.matmul(Kclw, profiles.clouds[0,:,0,0] ) #should come out to [nchan,29][29,1]=[nchan,1] kmatclw=np.expand_dims(kmatclw,axis=1) #convert to (nchan,1) kmat=np.append(kmat,kmatclw,axis=1) if ncloud>1: kmatice=np.matmul(Kice, profiles.clouds[0,:,1,0] ) #should come out to [nchan,29][29,1]=[nchan,1] kmatice=np.expand_dims(kmatice,axis=1) #convert to (nchan,1) kmat=np.append(kmat,kmatice,axis=1) #now (nchan,8) kt_rinv=np.matmul(np.transpose(kmat),rinv) #(8,nchan)(nchan,nchan)=(8,nchan) kt_rinv_k=np.matmul(kt_rinv,kmat) #(8,nchan)(nchan,8)=(8,8) #This is the x(n+1)=x(0)+.... #CBB: This was the "linear" way. (2/18/22) #CBB prod2=np.matmul(np.linalg.inv(bsmallinv+kt_rinv_k),kt_rinv) #CBB prod3=np.matmul(prod2,deltay) #CBB prod3.shape #CBB oldxarr=np.copy(xarr) #CBB xarr=oldxarr+prod3 #CBB ^This block is based on oldxarr (Rodgers 5.8) but is missing the #CBB attempting new solution with nonlinear term try: sx=np.linalg.inv(bsmallinv+kt_rinv_k) except: print('Matrix inversion error') ff.write('Matrix inversion error\n') status=-2 tiwv=-1 clwsum=-1 icesum=-1 break #prod1 is kt_rinv prod3=np.matmul(bsmallinv,(origxarr-xarr)) #sign OK? #CBB initially zero prod4=np.matmul(kt_rinv,deltay) #make sure deltay is updated prod5=np.matmul(sx,prod4+prod3) oldxarr=np.copy(xarr) #CBB adding stepfac 8/31 #CBB incrementing previous xarr is correct (Rodgers 5.8) xarr=oldxarr+prod5*stepfac #This updates the xarr xdiff=xarr-origxarr xpen=np.sqrt(np.sum(xdiff*xdiff)/np.size(xdiff)) #rmse file2.write(f' xarr: {xarr} \n') #if (it==0): # print('xarr xpen tbpen totpen di2') #print(xarr,' ',np.round(xpen,1),np.round(tbpen,1),np.round(xpen+tbpen,1),np.round(di2,2)) #print('it ',it,':',' xarr: ',xarr,' xpen: ',xpen, 'tbpen:' ,tbpen, 'di2: ',di2) print(f'it: {it} xarr: {xarr} ') print(f'xpen: {xpen:.2f} tbpen {tbpen:.2f} di2: {di2:.2f}') #print(' xarr:',xarr) #CBB adding additional check for extreme values in x vector if any(abs(xarr[0:5])>10): status=-5 #blew up in xarr check #pdb.set_trace() #CBB construct tprof,logqprof from xarr and eofs oldprof=np.copy(prof) #Maybe useful? #prof=bkprof+np.matmul(xarr,evecs) #This updates the bkprof (which ->T,logq) #CBB Updated to handle the cloud coeffs in the xarr prof=bkprof+np.matmul(xarr[0:5],evecs) #This updates the bkprof (which ->T,logq) #wrong prof=prof+np.matmul(xarr,evecs) #print('xarr is ',xarr) #CBB not needed since using log #if (xarr[5]<0): # #CLW tried to go negative, turn off liquid cloud # print('xarr[5] is ',xarr[5],'... Setting CLW to zero ******') #need to restar? # cloudstatus[0]=False # xarr[5]=0.0 # origxrr[5]=0.0 # bsmallinv[5,5]=1.0e-6 #Is this right, set to a small value? #if (xarr[6]<0): # #CLW tried to go negative, turn off liquid cloud # cloudstatus[1]=False # print('xarr[6] is ',xarr[6],'... Setting ICE to zero ******') # xarr[6]=0.0 # origxarr[6]=0.0 # bsmallinv[6,6]=1.0e-6 #Is this right, set to a small value? SWS=np.exp(xarr[5]) if ncloud>0: clwprof=np.exp(xarr[6])*clwbk #CBB these will be passed to SetProfile if ncloud>1: iceprof=np.exp(xarr[7])*icebk tprof=prof[0:nlay] logqprof=prof[nlay:2*nlay] #units remain kg/kg if (verbosity>4): print(' tprof is ',tprof[[0,1,2,-3,-2,-1]]) print(' ln qprof(g/kg) is ',logqprof[[0,1,2,-3,-2,-1]]) #CBB Code to print out current prof #print('after ',i+1,' iterations, prof is ...') #print(prof) #print('after ',i+1,' iterations, prof is ...',file=ff) #print(prof,file=ff) #Check for out of bounds to stop errors in CRTM if any(tprof<100) or any(tprof>399) or np.isnan(np.sum(prof)): print('TBpen history:',iter_tbpens[0:it]) print('Xpen history: ',iter_xpens[0:it]) #pdb.set_trace() #Example of "blow up" # TBpen 95, 134, 157 # Xpen 0, 58, 738... # so tbpen was already large iter_tbpens[it]=9999. status=-3 #blew up print('*** T profile out of bounds!',min(tprof),max(tprof)) #print('profile out of bounds!',min(tprof),max(tprof),file=ff) gg.write('profile blew up\n') #gg.write(f'{tprof: *tprof:.1f}\n') #gg.write('{:5.1f}'.format(tprof)+'\n') #gg.write(tprof) gg.write(myformat(tprof)+'\n') #CBB Save info so we can try this one again. #print(date,hindex,m,file=gg) #date not defined tiwv=-1 clwsum=-1 icesum=-1 break #CBB Some of these variables are known because they are global #i.e. defined in the top level. #print('SetProfile in ret1dvar loop, ncloud=',ncloud) profiles=SetProfile(NLVL,pprof,tprof,np.exp(logqprof),clwprof,iceprof,oprof,tsfc[jj,ii], SWS,yyyy,mm,dd,eia, ncloud=ncloud) #print('100*liquid cloud:',100*merraprof.clouds[0,:,0,0]) #print('100* ice cloud:',100*merraprof.clouds[0,:,1,0]) #crtmOb = pyCRTM() #Seg fault if omitted #crtmOb.sensor_id = sensor_id #crtmOb.profiles=profiles #if (i==2): # pdb.set_trace() #Now that new profile is calculated, repeat loop. #print('cloud after iteration ',it) ##print('cloud: ',profiles.clouds[0,:,0,0]*1000) if status>=0: if (verbosity>3): print('end of iterations',it) print(' tprof is ',tprof[[0,1,2,-3,-2,-1]]) print(' ln qprof(g/kg) is ',logqprof[[0,1,2,-3,-2,-1]]) #print('tprof is ',tprof) #print('ln qprof(g/kg) is ',logqprof) #print('tbpen',tbpen,' xpen',xpen,' total',tbpen+xpen) #penstring=' '+str(["{:8.2f}".format(x) for x in tbpens]) #outstring="%.2f"%ptlat+" %.2f"%ptlon+" %3i"%status+" %2i"%i+" %.2f"%origtbpen+penstring+"\n" #ff.write(outstring) #output file print('clw history',clwlog[0:it+1],' merra=',np.round(merraclwsum,2)) #print('ice history',icelog[0:it+1],' merra=',np.round(merraicesum,2)) print('sws history',swslog[0:it+1],' merra=',np.round(MSWS,2)) #CBB Put a chi-squared test here, Rodgers 5.32 (and 5.27) se=np.linalg.inv(rinv) bsmall=np.linalg.inv(bsmallinv) sdeltay=calc_sdeltay(se,kmat,bsmall) prod1=np.matmul(deltay,sdeltay) chisq=np.matmul(prod1,deltay) print('chisq(deltay)',np.round(chisq,2)) #(5.27), this should follow chi-sq distribution print(yyyy,mm,dd,hh,'%8.2f'%di2,'%8.2f'%tbpen,'%8.2f'%xpen,'%8.2f'%(tbpen+xpen),'%8.2f'%chisq, \ '%4i'%status,file=ff) #print('xmerra',xmerra) #CBB assign status if none reached yet if (tbpen+xpen < nchan+exsize): #average error <1 status=1 #good else: status=-4 #too many iterations #pdb.set_trace() #print() #if (di2<2): # status=0 #good #Do all this and write to file3 regardless of status clwsum=0.0 icesum=0.0 if ncloud>0: clwsum=np.sum(profiles.clouds[0,:,0,0]) #should be in mm icesum=0.0 if ncloud>1: icesum=np.sum(profiles.clouds[0,:,1,0]) #Get TIWV #NLVL is 30, Pi(30),P(29),q(29) #note: delpa=100.0*(pprof[1:]-pprof[0:-1]) delpa=100.0*(profiles.Pi[0,1:]-profiles.Pi[0,0:-1]) #29 levels, in Pa tiwv=0.001*np.dot(delpa,profiles.Q[0,:])/9.81 #in mm or kg/m2 #CBB fixed 9/14/2026 #print('clwsum,icesum',clwsum,icesum) #print('merra ',merraclwsum,merraicesum) print() merra_tiwv=0.001*np.dot(delpa,merraprof.Q[0,:])/9.81 merra_clwsum=np.sum(merraprof.clouds[0,:,0,0]) #MSWS is Merra wind speed #CBB added EIA file3.write(f"{status:3d} {ptlat:9.2f} {ptlon:9.2f} {eia:5.2f} {scanpos:3d} {tiwv:9.3f} {1000*clwsum:9.3f} {SWS:5.2f}\n") file3.write(f"Merra: {merra_tiwv:9.3f} {1000*merra_clwsum:9.3f} {MSWS:5.2f} {mtbpen:9.2f}\n") file3.write('TBpen:'+" ".join(f"{x:9.2f}" for x in iter_tbpens[0:it])+'\n') file3.write('Xpen: '+" ".join(f"{x:9.2f}" for x in iter_xpens[0:it])+'\n') file3.write('tbinit:'+" ".join(f"{x:9.2f}" for x in tbinit)+'\n') file3.write('tbvec:'+" ".join(f"{x:9.2f}" for x in tbvec)+'\n') file3.write('tbmerra:'+" ".join(f"{x:9.2f}" for x in tbmerra)+'\n') #file3.write('') print('status:',status ) #if tbpen>10: # pdb.set_trace() #dum=input('enter to continue, 9 to stop') #if dum=='9': # pdb.set_trace() return(status,prof,it,tbpen,xpen,mtbpen,mxpen,chisq,tiwv,clwsum,icesum,SWS,tbinit,merra_tiwv) #return(status,pprof,it) #ret1dvar def calc_sdeltay(se,kmat,sa): #Rodgers eq 5.27 prod1=np.matmul(kmat,sa) prod2=np.matmul(prod1,np.transpose(kmat)) sum1=prod2+se try: sdeltay=np.matmul(np.linalg.inv(sum1),se) except: print('maxtrix inversion error--singular matrix') print(sum1) pdb.set_trace() print() return(sdeltay) def profplot(xvec,plevs): #Set up a plot with pressure on a log scale from matplotlib import ticker fig=plt.figure() ax=plt.axes() ax.set_ylabel("Pressure") ax.set_yscale('log') ax.set_ylim(1000,100) yvec=plevs ax.plot(xvec,plevs) subs=[1,2,3,4,5,6,7,8,9] loc=ticker.LogLocator(base=10.,subs=subs) ax.yaxis.set_major_locator(loc) fmt=ticker.FormatStrFormatter("%g") ax.yaxis.set_major_formatter(fmt) ylabels=ax.get_yticklabels() return(fig,ax) def get_vertical_profile(data, x_coords, y_coords, target_x, target_y): """ Extracts a vertical profile from a 3D field (Z, Y, X). Args: data: 3D numpy array shaped (Z, Y, X). x_coords: 1D array of X coordinates (longitude/easting). y_coords: 1D array of Y coordinates (latitude/northing). target_x: The X coordinate for the profile. target_y: The Y coordinate for the profile. Returns: 1D array containing the interpolated temperature values for each vertical level. """ from scipy.interpolate import RegularGridInterpolator # Create the horizontal interpolator using the full 3D grid # We treat Z as a fixed dimension and interpolate over (Y, X) z_levels = data.shape[0] profile = np.zeros(z_levels) # Iterate through each vertical level to perform 2D bilinear interpolation for i in range(z_levels): # RegularGridInterpolator expects coordinates in the order of the data axes (Y, X) interp = RegularGridInterpolator((y_coords, x_coords), np.array(data[i, :, :])) profile[i] = interp([target_y, target_x])[0] #mysum=np.sum(interp) #if(np.isnan(mysum)): #check for any nan # print(profile[i]) # pdb.set_trace() # print() return profile def fix_tprof(tprof,pprofup): #Fill in nans below the surface #tprof is 29 points, bottom up (no t2m level) #pprofup is 29 points, bottom up (no psfc) return() def scanplot(tb,title,cmap='viridis',vmin=75,vmax=325): #copied from tbplot in fd_model_whymsie_debug #create a rectangular plot using (along track, across track) coordinates time_formatter=mdates.DateFormatter('%H:%M:%S') fig,ax=plt.subplots(figsize=(10,8)) #SDTV try: im=ax.imshow(tb, cmap=cmap,interpolation='nearest',extent=extent, aspect='auto',vmin=vmin,vmax=vmax) except: #This happens if cmap is not valid... pdb.set_trace() print() plt.title(title) plt.xlabel('time') plt.ylabel('scan position') #ax.xaxis.set_major_formatter(mticker.FuncFormatter(time_formatter_fn)) ax.xaxis.set_major_formatter(time_formatter) ax.set_xlabel("Time") #ax.tick_params(axis='x', labelrotation=-45) fig.colorbar(im,ax=ax,label='TB [K]',orientation='vertical',cmap=cmap) plt.tight_layout() return def myplot(lats,lons,var,title,bartitle,vmin,vmax,extent): # 3. Create the Map Plot fig = plt.figure(figsize=(10, 8)) # Use the PlateCarree projection (standard lat/lon) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # Add geographic features for context ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5) ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.5) ax.add_feature(cfeature.COASTLINE, linewidth=1) ax.add_feature(cfeature.STATES, linestyle=':') # Create a scatter plot of the data # We use c=statuses to color-code by the status variable print('##',len(lons),len(var) ) scatter = ax.scatter(lons, lats, c=var, cmap='viridis', transform=ccrs.PlateCarree(), s=3, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) # Add a colorbar to interpret the status colors plt.colorbar(scatter, label=bartitle, shrink=0.7) ax.set_extent(extent) # Add gridlines and coordinate labels gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl.top_labels = False gl.right_labels = False plt.title(title) plt.ion() plt.show() return(fig) if __name__ == "__main__": np.set_printoptions(precision=2,suppress=True) verbosity=1 #higher gives more messages bkopt='merra' fgopt='merra' cloudopt='Yes' #CBB change this to cloudskip option # bkopt='avg' (default): Use avgprof as the background profile # bkopt='merra' : Use Merra profile as the background profile # fgopt='avg' (default): Use avgprof as the first guess profile # fgopt='merra' : Use Merra profile as the first guess profile #cloudopt #'Yes' or 'No'...set 'Yes' to retrieve cloud #set 'No' to not retrieve clouds and skip cloudy (Merra) profiles #outdir='out/20241023_sec' outdir='out/20241023_new' mypath=Path(outdir) mypath.mkdir(exist_ok=True) file1 = open(outdir+"/log.txt", "w", encoding="utf-8") #summary line of each retrieval file2 = open(outdir+"/steps.txt", "w", encoding="utf-8") #more detailed info on each retrieval file3 = open(outdir+"/penalty_steps.txt", "w", encoding="utf-8") #more detailed info on each retrieval file4 = outdir+"/savefile" tbpen_denom=2.25 #(square of obs+fd model error) #Load bias corrections #Do I have these? bias_corr=False if bias_corr: bcdata=np.load(biasfile) bin_medians=bcdata['binmedians'] bin_bias=bcdata['binbias'] all_medians=bcdata['allmedians'] all_bias=bcdata['allbias'] #infile='save5plus/evecs5_2003_2017.npz' #this has vapor as log q infile='backgrounds/evecs5_201201_201512.npz' #vapor as log q? bin_avgprof, bin_count, bin_evecs, bin_evals, bin_sws, myplevs, bin_bkcov, avgclwprof, avgiceprof = \ read_bkstats(infile) #Set up retrieval ngood=0 nbad=0 nskip=0 ncloud=1 #sensor_id = 'gmi_gpm' #sensor='GMI' sensor_id = 'ampr_air' sensor='AMPR' chan_set=np.array([0,1,2,3,4,5,6,7]) #AMPR ff=open('test_1dvar.out','w',buffering=1) #output file ff.write(' di2 tbpen xpen total chisq') #ff.write('lat,lon,stat,n_it,oribtbpen,tbpen\n') #output file gg=open(sensor_id+'_details.txt','w',buffering=1) #File to write detailed results #gg.write('lat,lon,stat,n_it,oribtbpen,tbpen\n') #output file #Load the Olympex file (or other flight data) #fileName='OLYMPEX_20151123_0_1948_AMPR-CoSMIR_merger_conical_v3.nc' #fileName='OLYMPEX_20151123_0_1948_AMPR-CoSMIR_merger_conical_v3.nc' #fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241018_R1.nc' #fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241022_R1.nc' #Oklahoma? #fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241031_R1.nc' fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241023_R1.nc' #west of Mexico, mostly clear #fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241025_R1.nc' #over land #fileName='AMPR/AMPR_WH2yMSIE/whymsie-AMPR_ER2_20241030_R1.nc' #over land? ds=nc.Dataset(fileName) #print(ds.variables) nang=ds.dimensions['CrossTrackDim'].size #50 nscan=ds.dimensions['AlongTrackDim'].size #nx=np.size(ds.dimensions['AlongTrackDim']) #ny=np.size(ds.dimensions['CrossTrackDim']) nx=nscan ny=nang print('nscan=',nscan) frequencies=ds.variables['Frequency'] polarization=ds.variables['Channel'] #A,B,H,V scan_angles=ds.variables['ScanAngle'] ampr_lats=ds.variables['Lat'] #see also GPSLat ampr_lons=ds.variables['Lon'] #cosmir_tb=ds.variables['cosmir_tb'] #(9,1153,51) #cosmir_IA=ds.variables['cosmir_IA'] #(1153,51) #incidence angles #cosmir_lat=ds.variables['cosmir_lat'] #(1153,51) #cosmir_lon=ds.variables['cosmir_lon'] #(1153,51) ampr_tb=ds.variables['TB'] #(nchan,npol,nscan,50) This is V and H #dims are (4 pol, 4 freq, 4377, 50) #channel order is 10.7,19.35,37.1,85.5 ampr_tb_v=ampr_tb[3,:,:,:] ampr_tb_h=ampr_tb[2,:,:,:] ## #quick plot of flight path (show tbobs) ## #extent=[-125,-110,21,39] ## extent=[-119,-117,25,27] ## lats=np.array(ampr_lats).flatten() ## lons=np.array(ampr_lons).flatten() ## #tbobs1=np.array(tb[).flatten() ## tbobsv=ampr_tb[3,0,:].flatten() #3 is v, 2 is h, second index is channel index ## tbobsh=ampr_tb[2,1,:].flatten() ## fig=myplot(lats,lons,tbobsv,'TBobs 10V','--',110,150,extent) ## fig=myplot(lats,lons,tbobsh,'TBobs 10H','--',110,150,extent) ## plt.ion() ## plt.show() inc_angle=ds.variables['IncidenceAngle'] #(nscan,nang) #earth incidence angles #inc_angle is -76 to +76, but usually -44 to +44 #CRTM expects positive land_fraction=ds.variables['LandFraction'] QC=ds.variables['QC'] times=np.array([dt.datetime.fromtimestamp(ts) for ts in ds.variables['Time'][:].data]) ntimes=len(times) #print('AMPR means') #for j in range(3,2,-1): #pol # for i in range(0,4): #channel # print(np.nanmean(ampr_tb[j,i,:,:])) #npol=np.size(ds.dimensions['amprChannelDim']) #nband=np.size(ds.dimensions['amprBandDim']) nchan=len(chan_set) #should be npol*nband (=8) print('file time range:') #print(dt.datetime.fromtimestamp(1*ampr_time[0,0])) print(times[0]) print(times[-1]) t1=times[0] t2=times[-1] #Define custom time range (short) #t1=dt.datetime(2024,10,23,19,50,0) #t2=dt.datetime(2024,10,23,19,52,0) #t1=dt.datetime(2024,10,23,18,30,0) #t1=dt.datetime(2024,10,23,19,30,0) #t2=dt.datetime(2024,10,23,21,30,0) #t1=dt.datetime(2024,10,30,18,38,59) #t2=dt.datetime(2024,10,31,1,00,0) #timeslice=times[:,0].data timeslice=times mask=(timeslice>t1) & (timeslice100000: print('ngood',ngood) print('nbad',nbad) pdb.set_trace() print() #Set yyyy,mm,dd,hh from ampr_time #set incidence angle #set lat/lon ptlon=ampr_lons[x,y].data ptlat=ampr_lats[x,y].data ntries=ntries+1 eia=abs(inc_angle[x,y].data) #should be degrees #This is documented as angle from surface normal. #other variables in my old retrieval: time, spatialindex,tbs is_ocean=globe.is_ocean(ptlat,ptlon) if is_ocean==False: print('skip land point') continue #CBB skip to region of interest #if ptlat>28 or ptlat<22.5: #for 10/23 ##if ptlat>27 or ptlat<24: # print('outside of box') # continue file2.write(f"x {x:4d},y {y:2d},lat {ptlat:9.2f}, lon {ptlon:9.2f}, eia {eia:7.2f}\n") #Figure out hour of day #obstime=dt.datetime.fromtimestamp(1*ampr_time[x,y]) #obstime=dt.datetime(2015,11,23,0) obstime=times[x-x0] #obsyyyymmdd=obstime.strftime('%Y%m%d') yyyy=obstime.strftime('%Y') mm=obstime.strftime('%m') dd=obstime.strftime('%d') hh=obstime.strftime('%H') #Update MERRA data if hh is new if hh != hh_old: #Read a new hour if hour has changed print('readmerra_with_cloud',hh,hh_old) [NLVL,lats,lons,tsfc,t2m,q2m,u10m,v10m,plevs,tfield,qfield,clwfield,icefield,o3field,pfield]= \ readmerra_with_cloud(yyyy,mm,dd,hh) hh_old=hh #MERRA levels are bottom up; we flip to get top down plevs=np.flip(plevs) NLVL=29 xind=np.arange(np.shape(tfield)[2]) yind=np.arange(np.shape(tfield)[1]) #Read Merra here, for first guess, or validation. #Let's start by not using Merra in the retrieval but using it for validation #We will need it for SST and Qcon #Read Merra surface type #frocean,frland,frlandice=readmerrasfctype() #missing data file #CBB Read of 3d merra fields [merra_NLVL,merra_lats,merra_lons,merra_tsfc,merra_t2m,merra_q2m,merra_u10m,merra_v10m, merra_plevs,tfield,qfield,clwfield,icefield,o3field, pfield] = readmerra_with_cloud(yyyy,mm,dd,hh) #Screen out near-land points for now ii=merra_geti(ptlon) #The ii,jj from MERRA grid (nearest integer) jj=merra_getj(ptlat) iif=getifloat(ptlon) #The ii,jj from MERRA grid (float) jjf=getjfloat(ptlat) is_ocean1=globe.is_ocean(merra_lats[jj].values, merra_lons[ii].values) is_ocean2=globe.is_ocean(merra_lats[jj].values, merra_lons[ii+1].values) is_ocean3=globe.is_ocean(merra_lats[jj+1].values, merra_lons[ii].values) is_ocean4=globe.is_ocean(merra_lats[jj+1].values, merra_lons[ii+1].values) is_ocean=(is_ocean1 and is_ocean2 and is_ocean3 and is_ocean4) if is_ocean==False: print('skip near land point') continue print('verified 4 ocean corners') #Read QCON from Merra (SST from Carol Anne Clayson?? from Merra tfield now) qcon=readmerra_qcon(yyyy,mm,dd,hh) #MERRA levels are bottom up; we flip to get top down plevs=np.flip(merra_plevs.copy()) NLVL=29 tbvec=[ampr_tb_v[0,x,y],ampr_tb_h[0,x,y], ampr_tb_v[1,x,y],ampr_tb_h[1,x,y], ampr_tb_v[2,x,y],ampr_tb_h[2,x,y], ampr_tb_v[3,x,y],ampr_tb_h[3,x,y] ] #CBB: Good place to screen out land or sea ice print('') #CBB Extract individual profile from Merra data psfc=pfield[jj,ii] #CBB This is another piece of info from Merra psfc=pseudoindex(jjf,iif,pfield) #If we are near the coast, one of our corners can be elevated (e.g. Los Angeles) #LA has a corner point at 928 mb. #In a case like that, we would like to interpolate only from ocean points. #For now, let's skip points with any land corners if psfc < 951.0: #Maybe we should do this if any corner is <951? print('Surface pressure alert:',psfc) continue levskip=0 maxlev=NLVL-1 topskip=11 nlay=29 NLVL=nlay+1 levstart=11 levend=levstart+nlay tstart=2 #Because t not flipped yet tend=30 #CBB pprof=plevs[topskip:topskip+NLVL] #now 10:975 (30 levels) #pprof=plevs[topskip:topskip+nlay+1] #now 10:975 (30 levels) #P has already been flipped pprof=plevs.copy()[levstart:levend+1] #now 10:975 (30 levels) #keep the plus one pprof[NLVL-1]=psfc #Last layer is 950 to psfc (redefine 975 layer) #CBB Get first guess from MERRA? #1. Nearest Neighbor version #Mtproflev=np.array(tfield[tstart:tend+1,jj,ii]) #(11:40) or 11 to 39, with one appended #Mtproflev=np.insert(Mtproflev,0,merra_t2m[jj,ii]) #insert t2m at bottom (currently start) now 30 levs #Mtproflev=np.flip(Mtproflev) #2. Bilinear version Mtproflev=get_vertical_profile(tfield[tstart:tend+1],xind,yind,iif,jjf) #Fix for nan at surface when psfc<1000 #Mtproflev is at 29 levels from bottom up #pprof is 30 levels from top down (already flipped) pprofup=np.flip(pprof.copy())[1:] #bottom up, no surface level for fix_tprof if np.isnan(sum(Mtproflev)): pdb.set_trace() #should not be happening # tproflev=fix_tprof(Mtproflev,pprofup) myt2m=pseudoindex(jjf,iif,t2m) Mtproflev=np.insert(Mtproflev,0,myt2m) #insert t2m at bottom (currently start) now 30 levs Mtproflev=np.flip(Mtproflev) #t2map[y,x-xstart,y]=t2m[jj,ii] t2map[y,x-xstart]=myt2m #Add interpolation for Q, O3, cloud,wind #input MERRA qfield is in kg/kg, converts to g/kg # 1. Nearest neighbor #Mqproflev=np.array(qfield[tstart:tend+1,jj,ii]) #Mqproflev=np.insert(Mqproflev,0,merra_q2m[jj,ii]) #Mqproflev=1000*np.flip(Mqproflev) #now in g/kg #logMqproflev=np.log(Mqproflev) #2 Bilinear Mqproflev=get_vertical_profile(qfield[tstart:tend+1],xind,yind,iif,jjf) Mqproflev=np.insert(Mqproflev,0,pseudoindex(jjf,iif,q2m)) #insert q2m at bottom (currently start) now 30 levs Mqproflev=1000*np.flip(Mqproflev) #converts to g/kg logMqproflev=np.log(Mqproflev) #1 Nearest neighbor for ozone #oprof=np.array(o3field[tstart:tend+1,jj,ii]) #oprof=np.insert(oprof,0,oprof[0]) #oprof=np.flip(oprof) #oprof=0.5*(oprof[0:-1]+oprof[1:]) # #myfill=np.repeat(oprof[-1],levskip) # #oprof=np.insert(oprof,nlay+1,myfill) #2 Bilinear for ozone #print('o3field original',o3field[tstart:tend+1,jj,ii]*1000000) oprof=get_vertical_profile(o3field[tstart:tend+1],xind,yind,iif,jjf) #print('new o3field',oprof*1000000) oprof=np.insert(oprof,0,oprof[0]) #repeat bottom value for near-surface oprof=np.flip(oprof) #oprof comes in as mass mixing ratio (kg/kg) and needs to be in volume mixing ratio (mol 03/mol air) #Multiply by MW of air/MW of O3. May also need factor of 1M to get in ppmv. **** oprof=oprof*28.97/24.0 #1. Nearest neighbor for cloud #Mclwproflev=np.array(clwfield[tstart:tend+1,jj,ii]) # #Mclwproflev=np.flip(Mclwproflev) #Miceproflev=np.array(icefield[tstart:tend+1,jj,ii]) # #Miceproflev=np.flip(Miceproflev) #2. Bilinear for cloud #print('clwfield at ii/jj',clwfield[tstart:tend+1,jj,ii]*1000000) Mclwproflev=get_vertical_profile(clwfield[tstart:tend+1],xind,yind,iif,jjf) #print('updated Mclwproflev',Mclwproflev*1000000) Mclwproflev=np.flip(Mclwproflev) Miceproflev=get_vertical_profile(icefield[tstart:tend+1],xind,yind,iif,jjf) Miceproflev=np.flip(Miceproflev) #print('oprof*1M is ',oprof*1.0e6) #plevs=1000,975,950,925... #psfc=1010. n=0. nlay=NVL #psfc=990. n=1 nlay=NLVL-1 #psfc=962. n=3 #The profiles passed will be #The layer from psfc to plevs(levskip) #The layer from plevs(levskip) to plevs(levskip+1) #... #ending at plevs(NLVL) (output level nlay-1) #Nearest neighbor #MSWS = np.sqrt(merra_u10m[jj,ii]**2+merra_v10m[jj,ii]**2) #m/s #Fix missing levels here (e.g. when psfc<1000) if np.isnan(sum(Mtproflev)): print('nan in profile') print('psfc:',psfc) print('Mtproflev:',Mtproflev) print('pfield:',pfield[jj:jj+2, ii:ii+2]) pdb.set_trace() print() #Bilinear u10pt=pseudoindex(jjf,iif,u10m) v10pt=pseudoindex(jjf,iif,v10m) MSWS = np.sqrt(u10pt**2+v10pt**2) #m/s mytsfc=pseudoindex(jjf,iif,tsfc) #CBB updated to here #First set background profile from bkopt. Then set fg from fgopt, overwriting profile. #We have tprof,logqprof,clwprof,iceprof from reading merra #First process the merra prof because we need it whether or not it is our bk/first guess Mtprof=0.5*(Mtproflev[0:-1]+Mtproflev[1:]) #pdb.set_trace() if(np.isnan(np.sum(Mtprof))): #check for any nan print('error in Mtprof',Mtprof) status=-5 #skipped logMqprof=0.5*(logMqproflev[0:-1]+logMqproflev[1:]) #if cloudopt=='Yes': #CLW/Ice do not have a surface value so extend the last level before averaging Mclwproflev=np.append(Mclwproflev,Mclwproflev[-1]) Mclwprof=0.5*(Mclwproflev[0:-1]+Mclwproflev[1:]) Miceproflev=np.append(Miceproflev,Miceproflev[-1]) Miceprof=0.5*(Miceproflev[0:-1]+Miceproflev[1:]) #else: # clwprof=np.zeros(29) # iceprof=np.zeros(29) #CBB adding this here to compare with fd_model_whymsie #CBB fdmodel (wrapper subroutine) and runDirect (on merraprof) give # the same results. Not sure which is better/simpler/more efficient. # fdmodel is a straightforward wrapper for runDirect tbmerra=fdmodel(sensor_id,chan_set,ptlat,ptlon, \ pprof,Mtprof,logMqprof,Mclwprof,Miceprof,mytsfc,MSWS,eia) tbmerragrid[y,x-xstart,:]=tbmerra #channel order is 10.7,19.35,37.1,85.5 (each V, then H) print('tb from merra profile in main:',tbmerra) print('initial SetProfile with Merra, ncloud=',ncloud) #SetProfile called now to get merraprof, merraclwsum (only needed for validation?) merraprof=SetProfile( NLVL, pprof, Mtprof, np.exp(logMqprof), Mclwprof, Miceprof, oprof, mytsfc, MSWS, yyyy,mm,dd,eia,ncloud=ncloud) #if y==ny-10: # pdb.set_trace() print('cloud in merraprof') print(merraprof.clouds[0,:,0,0]*1000) clwbk=Mclwprof #CBB corrected error here. 6/10/22 icebk=Miceprof if ncloud>0: merraclwsum=np.sum(merraprof.clouds[0,:,0,0]) print('100*liquid cloud:',100*merraprof.clouds[0,:,0,0]) clwsum=merraclwsum else: clwsum=0 merraclwsum=0. if ncloud>1: merraicesum=np.sum(merraprof.clouds[0,:,1,0]) print('100* ice cloud:',100*merraprof.clouds[0,:,1,0]) icesum=merraicesum else: icesum=0 merraicesum=0. print('CLWSUM,ICESUM',clwsum,icesum) #profiles=merraprof #The initial bins are determined by Merra data mysst=float(merra_tsfc[jj,ii]) qcon0=float(qcon[jj,ii]) myqcon=pseudoindex(jjf,iif,qcon) print('qcon:',qcon0,myqcon) isst=int((mysst-271.15)/4) jqcon=int(((myqcon+2)/0.4)) #Set initial values based on bin mean values tprof=bin_avgprof[jqcon,isst,0:29] logqprof=bin_avgprof[jqcon,isst,29:58] bkcov=bin_bkcov[jqcon,isst,:,:] evecs=bin_evecs[jqcon,isst,:,:] evals=bin_evals[jqcon,isst,:] mysws=bin_sws[jqcon,isst] #This is the first guess SWS ## qprof=np.exp(logqprof) ## #set oprof to mean oprof? ## #SWS: Use Merra? Or just mean? Do we have SWS by bin? (Currently Merra) ## profile=SetProfile( NLVL, pprof, tprof, qprof, clwprof, iceprof, oprof, mytsfc, SWS, yyyy,mm,dd) xarr=np.zeros(6+ncloud) #State vector: coefficients of each EOF (5), plus SWS and maybe clw/ice profiles (2) xarr[5]=np.log(mysws) #xarr[6]=0.0 #log of clw coefficient (zero by default) #xarr[7]=0.0 #log of ice coefficient (zero by default) clwbk=avgclwprof[jqcon,isst,:] #Mean for this bin icebk=avgiceprof[jqcon,isst,:] #CBB This is old logic. If cloudy, use the Merra profile for cloud background, otherwise #use the global (or bin) mean. Now we are just starting with the bin mean. ## if clwsum<.001: #what is a good threshold? one val is .0057 for clwsum and .044 icesum ## clwbk=avgclwprof[jqcon,isst,:] ## xarr[6]=-5.0 #ln(coeff) for cloud water ## print('near zero clw in merra prof, using mean clw profile') ## if icesum<.001: ## icebk=avgiceprof[jqcon,isst,:] ## xarr[7]=-5.0 #ln(coeff) for cloud ice ## print('near zero ice in merra prof, using mean ice profile') ## CBB Took out loop for cloud options. ## for itry in range(0,1): #CBB Loop can include both cloudstatus options if desired. cloudstatus=[False,False] print('calling ret1dvar with cloudstatus ',cloudstatus) print('and tbvec',tbvec) #CBB Does "profile[s] and merraprof need to be set here?? Seems redundant but verify. #taking out passing of profile[s] #CBB 5/26/26 taking out cloudstatus for now (was last argument, a 2-length boolean) #CBB 9/11/26 Do the Merra forward model before ret1dvar # #### print('initialize CRTM') #### crtmOb = pyCRTM() #Seg fault if omitted #### #CBB Aerosol profile removed to avoid "phase errors" #### # Depends on presence or absence of aerosols and clouds #### #crtmOb.AerosolCoeff_File='AerosolCoeff.bin' #CBB TEST #### crtmOb.sensor_id = sensor_id #### crtmOb.profiles = merraprof #### crtmOb.sensor_id = sensor_id #### crtmOb.nThreads = 1 #### print('crtmOb.loadInst()................................') #### crtmOb.loadInst() #load the instrument, set number of channels #### print('Total channels loaded:',crtmOb.nChanTotal) #### crtmOb.channelSubset = chan_set+1 #### print('call runDirect................................') #### crtmOb.runDirect() #CBB This is what runs the Forward Model #### #merratb = crtmOb.Bt[0,chan_set] #### tbmerra2 = crtmOb.Bt[0,:] #was merratb (WRONG name), this is from first guess scanpos=y #for printout in ret1dvar status,prof,it,tbpen,xpen,mtbpen,mxpen,chisq,tiwv,clwsum,icesum,SWS,tbinit,merra_tiwv = \ ret1dvar( sensor_id,chan_set,ptlat,ptlon, \ tbvec,xarr,bkcov,evecs,evals,pprof,tprof,logqprof,clwbk,icebk,mysst,mysws,ff, \ merraprof,Mtprof,logMqprof,eia,ncloud,tbmerra) tbinitgrid[y,x-xstart,:]=tbinit tbobsgrid[y,x-xstart,:]=tbvec tbpengrid[y,x-xstart]=tbpen xpengrid[y,x-xstart]=xpen mtbpengrid[y,x-xstart]=mtbpen statusgrid[y,x-xstart]=status #nprof=1 flightdata['lat'][y,x-xstart]=ptlat flightdata['lon'][y,x-xstart]=ptlon flightdata['eia'][y,x-xstart]=eia flightdata['scanpos'][y,x-xstart]=scanpos flightdata['tbinit'][y,x-xstart,:]=tbinit flightdata['tbobs'][y,x-xstart,:]=tbvec flightdata['tbmerra'][y,x-xstart,:]=tbmerra flightdata['clw'][y,x-xstart]=clwsum flightdata['sws'][y,x-xstart]=SWS flightdata['iwv'][y,x-xstart]=tiwv flightdata['m_clw'][y,x-xstart]=merraclwsum flightdata['m_sws'][y,x-xstart]=MSWS flightdata['m_iwv'][y,x-xstart]=merra_tiwv flightdata['xpen'][y,x-xstart]=xpen flightdata['tbpen'][y,x-xstart]=tbpen flightdata['mtbpen'][y,x-xstart]=mtbpen flightdata['status'][y,x-xstart]=status if(status>=0): print('GOOD RETRIEVAL after ',it+1,' iterations',ngood,nbad,nskip) ngood=ngood+1 #arg=input("Press enter, or type 9 to stop.") #if arg=='9': # pdb.set_trace() file1.write(f"{ptlat:9.2f} {ptlon:9.2f} {status:3d} {tiwv:9.3f} {clwsum:9.3f} {mysws:5.2f}\n") file2.write(f"{ptlat:9.2f} {ptlon:9.2f} {status:3d} {tiwv:9.3f} {clwsum:9.3f} {mysws:5.2f}\n") file2.write('') print('.............................................') continue #exit loop elif (status==-4): print('SKIPPING') nskip=nskip+1 continue else: print('BAD RETRIEVAL',status,ngood,nbad,nskip) nbad=nbad+1 #arg=input("Press enter, or type 9 to stop.") #if arg=='9': # pdb.set_trace() TIWV=-1 clwsum=-1 mysws=-1 file1.write(f"{ptlat:9.2f} {ptlon:9.2f} status:{status:3d}\n" ) file2.write(f"{ptlat:9.2f} {ptlon:9.2f} status:{status:3d}\n" ) file2.write('') print('Try again....................................') continue file1.close() with open(file4,'wb') as file: pickle.dump(flightdata,file) print('ngood',ngood) print('nbad',nbad) #Plot these #extent=[t1,t2,0,ny-1] tbdiffs=tbinitgrid[:,:,0]-tbobsgrid[:,:,0] scanplot(tbdiffs,'Ch 1 initial difference', cmap='seismic',vmin=-20,vmax=20) scanplot(tbobsgrid[:,:,0],'Ch 1 TBobs', cmap='viridis',vmin=110,vmax=150) scanplot(tbinitgrid[:,:,0],'Ch 1 TB (first guess)', cmap='viridis',vmin=110,vmax=150) tbdiffs=tbinitgrid[:,:,2]-tbobsgrid[:,:,2] scanplot(tbdiffs,'Ch 3 initial difference', cmap='seismic',vmin=-20,vmax=20) #tbdiffs=tbinitgrid[:,:,4]-tbobsgrid[:,:,4] #scanplot(tbdiffs,'Ch 5 initial difference', cmap='seismic',vmin=-20,vmax=20) #tbdiffs=tbinitgrid[:,:,6]-tbobsgrid[:,:,6] #scanplot(tbdiffs,'Ch 7 initial difference', cmap='seismic',vmin=-20,vmax=20) scanplot(tbpengrid,'TBpen', cmap='viridis',vmin=0,vmax=5) scanplot(xpengrid,'Xpen', cmap='viridis',vmin=0,vmax=7) scanplot(mtbpengrid,'Merra TBpen', cmap='viridis',vmin=0,vmax=20) scanplot(statusgrid,'Status', cmap='viridis',vmin=-5,vmax=5) #Test extracting gridded data from 2d array of dicts scanplot(flightdata['sws'],'SWS', cmap='viridis',vmin=0,vmax=10) scanplot(flightdata['m_sws'],'Merra SWS', cmap='viridis',vmin=0,vmax=10) extent=[-119,-117,25,27] fig=myplot(lats,lons,flightdata['sws'],'SWS','m/s',0,10,extent) fig=myplot(lats,lons,flightdata['m_sws'],'SWS','m/s',0,10,extent) print('end of program') pdb.set_trace() print() #[NLVL,Mlats,Mlons,Mtsfc,Mt2m,Mq2m,Mu10m,Mv10m,Mplevs,Mtfield,Mqfield,Mclwfield,Micefield,Mo3field,Mpfield]= \ #Kmat test code # #qprof=np.exp(logqprof) # #profiles.T[n,:] = tprof # [0:29] # #profiles.Q[n,:] = qprof # # #kEmissivity = crtmOb.surfEmisRefl[0,:] #CBB OK? Not really needed unless in state vector. # #Calculate the next profile vector # #CBB This is how to do it using the full profile vector # # #CBB: kmat has to be calculated from kTb######################vvvvvv # #kmat=kTb #CBB: check orientation, also this needs a factor added # #ktb is (22,58) # #kmat needs to be (22,5) # #kmat=(dTb/dx)=(dtB/dprof)(dprof/dx)=(kTb)(dprof/dx) # #How to calc last term? # #prof=sum(xarr*evecs) # #dprof/dx=evecs # #CBB Old code, probably not right, I think this uses the "nonlinear" version #prod1=np.matmul(kmat,(xarr-oldxarr)) #prod2=np.matmul((bkcovinv+kt_rinv_k),kt_rinv) #prod2.shape #prod3=np.matmul(prod2,(yobs-Tbcalc[0,:]-prod1)) #prod3.shape #oldxarr=xarr #xarr=origxarr+prod3 #Notes #crtmOb.frequencyGHz has channel frequencies (what about polarizations?) #P,Pi,T,Q #cloud info: cloudType[ntype], cloudFraction[nlev] #Jacobians TK, QK, O3K, SkinK, SurfaceEmisK, WindSpeed