Skip to content
Snippets Groups Projects
Commit e5716b9a authored by brunner's avatar brunner
Browse files

lambdas are not 0

parent a081bb37
No related branches found
No related tags found
No related merge requests found
...@@ -347,11 +347,23 @@ class Optimizer(object): ...@@ -347,11 +347,23 @@ class Optimizer(object):
alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n])) alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n]))
# if all(self.x[:] + self.KG[:] * res >= 0.):
self.x[:] = self.x + self.KG[:] * res self.x[:] = self.x + self.KG[:] * res
self.x[self.x<0.] = 0. # cut off negative values, COSMO don't like negative fluxes # pavle self.x[self.x<0.] = 0. # cut off negative values, COSMO don't like negative fluxes # pavle
for r in range(self.nmembers): for r in range(self.nmembers):
self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r]) X_prime_temp = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
X_prime_temp[X_prime_temp+self.x < 0.] = 0.
# if all(self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r]) + self.x[:] >= 0.):
self.X_prime[:, r] = X_prime_temp
#self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
# continue
# else:
# self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])
# for p in range(self.nparams):
# if self.X_prime[p, r]+self.x[p]<0.:
# self.X_prime[p, r]=0.
#WP !!!! Very important to first do all obervations from n=1 through the end, and only then update 1,...,n. The current observation #WP !!!! Very important to first do all obervations from n=1 through the end, and only then update 1,...,n. The current observation
#WP should always be updated last because it features in the loop of the adjustments !!!! #WP should always be updated last because it features in the loop of the adjustments !!!!
......
...@@ -266,7 +266,7 @@ def prepare_state(dacycle, statevector): ...@@ -266,7 +266,7 @@ def prepare_state(dacycle, statevector):
current_sv = os.path.join(dacycle['dir.restart'], 'savestate_%s.nc' % dacycle['time.start'].strftime('%Y%m%d')) current_sv = os.path.join(dacycle['dir.restart'], 'savestate_%s.nc' % dacycle['time.start'].strftime('%Y%m%d'))
statevector.write_to_file(current_sv, 'prior') # write prior info statevector.write_to_file(current_sv, 'prior') # write prior info
nlag = int(dacycle['time.nlag']) # nlag = int(dacycle['time.nlag'])
# for l in range(0,nlag): # pavle added from here # for l in range(0,nlag): # pavle added from here
## statevector.write_members_to_file(l, dacycle['restartmap.dir']) ## statevector.write_members_to_file(l, dacycle['restartmap.dir'])
# statevector.write_members_for_cosmo(l, dacycle['restartmap.dir']) # statevector.write_members_for_cosmo(l, dacycle['restartmap.dir'])
......
...@@ -290,8 +290,7 @@ class StateVector(object): ...@@ -290,8 +290,7 @@ class StateVector(object):
newmember.param_values = newmean.flatten() # no deviations newmember.param_values = newmean.flatten() # no deviations
self.ensemble_members[lag].append(newmember) self.ensemble_members[lag].append(newmember)
# ifile = Dataset('/scratch/snx3000/parsenov/40_9reg/bckp/output/20130401/lambda.nc', mode='r') ifile = Dataset('/scratch/snx3000/parsenov/octe/optimizer.20130401.nc', mode='r')
ifile = Dataset('/scratch/snx3000/parsenov/9reg/optimizer.20130401.nc', mode='r')
# par = ifile.variables['lambda'][:] # par = ifile.variables['lambda'][:]
par = ifile.variables['statevectordeviations_prior'][:] par = ifile.variables['statevectordeviations_prior'][:]
par = par.reshape(181,2,21) par = par.reshape(181,2,21)
...@@ -313,8 +312,8 @@ class StateVector(object): ...@@ -313,8 +312,8 @@ class StateVector(object):
# newmember.param_values[newmember.param_values>2.] = 2. # newmember.param_values[newmember.param_values>2.] = 2.
newmember.param_values = par[:, lag, member] newmember.param_values = par[:, lag, member]
# newmember.param_values = par[lag, member, :] # newmember.param_values = par[lag, member, :]
if member==1: # if member==1:
print(par[:, lag, member]) # print(par[:, lag, member])
self.ensemble_members[lag].append(newmember) self.ensemble_members[lag].append(newmember)
logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1))) logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))
......
...@@ -38,7 +38,7 @@ abs.time.start : 2013-04-01 00:00:00 ...@@ -38,7 +38,7 @@ abs.time.start : 2013-04-01 00:00:00
! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE ! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE
time.restart : F time.restart : F
da.restart.tstamp : 2013-01-01 00:00:00 da.restart.tstamp : 2013-04-01 00:00:00
! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1 ! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment