Skip to content

GitLab

Commit 129a712d authored by Audrey Fortems-Cheiney's avatar Audrey Fortems-Cheiney
Browse files

sim taken into account for TL

parent fde88e71
Hide whitespace changes
Inline Side-by-side
model_sources/chimereGES/src/initio/inidoms.f90
View file @ 129a712d
......@@ -42,7 +42,7 @@ subroutine inidoms
nobs_glo=0
end if
if(ad.eq.0)allocate ( tabobs_glo(nobs_glo,11) )
if(ad.eq.0)allocate ( tabobs_glo(nobs_glo,12) )
if(ad.eq.1)allocate ( tabobs_glo(nobs_glo,12) )
tabobs_glo(:, :) = 0.
......@@ -76,7 +76,7 @@ subroutine inidoms
enddo
do ip=1,ndoms
print*,'tile number',ip,' number of obs=',nobs(ip)
if(ad.eq.0)allocate ( tabobs(nobs(ip),12) )
if(ad.eq.0)allocate ( tabobs(nobs(ip),13) )
if(ad.eq.1) allocate ( tabobs(nobs(ip),13) )
nb=0
do i=1,nobs_glo
......@@ -90,7 +90,7 @@ subroutine inidoms
if(ad.eq.0) then
tabobs(nb,1:7)=tabobs_glo(i,:)
tabobs(nb,8)=dble(i) ! for reconstitution in the end
tabobs(nb,9:12)=tabobs_glo(i,8:11)
tabobs(nb,9:13)=tabobs_glo(i,8:12)
endif
if(ad.eq.1) then
tabobs(nb,1:8)=tabobs_glo(i,:)
......@@ -102,7 +102,7 @@ subroutine inidoms
enddo
! send to workers
call mpi_send(nobs(ip),1,mpi_integer,ip, ias_nobs, mpi_comm_world,ierr)
if(ad.eq.0)call mpi_send(tabobs,nobs(ip)*12,mpi_double_precision,ip, &
if(ad.eq.0)call mpi_send(tabobs,nobs(ip)*13,mpi_double_precision,ip, &
ias_tabobs,mpi_comm_world,ierr)
if(ad.eq.1)call mpi_send(tabobs,nobs(ip)*13,mpi_double_precision,ip, &
ias_tabobs,mpi_comm_world,ierr)
......
model_sources/chimereGES/src_tl/main/integrun_tl.f90
View file @ 129a712d
......@@ -80,19 +80,19 @@ subroutine integrun_tl
!IP receiving and merging obs
do ip=1,ndoms
allocate ( tabobs(nobs(ip),12) )
call mpi_recv(tabobs,nobs(ip)*12,mpi_double_precision,ip, &
allocate ( tabobs(nobs(ip),13) )
call mpi_recv(tabobs,nobs(ip)*13,mpi_double_precision,ip, &
iar_tabobs,mpi_comm_world,mpi_status_ignore,ierr)
do i=1,nobs(ip)
nb=tabobs(i,8)
tabobs_glo(nb,1:7) = tabobs(i,1:7)
tabobs_glo(nb,8:11) = tabobs(i,9:12)
tabobs_glo(nb,8:12) = tabobs(i,9:13)
enddo
deallocate(tabobs)
enddo
! writing
do i=1,nobs_glo
write(2,1001)int(tabobs_glo(i,1:5)),species(int(tabobs_glo(i,6)))%name,tabobs_glo(i,7:11)
write(2,1001)int(tabobs_glo(i,1:5)),species(int(tabobs_glo(i,6)))%name,tabobs_glo(i,7:12)
enddo
close(2)
......
model_sources/chimereGES/src_tl/model/iniworker_tl.f90
View file @ 129a712d
......@@ -17,8 +17,8 @@ subroutine iniworker_tl
call mpi_recv(nobs,1,mpi_integer,0,ias_nobs, &
mpi_comm_world,mpi_status_ignore,ierr)
print*,'worker numero',rank,' nb obs=',nobs
allocate(tabobs(nobs,12))
call mpi_recv(tabobs,nobs*12,mpi_double_precision,0,ias_tabobs, &
allocate(tabobs(nobs,13))
call mpi_recv(tabobs,nobs*13,mpi_double_precision,0,ias_tabobs, &
mpi_comm_world,mpi_status_ignore,ierr)
! fin IP
call mpi_barrier(mpi_comm_world,ierr)
......
model_sources/chimereGES/src_tl/model/worker_tl.f90
View file @ 129a712d
......@@ -91,10 +91,12 @@ subroutine worker_tl
if (species(ns)%fspec > dzero) then
! conversion molec.cm-3 to the unit specified by the user
tabobs(i,7)= species(ns)%fspec*conc_tl(ns,izo,ime,ivert)
tabobs(i, 13) = species(ns)%fspec*conc(ns,izo,ime,ivert)
else
! conversion molec.cm-3 to ppb
fppb = 1d9/airmloc(izo,ime,ivert)
tabobs(i,7) = conc_tl(ns,izo,ime,ivert) * fppb
tabobs(i, 13) = conc(ns,izo, ime, ivert) * fppb
endif
if (ivert.eq.1) then
......@@ -114,8 +116,6 @@ subroutine worker_tl
tabobs(i, 10) = ddp !in Pa
tabobs(i, 11) = airmloc(izo, ime, ivert) !in molec.cm-3
tabobs(i, 12) = thlayloc(izo, ime, ivert) !in cm
endif ! ihourrun,np,species
enddo ! nobs
......@@ -140,7 +140,7 @@ subroutine worker_tl
enddo ! nh=1,nhourrun
!IP sending tabobs to the master
if (ad.eq.0) call mpi_send(tabobs, nobs * 12, mpi_double_precision, 0, &
if (ad.eq.0) call mpi_send(tabobs, nobs * 13, mpi_double_precision, 0, &
iar_tabobs, mpi_comm_world, ierr)
if(ad.eq.1)call mpi_send(tabobs, nobs * 13, mpi_double_precision, 0, &
iar_tabobs, mpi_comm_world, ierr)
......
pycif/plugins/models/chimere/outputs2native.py
View file @ 129a712d
......@@ -41,6 +41,7 @@ def outputs2native(self, data2dump, input_type,
dp = pd.read_csv(sim_file, delim_whitespace=True, header=None)[8]
airm = pd.read_csv(sim_file, delim_whitespace=True, header=None)[9]
hlay = pd.read_csv(sim_file, delim_whitespace=True, header=None)[10]
simfwd = pd.read_csv(sim_file, delim_whitespace=True, header=None)[11]
# Loop over observations in active species
mask = self.dataobs['parameter'].str.upper() \
.isin(self.chemistry.acspecies.attributes)
......@@ -65,6 +66,9 @@ def outputs2native(self, data2dump, input_type,
self.dataobs.loc[mask, 'hlay'] = \
[hlay.iloc[k:k + dt].sum()
for k, dt in zip(inds, self.dataobs.loc[mask, 'dtstep'])]
self.dataobs.loc[mask, 'simfwd'] = \
[simfwd.iloc[k:k + dt].sum()
for k, dt in zip(inds, self.dataobs.loc[mask, 'dtstep'])]
return self.dataobs
......
pycif/plugins/tranforms/satellites/native2obsvect.py
View file @ 129a712d
......@@ -33,25 +33,32 @@ def native2obsvect(transf, xmod, mapper, mod_input,
iq1 = y0['iq1']
nblinfo = y0['nblinfo']
list_satIDs = iq1.loc[y0['is_satellite']].unique()
ds_p = datastore.set_index('indorig')[['pressure', 'dp', 'airm', 'hlay']]
if 'sim' in datastore:
ds_p = datastore.set_index('indorig')[['pressure', 'dp', 'airm', 'hlay']]
if 'sim_tl' in datastore:
ds_p = datastore.set_index('indorig')[['pressure', 'dp', 'airm', 'hlay','simfwd']]
print('AAAAAA', datastore)
for satID in list_satIDs:
satmask = iq1 == satID
nobs = np.sum(satmask)
nblloc = nblinfo.loc[satmask].values - 1
print('satID',satID)
# Stacking output datastore into levels * nobs
native_ind_stack = np.flatnonzero(ref_indexes)[satmask] \
+ np.arange(transf.model.domain.nlev)[:, np.newaxis]
# If all nans in datasim, meaning that the species was not simulated
sim = datastore['sim'].values[native_ind_stack]
if 'sim' in datastore:
sim = datastore['sim'].values[native_ind_stack]
if 'sim_tl' in datastore:
sim = datastore['sim_tl'].values[native_ind_stack]
simfwd = datastore['simfwd'].values[native_ind_stack]
if not np.any(~np.isnan(sim)):
continue
# Grouping all data from this satellite
datasim = xr.Dataset(
if 'sim' in datastore:
datasim = xr.Dataset(
{'pressure': (['level', 'index'],
np.log(ds_p['pressure'].values[native_ind_stack])),
'dp': (['level', 'index'], ds_p['dp'].values[native_ind_stack]),
......@@ -62,12 +69,22 @@ def native2obsvect(transf, xmod, mapper, mod_input,
'sim': (['level', 'index'], sim)},
coords={'index': nblloc,
'level': np.arange(transf.model.domain.nlev)}
)
)
if 'sim_tl' in datastore:
datasim['sim_tl'] = (['level', 'index'],
datastore['sim_tl'].values[native_ind_stack])
datasim = xr.Dataset(
{'pressure': (['level', 'index'],
np.log(ds_p['pressure'].values[native_ind_stack])),
'dp': (['level', 'index'], ds_p['dp'].values[native_ind_stack]),
'airm': (
['level', 'index'], ds_p['airm'].values[native_ind_stack]),
'hlay': (
['level', 'index'], ds_p['hlay'].values[native_ind_stack]),
'sim': (['level', 'index'], sim),
'simfwd': (['level', 'index'], simfwd)},
coords={'index': nblloc,
'level': np.arange(transf.model.domain.nlev)}
)
# convert CHIMERE fields to the correct unit
# from ppb to molec.cm-2 if the satellite product is a column
if getattr(transf.available_satellites, satID).product == 'column':
......@@ -78,12 +95,13 @@ def native2obsvect(transf, xmod, mapper, mod_input,
file_aks = ddi.strftime(
'{}/obsvect/satellites/infos_{}%Y%m%d%H%M.nc'
.format(transf.workdir, satID))
print(file_aks)
print('infos file',file_aks)
try:
#try:
if True:
sat_aks = read_datastore(file_aks).to_xarray()
except IOError:
#except IOError:
# Assumes total columns?
# datastore['qdp'] = datastore['sim'] * datastore['dp']
# groups = datastore.groupby(['indorig'])
......@@ -95,10 +113,10 @@ def native2obsvect(transf, xmod, mapper, mod_input,
# groups = datastore.groupby(['indorig'])
# y0.loc[:, 'sim_tl'] += \
# groups['qdp'].sum().values / groups['dp'].sum().values
continue
# continue
aks = sat_aks['ak'][nblloc, ::-1].T
print('aks',aks)
if getattr(transf.available_satellites, satID).pressure == 'Pa':
pavgs = sat_aks['pavg'][nblloc, ::-1].T
else:
......@@ -118,6 +136,7 @@ def native2obsvect(transf, xmod, mapper, mod_input,
'level': np.arange(aks.level.size)})
qa0lus = sat_aks['qa0lu'][nblloc, ::-1].T
print('qa0lus',qa0lus)
# Interpolating simulated values to averaging kernel pressures
# Doing it by chunk to fasten the process
......@@ -126,7 +145,6 @@ def native2obsvect(transf, xmod, mapper, mod_input,
# of scipy automatic parallelisation
# 50 chunks seems to be fairly efficient
sim_ak = 0. * pavgs_mid
sim_ak_tl = 0. * pavgs_mid
chunks = np.linspace(0, nobs, num=50, dtype=int)
cropstrato = getattr(getattr(transf.available_satellites, satID),
'cropstrato', 0)
......@@ -141,15 +159,18 @@ def native2obsvect(transf, xmod, mapper, mod_input,
pavgs_mid[:, k1:k2].values, cropstrato)
# print('pppppppp', xlow,xhigh, alphalow,alphahigh)
# Applying coefficients
meshout = np.array(pavgs_mid.shape[0] * [list(range(k2 - k1))])
if 'sim_tl' in datastore:
simfwd = datasim['simfwd'][:, k1:k2].values
sim = datasim['sim'][:, k1:k2].values
sim_ak[:, k1:k2] = alphalow * sim[xlow, meshout] \
+ alphahigh * sim[xhigh, meshout]
if 'sim_tl' in datasim:
sim_tl = datasim['sim_tl'][:, k1:k2].values
sim_ak_tl[:, k1:k2] = alphalow * sim_tl[xlow, meshout] \
+ alphahigh * sim_tl[xhigh, meshout]
if 'sim_tl' in datastore:
sim_fwd_ak[:, k1:k2] = alphalow * simfwd[xlow, meshout] \
+ alphahigh * simfwd[xhigh, meshout]
# # Correction with the pressure thickness
# # WARNING: there is an inconsistency in the number of levels
......@@ -170,14 +191,15 @@ def native2obsvect(transf, xmod, mapper, mod_input,
print('nbformula:', nbformula)
print('chosenlev:', chosenlevel)
y0.loc[satmask, 'sim'] = apply_ak(sim_ak, dpavgs, aks.values, nbformula,
if 'sim' in datastore:
y0.loc[satmask, 'sim'] = apply_ak(sim_ak, dpavgs, aks.values, nbformula,
qa0lus.values, chosenlevel)
print('ffffffffffffffff', y0.loc[satmask, 'sim'])
if 'sim_tl' in datasim:
y0.loc[satmask, 'sim_tl'] = apply_ak_tl(sim_ak_tl, dpavgs,
if 'sim_tl' in datastore:
y0.loc[satmask, 'sim_tl'] = apply_ak_tl(sim_ak, dpavgs,
aks.values, nbformula,
qa0lus.values, chosenlevel,
sim_ak)
sim_fwd_ak)
# Save forward datastore for later use by adjoint
file_monit = ddi.strftime('{}/chain/monit_%Y%m%d%H%M.nc'
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment