[Dart-dev] [3477] DART/trunk/models/MITgcm_ocean:
Given a directory with N snapshot files and a DART input.nml file,
nancy at ucar.edu
nancy at ucar.edu
Mon Jul 28 16:58:45 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080728/3c30ad3b/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/MITgcm_ocean/inputs/data
===================================================================
--- DART/trunk/models/MITgcm_ocean/inputs/data 2008-07-28 17:50:02 UTC (rev 3476)
+++ DART/trunk/models/MITgcm_ocean/inputs/data 2008-07-28 22:58:45 UTC (rev 3477)
@@ -51,25 +51,27 @@
&
# Time stepping parameters
-# dumpFreq = 86400.,
-# taveFreq = 86400.,
-# to get 20 ensemble members ... from a 14 day run ...
-# I'm going to dump more often than every day.
-# endTime = 1209600., 14 days
-# endTime = 172800., two days
-#dumpFreq = 43200.,
-#taveFreq = 43200.,
+# endTime = 172800., 2 days
+# endTime = 1209600., 14 days
+# endTime = 1728000., 20 days
+# to get N ensemble members ... from a 20 day run ...
+# simply change the output frequency to generate the appropriate number
+# of snapshot files and then use MakeInitialEnsemble.csh
+# dumpFreq = 86400., (daily assimilation cycle)
+# taveFreq = 86400., (daily assimilation cycle)
+# dumpFreq = 43200., (12-hourly assimilation cycle)
+# taveFreq = 43200., (12-hourly assimilation cycle)
&PARM03
startTime = 0.,
- endTime = 172800.,
+ endTime = 86400.,
deltaTmom = 900.,
deltaTtracer = 900.,
deltaTClock = 900.,
abEps = 0.1,
chkptFreq = 0.,
pChkptFreq = 1209600.,
- dumpFreq = 172800.,
- taveFreq = 172800.,
+ dumpFreq = 86400.,
+ taveFreq = 86400.,
tauThetaClimRelax = 0.,
tauSaltClimRelax = 2592000.,
monitorFreq=259200000.0,
Added: DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m
===================================================================
--- DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m (rev 0)
+++ DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m 2008-07-28 22:58:45 UTC (rev 3477)
@@ -0,0 +1,79 @@
+function CheckICS()
+% CheckICS
+%
+%
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+%-------------------------------------------------------------------------------
+dsize = [256 225 40];
+%-------------------------------------------------------------------------------
+
+mitbase = '/ptmp/thoar/MITgcm/ics';
+fname = sprintf('%s/filter_ics',mitbase);
+ic00 = rdinit(fname,dsize);
+orgmask = find( ic00.data == 0.0 ); % should be dry locations
+
+if ( 1 == 1 )
+
+ for i = 1:40
+ fname = sprintf('%s/ens_mem_%03d',mitbase,i);
+ copy = rdinit(fname,dsize);
+ dry = copy.data(orgmask); % should all be dry ... 0.0
+ tdiff = sum(dry(:));
+ disp(sprintf('copy %02d total is %bx aka %e',i,tdiff,tdiff))
+
+ mymask = find(copy.data == 0.0);
+ disp(sprintf('mask %i length %06d mask 1 length %06d',i,length(orgmask),length(mymask)))
+ if (length(orgmask) ~= length(mymask))
+ error(sprintf('mask %i length unequal to mask 1 length',i))
+ end
+
+% bob = copy.data(:) - ic00.data(:);
+% disp(sprintf('tot diff between ens_mem_%03d and ens_mem_001 %f',i,sum(bob(isfinite(bob)))))
+
+ disp(sprintf('timestamp of ens_mem_%03d %f %f\n',i,copy.time(1),copy.time(2)))
+
+ end
+
+else
+
+
+end
+
+%======================================================================
+
+function y = rdinit(fname,dsize)
+
+if (exist(fname,'file'))
+ % disp(sprintf('Opening %s',fname))
+else
+ error('Opening %s',fname)
+end
+
+% A DART Initial Conditions file
+fid = fopen(fname,'rb','ieee-be');
+h1 = fread(fid,1,'int32');
+t = fread(fid,2,'int32');
+h2 = fread(fid,1,'int32');
+
+h3 = fread(fid, 1,'int32');
+[x count] = fread(fid,prod(dsize),'float64');
+h4 = fread(fid, 1,'int32');
+fclose(fid);
+
+if (count ~= prod(dsize))
+ error('only read %d of %d items from %s',count,prod(dsize),fname)
+end
+
+y.data = reshape(x, dsize);
+y.time = t;
Property changes on: DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Added: DART/trunk/models/MITgcm_ocean/matlab/Check_trans_pv_sv.m
===================================================================
--- DART/trunk/models/MITgcm_ocean/matlab/Check_trans_pv_sv.m (rev 0)
+++ DART/trunk/models/MITgcm_ocean/matlab/Check_trans_pv_sv.m 2008-07-28 22:58:45 UTC (rev 3477)
@@ -0,0 +1,164 @@
+function Check_trans_pv_sv
+%
+% data.cal holds the starting time information
+%
+
+dt = 900;
+startDate_1=19960101;
+startDate_2=60000;
+tbase = datenum(1996,1,1,6,0,0);
+
+timestep = 40992; % data70levels
+timestep = 1344; % data40levels
+timestep = 96; % real data 96*900 = 86400 = 1day
+
+mitdir = '/ptmp/thoar/MITgcm/adv1day_assim/advance_temp0';
+
+S = rdmds(sprintf( '%s/S.%10.10d',mitdir,timestep));
+T = rdmds(sprintf( '%s/T.%10.10d',mitdir,timestep));
+U = rdmds(sprintf( '%s/U.%10.10d',mitdir,timestep));
+V = rdmds(sprintf( '%s/V.%10.10d',mitdir,timestep));
+SSH = rdmds(sprintf('%s/Eta.%10.10d',mitdir,timestep));
+
+modelsize = prod(size(S)) + prod(size(T)) + ...
+ prod(size(U)) + prod(size(V)) + prod(size(SSH));
+
+fsize = 4+4+4+4 + 4+modelsize*8+4;
+disp(sprintf('with a modelsize of %d the file size should be %d bytes', ...
+ modelsize,fsize))
+
+toffset = (timestep*dt)/86400;
+
+[nx ny nz] = size(S);
+
+fname = 'assim_model_state_ud';
+fname = '/ptmp/thoar/MITgcm/adv1day_assim/assim_model_state_ud.0020';
+fid = fopen(fname,'rb','ieee-be');
+
+trec1 = fread(fid,1,'int32');
+seconds = fread(fid,1,'int32');
+days = fread(fid,1,'int32');
+trecN = fread(fid,1,'int32');
+
+if (trec1 ~= trecN)
+ error('first record mismatch')
+end
+
+rec1 = fread(fid, 1,'int32');
+datmat = fread(fid,modelsize,'float64');
+recN = fread(fid, 1,'int32');
+fclose(fid);
+
+if (rec1 ~= recN)
+ error('second record mismatch')
+end
+
+disp(sprintf('time tag is days/seconds : %d %d',days,seconds))
+
+offset = prod(size(S)); ind1 = 1; ind2 = offset;
+myS = reshape(datmat(ind1:ind2),size(S));
+d = myS - S;
+disp(sprintf('S diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(T)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myT = reshape(datmat(ind1:ind2),size(T));
+d = myT - T;
+disp(sprintf('T diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(U)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myU = reshape(datmat(ind1:ind2),size(U));
+d = myU - U;
+disp(sprintf('U diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(V)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+myV = reshape(datmat(ind1:ind2),size(V));
+d = myV - V;
+disp(sprintf('V diffs are %f %f',min(d(:)),max(d(:))))
+
+offset = prod(size(SSH)); ind1 = ind2+1; ind2 = ind1 + offset - 1;
+mySSH = reshape(datmat(ind1:ind2),size(SSH));
+d = mySSH - SSH;
+disp(sprintf('SSH diffs are %f %f',min(d(:)),max(d(:))))
+
+clear datmat myS myT myU myV mySSH d
+
+%----------------------------------------------------------------------
+% The perfect_ics file contains precisely two records.
+% record 1 is two integers defining the valid time of the state vector.
+% record 2 is the state vector.
+% This is exactly the same as the 'assim_model_state_ud' file.
+% (except they normally are for different times ...)
+%----------------------------------------------------------------------
+
+tbase = datenum(1601,1,1,0,0,0); % this is zero in the DART(gregorian) world
+t_one = datenum(1996,1,1,0,0,0); % valid time of the 'gom' files.
+
+toffset = t_one - tbase;
+days = floor(toffset);
+seconds = round(toffset - days)*86400;
+
+disp(sprintf('Creating a ''perfect_ics'' file with elements of shape %d %d %d',nx,ny,nz))
+
+disp(sprintf('prod(size(S)) is %d', prod(size(S))))
+disp(sprintf('prod(size(T)) is %d', prod(size(T))))
+disp(sprintf('prod(size(U)) is %d', prod(size(U))))
+disp(sprintf('prod(size(V)) is %d', prod(size(V))))
+disp(sprintf('prod(size(SSH)) is %d',prod(size(SSH))))
+
+nitems = prod(size(S))+ prod(size(T)) + prod(size(U)) + prod(size(V)) + prod(size(SSH));
+disp(sprintf('total restart size should %d',nitems*8 + 8 + 16))
+
+fid = fopen('gom_S_199601.bin','rb','ieee-be');
+[Sics,count] = fread(fid,prod(size(S)),'float32');
+fclose(fid);
+if (count ~= prod(size(S)))
+ error(sprintf('S record length wrong %d %d',count,prod(size(S))))
+end
+
+fid = fopen('gom_T_199601.bin','rb','ieee-be');
+[Tics,count] = fread(fid,prod(size(T)),'float32');
+fclose(fid);
+if (count ~= prod(size(T)))
+ error(sprintf('T record length wrong %d %d',count,prod(size(T))))
+end
+
+fid = fopen('gom_U_199601.bin','rb','ieee-be');
+[Uics,count] = fread(fid,prod(size(U)),'float32');
+fclose(fid);
+if (count ~= prod(size(U)))
+ error(sprintf('U record length wrong %d %d',count,prod(size(U))))
+end
+
+fid = fopen('gom_V_199601.bin','rb','ieee-be');
+[Vics,count] = fread(fid,prod(size(V)),'float32');
+fclose(fid);
+if (count ~= prod(size(V)))
+ error(sprintf('V record length wrong %d %d',count,prod(size(V))))
+end
+
+fid = fopen('gom_H_199601.bin','rb','ieee-be');
+[SSHics,count] = fread(fid,prod(size(SSH)),'float32');
+fclose(fid);
+if (count ~= prod(size(SSH)))
+ error(sprintf('SSH record length wrong %d %d',count,prod(size(SSH))))
+end
+
+datvec = [Sics; Tics; Uics; Vics; SSHics];
+disp(sprintf('total model size is %d',length(datvec)*8))
+
+fid = fopen('perfect_ics','wb','ieee-be');
+fwrite(fid, trec1,'int32');
+fwrite(fid,seconds,'int32');
+fwrite(fid, days,'int32');
+fwrite(fid, trecN,'int32');
+
+fwrite(fid, rec1,'int32');
+fwrite(fid, datvec,'float64');
+fwrite(fid, recN,'int32');
+fclose(fid);
+
+fid = fopen('perfect_ics.txt','wt');
+fprintf(fid,'%d %d\n',seconds,days);
+fprintf(fid,'%.15e\n',datvec);
+fclose(fid);
+
Property changes on: DART/trunk/models/MITgcm_ocean/matlab/Check_trans_pv_sv.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Added: DART/trunk/models/MITgcm_ocean/matlab/Check_trans_sv_pv.m
===================================================================
--- DART/trunk/models/MITgcm_ocean/matlab/Check_trans_sv_pv.m (rev 0)
+++ DART/trunk/models/MITgcm_ocean/matlab/Check_trans_sv_pv.m 2008-07-28 22:58:45 UTC (rev 3477)
@@ -0,0 +1,60 @@
+function Check_trans_sv_pv
+%
+%
+%time for timestep 0 day=128565, sec=21600
+%time for timestep 0 1953 Jan 01 06:00:00
+%time for 0000001344 day=128992, sec=21600
+%time for 0000001344 1954 Mar 04 06:00:00
+%
+
+dt = 900;
+startDate_1=19960101;
+startDate_2=60000;
+tbase = datenum(1996,1,1,6,0,0);
+
+S = rdmds( '~/MITgcm/data40levels/S.0000001344');
+T = rdmds( '~/MITgcm/data40levels/T.0000001344');
+U = rdmds( '~/MITgcm/data40levels/U.0000001344');
+V = rdmds( '~/MITgcm/data40levels/V.0000001344');
+SSH = rdmds('~/MITgcm/data40levels/Eta.0000001344');
+
+modelsize = prod(size(S)) + prod(size(T)) + ...
+ prod(size(U)) + prod(size(V)) + prod(size(SSH));
+
+toffset = (01344*dt)/86400
+
+[nx ny nz] = size(S);
+
+fid = fopen('assim_model_state_ud','rb','ieee-be');
+
+rec1 = fread(fid,1,'int32');
+seconds = fread(fid,1,'int32');
+days = fread(fid,1,'int32');
+recN = fread(fid,1,'int32');
+
+if (rec1 ~= recN)
+ error('first record mismatch')
+end
+disp(sprintf('time tag is days/seconds : %d %d',days,seconds))
+
+recN = fread(fid, 1,'int32');
+datmat = fread(fid,modelsize,'float64');
+recN = fread(fid, 1,'int32');
+fclose(fid);
+
+fid = fopen('assim_model_state_ic','wb','ieee-be');
+fwrite(fid,rec1,'int32');
+fwrite(fid,seconds+900,'int32');
+fwrite(fid,days+2,'int32');
+fwrite(fid,rec1,'int32');
+
+fwrite(fid,rec1,'int32');
+fwrite(fid,seconds,'int32');
+fwrite(fid,days,'int32');
+fwrite(fid,rec1,'int32');
+
+fwrite(fid, recN,'int32');
+fwrite(fid,datmat,'float64');
+fwrite(fid, recN,'int32');
+fclose(fid);
+
Property changes on: DART/trunk/models/MITgcm_ocean/matlab/Check_trans_sv_pv.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Modified: DART/trunk/models/MITgcm_ocean/shell_scripts/MakeInitialEnsemble.csh
===================================================================
--- DART/trunk/models/MITgcm_ocean/shell_scripts/MakeInitialEnsemble.csh 2008-07-28 17:50:02 UTC (rev 3476)
+++ DART/trunk/models/MITgcm_ocean/shell_scripts/MakeInitialEnsemble.csh 2008-07-28 22:58:45 UTC (rev 3477)
@@ -8,12 +8,14 @@
# $Id$
#
#=============================================================================
-# So the background is that we ran the model for 14 days with a 900 second
-# timestep - starting 1996 01 01 ??Z
+# So the background is that I want 40 snapshot files to create an initial
+# conditions file with 40 ensemble members.
#
-# I wanted 20 ensemble members.
+# For lack of a better idea: I ran the model starting from 1996 01 01 ??Z
+# for 20 days with a 900 second timestep and a dumpFreq of 12 hours.
+# This actually generates 41 snapshot files, but who's counting.
+# I just deleted the LAST snapshot before running this script.
#
-# I output snapshot files every 12 hours to generate 29 sets of snapshots.
# We're just going to make the first snapshot into ensemble member 1.
# Second snapshot will be ensemble member 2 ... I'm just not going to use
# the last snapshots. This may not be the greatest idea - I don't know
@@ -32,7 +34,7 @@
# to be whatever is in the namelist. In this case, we want the time
# tag to be 1996 01 01 00Z ... in DART-speak ... 144270 days 0 seconds.
#
-# TJH - 21 May 2008
+# TJH - 28 July 2008
# WARNING: You should run this in an 'empty' directory
# WARNING: You should run this in an 'empty' directory
@@ -40,7 +42,7 @@
# Get input for this execution
-set SNAPSHOTDIR = /fs/image/home/thoar/SVN/DART/models/MITgcm_ocean/run1
+set SNAPSHOTDIR = /ptmp/thoar/MITgcm/nodart/fortnight
set DARTEXE = /fs/image/home/thoar/SVN/DART/models/MITgcm_ocean/work
set DARTNML = /fs/image/home/thoar/SVN/DART/models/MITgcm_ocean/work/input.nml
@@ -80,6 +82,7 @@
echo "Using DART trans_pv_sv from ${DARTEXE}"
else
echo "ERROR: ${DARTEXE}/trans_pv_sv does not exist."
+ echo " Build ${DARTEXE}/trans_pv_sv"
set checkstat = 1
endif
@@ -87,6 +90,7 @@
echo "Using DART restart_file_utility from ${DARTEXE}"
else
echo "ERROR: ${DARTEXE}/restart_file_utility does not exist."
+ echo " Build ${DARTEXE}/restart_file_utility"
set checkstat = 1
endif
@@ -236,5 +240,3 @@
#-------------------------------------------------------------------------
\rm -f ens_mem_*
-
-
More information about the Dart-dev
mailing list