matlab_class/gcmfaces_diags/diags_grid_parms.m

function []=diags_grid_parms(dirModel,listTimes,doInteractive);
%object :      load grid, set params, and save both to dirMat
%input :       listTimes is the time list obtained from diags_list_times
%(optional)    doInteractive=1 allows users to specify parameters interactively
%                     doInteractive = 0 (default) uses ECCO v4 parameters
%                     and omits budgets and model-data misfits analyses

%global variables
gcmfaces_global;
global myparms;

%load grid
if isempty(dir('GRID'))&isempty(dir('nctiles_grid'))&...
   isempty(dir([dirModel filesep 'GRID']))&...
   isempty(dir([dirModel filesep 'nctiles_grid']));
    dirGrid=input('grid directory?\n');
    fprintf('\nFor the ECCO v4 LLC90 grid, the following parameters\n');
    fprintf('apply: nF=5; frmt=''compact''; memoryLimit=0; \n\n');
    nF=input('Number of faces? (nF=1, 4, 5 or 6)\n');
    frmt=input('File format? (frmt=''straight'', ''cube'' or ''compact'')\n');
    memoryLimit=input('memoryLimit? (0=load full grid, 1=less, 2=even less)\n');
elseif ~isempty(dir('GRID'));
    dirGrid=['GRID' filesep];
    nF=5;
    frmt='compact';
    memoryLimit=0;
elseif ~isempty(dir('nctiles_grid'));
    dirGrid=['nctiles_grid' filesep];
    nF=5;
    frmt='nctiles';
    memoryLimit=0;
elseif ~isempty(dir([dirModel filesep 'GRID']));
    dirGrid=[dirModel filesep 'GRID' filesep];
    nF=5;
    frmt='compact';
    memoryLimit=0;
elseif ~isempty(dir([dirModel filesep 'nctiles_grid']));
    dirGrid=[dirModel filesep 'nctiles_grid' filesep];
    nF=5;
    frmt='nctiles';
    memoryLimit=0;
end;

if doInteractive;
    nF=input('number of faces? (1, 4, 5or 6)\n'); 
    frmt=input('file format ? (''straight'', ''cube'' or ''compact'')\n');
    memoryLimit=input('memoryLimit ? (0=load full grid, 1=load less, 2=load even less)\n');
end;

grid_load(dirGrid,nF,frmt,memoryLimit);

%add definition of zonal and transport lines to mygrid
gcmfaces_lines_zonal;
[lonPairs,latPairs,names]=gcmfaces_lines_pairs;
gcmfaces_lines_transp(lonPairs,latPairs,names);

%set default for model run parameters
if doInteractive;
    choiceParams=input(['choice of default parameters? (1=ecco v4, ' ...
                 '2=core, 3=ecco2 adjoint, 4=ecco v3, 5=core-Jeff)\n']);
else;
    choiceParams=1;
end;
myparms=default_parms(choiceParams);

%allow user to change model params if necessary
myparms=review_parms(myparms,listTimes,doInteractive);

function [parms]=default_parms(choiceParams);
%set model parameters to default (ecco_v4)

if choiceParams==1|choiceParams==2|choiceParams==4|choiceParams==5;
parms.yearFirst=1992; %first year covered by model integration
parms.yearLast =2011; %last year covered by model integration
parms.yearInAve=[parms.yearFirst parms.yearLast]; %period for time averages and variance computations
parms.timeStep =3600; %model time step for tracers
parms.iceModel =1;%0=use freezing point   1=use pkg/seaice   2=use pkg/thsice
parms.useRFWF  =1;%1=real fresh water flux 0=virtual salt flux
parms.useNLFS  =2;%2=rstar 1=nlfs 0=linear free surface
parms.rhoconst =1029; %sea water density
parms.rcp      =3994*parms.rhoconst; % sea water rho X heat capacity
parms.rhoi     = 910; %sea ice density
parms.rhosn    = 330; %snow density
parms.flami    = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
parms.flamb    = 2.50e06; % latent heat of evaporation (J/kg)
parms.SIsal0   =4;
if choiceParams==2;
  parms.yearFirst=1948; %first year covered by model integration
  parms.yearLast =2007; %last year covered by model integration
  parms.yearInAve = [1948 2007];
end;
if choiceParams==4;
  parms.useRFWF  =0;%1=real fresh water flux 0=virtual salt flux
  parms.useNLFS  =0;%2=rstar 1=nlfs 0=linear free surface
end;
if choiceParams==5;
 parms.yearFirst=2006; %first year covered by model integration
 parms.yearLast =2305; %last year covered by model integration
 parms.yearInAve = [2006 2305];
end;
end;

if choiceParams==3;
        parms.yearFirst=2004; %first year covered by model integration
        parms.yearLast =2005; %last year covered by model integration
        parms.yearInAve=[parms.yearFirst parms.yearLast]; %period for time averages and variance computations
        parms.timeStep =1200; %model time step for tracers
        parms.iceModel =1;%0=use freezing point   1=use pkg/seaice   2=use pkg/thsice
        parms.useRFWF  =0;%1=real fresh water flux 0=virtual salt flux
        parms.useNLFS  =0;%2=rstar 1=nlfs 0=linear free surface
        parms.rhoconst =1027.5; %sea water density
        parms.rcp      =3994*parms.rhoconst; % sea water rho X heat capacity
        parms.rhoi     = 910; %sea ice density
        parms.rhosn    = 330; %snow density
        parms.flami    = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
        parms.flamb    = 2.50e06; % latent heat of evaporation (J/kg)
        parms.SIsal0   = 0;
end;

function [parms]=review_parms(parms,listTimes,doInteractive);
%review model parameters, correct them if needed, and check a couple more things

test1=1;%so that we print params at least once
while test1;
    fprintf('\n\n');
    gcmfaces_msg('model parameters summary','==== ');
    
    tmp1=sprintf('parms.yearFirst  = %i (first year covered by model integration)',parms.yearFirst); gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.yearLast   = %i (first year covered by model integration)',parms.yearLast);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.yearInAve  = [%i %i] (time mean and variance years)',parms.yearInAve);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.timeStep   = %i (model time step for tracers)',parms.timeStep);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.iceModel   = %i (0=freezing point  1=pkg/seaice  2=pkg/thsice)',parms.iceModel);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.useRFWF    = %i (1=real fresh water flux 0=virtual salt flux)',parms.useRFWF);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.useNLFS    = %i; (2=rstar 1=nlfs 0=linear free surface)',parms.useNLFS);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.rhoconst   = %0.6g (sea water density)',parms.rhoconst);  gcmfaces_msg(tmp1,'== ');
    tmp1=sprintf('parms.rcp        = %0.6g (sea water density X heat capacity)',parms.rcp);  gcmfaces_msg(tmp1,'== ');
    if parms.iceModel==1;
        tmp1=sprintf('parms.rhoi       = %0.6g (sea ice density)',parms.rhoi);  gcmfaces_msg(tmp1,'== ');
        tmp1=sprintf('parms.rhosn      = %0.6g (snow density)',parms.rhosn);  gcmfaces_msg(tmp1,'== ');
        tmp1=sprintf('parms.flami      = %0.6g (latent heat of fusion of ice/snow in J/kg)',parms.flami);  gcmfaces_msg(tmp1,'== ');
        tmp1=sprintf('parms.flamb      = %0.6g (latent heat of evaporation in J/kg)',parms.flamb);  gcmfaces_msg(tmp1,'== ');
        tmp1=sprintf('parms.SIsal0     = %0.6g (sea ice constant salinity)',parms.SIsal0);  gcmfaces_msg(tmp1,'== ');
        %tmp1=sprintf('',);  gcmfaces_msg(tmp1,'== ');
    else;
        error('only parms.iceModel=1 is currently treated\n');
    end;
    
    if doInteractive;
      gcmfaces_msg('to change a param type e.g. ''parms.yearFirst=1;'' or hit return if all params are ok. Change a param?','==== ');
      tmp1=input('');
      test1=~isempty(tmp1); %so that we change param and iterate process
      if test1; eval(tmp1); end;
    else;
      test1=[];
    end;    
end;

%determine a few more things about the diagnostic time axis
fprintf('\n\n');
parms.diagsNbRec=length(listTimes);
test1=median(diff(listTimes)*parms.timeStep/86400);
if abs(test1-30.5)<1; parms.diagsAreMonthly=1; else; parms.diagsAreMonthly=0; end;
if abs(test1-365.25)<1; parms.diagsAreAnnual=1; else; parms.diagsAreAnnual=0; end;
if doInteractive;
   tmp1=sprintf('parms.diagsNbRec       = %i (number of records, based on model output files)',parms.diagsNbRec); gcmfaces_msg(tmp1,'== ');
   tmp1=sprintf('parms.diagsAreMonthly  = %i (0/1 = false/true; based on output frequency)',parms.diagsAreMonthly); gcmfaces_msg(tmp1,'== ');
   tmp1=sprintf('parms.diagsAreAnnual   = %i (0/1 = false/true; based on output frequency)',parms.diagsAreAnnual); gcmfaces_msg(tmp1,'== ');
   gcmfaces_msg('hit return if this seems correct otherwise stop here','== '); test0=input(''); if ~isempty(test0); error('likely dir problem'); end;
end;

listTimes2=parms.yearFirst+listTimes*parms.timeStep/86400/365.25;%this approximation of course makes things simpler
tmp1=-0.5*diff(listTimes,1,1)*parms.timeStep/86400/365.25; tmp1=[median(tmp1);tmp1];
listTimes2=listTimes2+tmp1;%this converts the enddate to the middate of pkg/diags
ii=find(listTimes2>=parms.yearInAve(1)&listTimes2<=parms.yearInAve(2)+1);
if parms.diagsAreMonthly;%then restrict to full years
    ni=floor(length(ii)/12)*12; 
    if ni>0; 
      parms.recInAve=[ii(1) ii(floor(ni))];
    else;
      parms.recInAve=[ii(1) ii(end)];
    end;
elseif ~isempty(ii);
    parms.recInAve=[ii(1) ii(end)];
else;
    parms.recInAve=[1 1];
end;
if doInteractive; 
  tmp1=sprintf('parms.recInAve  = [%i %i] (time mean and variance records)',parms.recInAve);  gcmfaces_msg(tmp1,'== ');
end;

fprintf('\n\n');

1	gforget	1.12	function []=diags_grid_parms(dirModel,listTimes,doInteractive);
2	gforget	1.8	%object : load grid, set params, and save both to dirMat
3			%input : listTimes is the time list obtained from diags_list_times
4			%(optional) doInteractive=1 allows users to specify parameters interactively
5			% doInteractive = 0 (default) uses ECCO v4 parameters
6			% and omits budgets and model-data misfits analyses
7	gforget	1.1
8			%global variables
9			gcmfaces_global;
10			global myparms;
11
12			%load grid
13	gforget	1.12	if isempty(dir('GRID'))&isempty(dir('nctiles_grid'))&...
14			isempty(dir([dirModel filesep 'GRID']))&...
15			isempty(dir([dirModel filesep 'nctiles_grid']));
16	gforget	1.1	dirGrid=input('grid directory?\n');
17	gforget	1.13	fprintf('\nFor the ECCO v4 LLC90 grid, the following parameters\n');
18			fprintf('apply: nF=5; frmt=''compact''; memoryLimit=0; \n\n');
19			nF=input('Number of faces? (nF=1, 4, 5 or 6)\n');
20			frmt=input('File format? (frmt=''straight'', ''cube'' or ''compact'')\n');
21			memoryLimit=input('memoryLimit? (0=load full grid, 1=less, 2=even less)\n');
22	gforget	1.10	elseif ~isempty(dir('GRID'));
23	gforget	1.11	dirGrid=['GRID' filesep];
24	gforget	1.10	nF=5;
25			frmt='compact';
26			memoryLimit=0;
27			elseif ~isempty(dir('nctiles_grid'));
28	gforget	1.11	dirGrid=['nctiles_grid' filesep];
29	gforget	1.10	nF=5;
30			frmt='nctiles';
31			memoryLimit=0;
32	gforget	1.12	elseif ~isempty(dir([dirModel filesep 'GRID']));
33			dirGrid=[dirModel filesep 'GRID' filesep];
34			nF=5;
35			frmt='compact';
36			memoryLimit=0;
37			elseif ~isempty(dir([dirModel filesep 'nctiles_grid']));
38			dirGrid=[dirModel filesep 'nctiles_grid' filesep];
39			nF=5;
40			frmt='nctiles';
41			memoryLimit=0;
42	gforget	1.1	end;
43	gforget	1.2
44	gforget	1.8	if doInteractive;
45			nF=input('number of faces? (1, 4, 5or 6)\n');
46			frmt=input('file format ? (''straight'', ''cube'' or ''compact'')\n');
47			memoryLimit=input('memoryLimit ? (0=load full grid, 1=load less, 2=load even less)\n');
48			end;
49	gforget	1.2
50			grid_load(dirGrid,nF,frmt,memoryLimit);
51	gforget	1.1
52	gforget	1.14	%add definition of zonal and transport lines to mygrid
53			gcmfaces_lines_zonal;
54			[lonPairs,latPairs,names]=gcmfaces_lines_pairs;
55			gcmfaces_lines_transp(lonPairs,latPairs,names);
56
57	gforget	1.1	%set default for model run parameters
58	gforget	1.8	if doInteractive;
59			choiceParams=input(['choice of default parameters? (1=ecco v4, ' ...
60			'2=core, 3=ecco2 adjoint, 4=ecco v3, 5=core-Jeff)\n']);
61			else;
62			choiceParams=1;
63			end;
64	gforget	1.2	myparms=default_parms(choiceParams);
65	gforget	1.1
66			%allow user to change model params if necessary
67	gforget	1.8	myparms=review_parms(myparms,listTimes,doInteractive);
68	gforget	1.1
69	gforget	1.2	function [parms]=default_parms(choiceParams);
70	gforget	1.1	%set model parameters to default (ecco_v4)
71
72	gforget	1.6	if choiceParams==1\|choiceParams==2\|choiceParams==4\|choiceParams==5;
73	gforget	1.3	parms.yearFirst=1992; %first year covered by model integration
74	gforget	1.5	parms.yearLast =2011; %last year covered by model integration
75	gforget	1.1	parms.yearInAve=[parms.yearFirst parms.yearLast]; %period for time averages and variance computations
76	gforget	1.3	parms.timeStep =3600; %model time step for tracers
77	gforget	1.1	parms.iceModel =1;%0=use freezing point 1=use pkg/seaice 2=use pkg/thsice
78	gforget	1.3	parms.useRFWF =1;%1=real fresh water flux 0=virtual salt flux
79			parms.useNLFS =2;%2=rstar 1=nlfs 0=linear free surface
80			parms.rhoconst =1029; %sea water density
81	gforget	1.1	parms.rcp =3994*parms.rhoconst; % sea water rho X heat capacity
82			parms.rhoi = 910; %sea ice density
83			parms.rhosn = 330; %snow density
84			parms.flami = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
85			parms.flamb = 2.50e06; % latent heat of evaporation (J/kg)
86	gforget	1.3	parms.SIsal0 =4;
87	gforget	1.2	if choiceParams==2;
88			parms.yearFirst=1948; %first year covered by model integration
89			parms.yearLast =2007; %last year covered by model integration
90	gforget	1.7	parms.yearInAve = [1948 2007];
91	gforget	1.2	end;
92	gforget	1.4	if choiceParams==4;
93			parms.useRFWF =0;%1=real fresh water flux 0=virtual salt flux
94			parms.useNLFS =0;%2=rstar 1=nlfs 0=linear free surface
95			end;
96	gforget	1.6	if choiceParams==5;
97			parms.yearFirst=2006; %first year covered by model integration
98			parms.yearLast =2305; %last year covered by model integration
99			parms.yearInAve = [2006 2305];
100			end;
101	gforget	1.2	end;
102
103			if choiceParams==3;
104			parms.yearFirst=2004; %first year covered by model integration
105			parms.yearLast =2005; %last year covered by model integration
106			parms.yearInAve=[parms.yearFirst parms.yearLast]; %period for time averages and variance computations
107			parms.timeStep =1200; %model time step for tracers
108			parms.iceModel =1;%0=use freezing point 1=use pkg/seaice 2=use pkg/thsice
109			parms.useRFWF =0;%1=real fresh water flux 0=virtual salt flux
110			parms.useNLFS =0;%2=rstar 1=nlfs 0=linear free surface
111			parms.rhoconst =1027.5; %sea water density
112			parms.rcp =3994*parms.rhoconst; % sea water rho X heat capacity
113			parms.rhoi = 910; %sea ice density
114			parms.rhosn = 330; %snow density
115			parms.flami = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
116			parms.flamb = 2.50e06; % latent heat of evaporation (J/kg)
117			parms.SIsal0 = 0;
118			end;
119	gforget	1.1
120	gforget	1.8	function [parms]=review_parms(parms,listTimes,doInteractive);
121	gforget	1.1	%review model parameters, correct them if needed, and check a couple more things
122
123			test1=1;%so that we print params at least once
124			while test1;
125			fprintf('\n\n');
126			gcmfaces_msg('model parameters summary','==== ');
127
128			tmp1=sprintf('parms.yearFirst = %i (first year covered by model integration)',parms.yearFirst); gcmfaces_msg(tmp1,'== ');
129			tmp1=sprintf('parms.yearLast = %i (first year covered by model integration)',parms.yearLast); gcmfaces_msg(tmp1,'== ');
130			tmp1=sprintf('parms.yearInAve = [%i %i] (time mean and variance years)',parms.yearInAve); gcmfaces_msg(tmp1,'== ');
131			tmp1=sprintf('parms.timeStep = %i (model time step for tracers)',parms.timeStep); gcmfaces_msg(tmp1,'== ');
132			tmp1=sprintf('parms.iceModel = %i (0=freezing point 1=pkg/seaice 2=pkg/thsice)',parms.iceModel); gcmfaces_msg(tmp1,'== ');
133			tmp1=sprintf('parms.useRFWF = %i (1=real fresh water flux 0=virtual salt flux)',parms.useRFWF); gcmfaces_msg(tmp1,'== ');
134			tmp1=sprintf('parms.useNLFS = %i; (2=rstar 1=nlfs 0=linear free surface)',parms.useNLFS); gcmfaces_msg(tmp1,'== ');
135			tmp1=sprintf('parms.rhoconst = %0.6g (sea water density)',parms.rhoconst); gcmfaces_msg(tmp1,'== ');
136			tmp1=sprintf('parms.rcp = %0.6g (sea water density X heat capacity)',parms.rcp); gcmfaces_msg(tmp1,'== ');
137			if parms.iceModel==1;
138			tmp1=sprintf('parms.rhoi = %0.6g (sea ice density)',parms.rhoi); gcmfaces_msg(tmp1,'== ');
139			tmp1=sprintf('parms.rhosn = %0.6g (snow density)',parms.rhosn); gcmfaces_msg(tmp1,'== ');
140			tmp1=sprintf('parms.flami = %0.6g (latent heat of fusion of ice/snow in J/kg)',parms.flami); gcmfaces_msg(tmp1,'== ');
141			tmp1=sprintf('parms.flamb = %0.6g (latent heat of evaporation in J/kg)',parms.flamb); gcmfaces_msg(tmp1,'== ');
142			tmp1=sprintf('parms.SIsal0 = %0.6g (sea ice constant salinity)',parms.SIsal0); gcmfaces_msg(tmp1,'== ');
143			%tmp1=sprintf('',); gcmfaces_msg(tmp1,'== ');
144			else;
145			error('only parms.iceModel=1 is currently treated\n');
146			end;
147
148	gforget	1.8	if doInteractive;
149			gcmfaces_msg('to change a param type e.g. ''parms.yearFirst=1;'' or hit return if all params are ok. Change a param?','==== ');
150			tmp1=input('');
151			test1=~isempty(tmp1); %so that we change param and iterate process
152			if test1; eval(tmp1); end;
153			else;
154			test1=[];
155			end;
156	gforget	1.1	end;
157
158			%determine a few more things about the diagnostic time axis
159			fprintf('\n\n');
160			parms.diagsNbRec=length(listTimes);
161			test1=median(diff(listTimes)*parms.timeStep/86400);
162			if abs(test1-30.5)<1; parms.diagsAreMonthly=1; else; parms.diagsAreMonthly=0; end;
163			if abs(test1-365.25)<1; parms.diagsAreAnnual=1; else; parms.diagsAreAnnual=0; end;
164	gforget	1.8	if doInteractive;
165	gforget	1.9	tmp1=sprintf('parms.diagsNbRec = %i (number of records, based on model output files)',parms.diagsNbRec); gcmfaces_msg(tmp1,'== ');
166			tmp1=sprintf('parms.diagsAreMonthly = %i (0/1 = false/true; based on output frequency)',parms.diagsAreMonthly); gcmfaces_msg(tmp1,'== ');
167			tmp1=sprintf('parms.diagsAreAnnual = %i (0/1 = false/true; based on output frequency)',parms.diagsAreAnnual); gcmfaces_msg(tmp1,'== ');
168	gforget	1.8	gcmfaces_msg('hit return if this seems correct otherwise stop here','== '); test0=input(''); if ~isempty(test0); error('likely dir problem'); end;
169			end;
170	gforget	1.1
171			listTimes2=parms.yearFirst+listTimes*parms.timeStep/86400/365.25;%this approximation of course makes things simpler
172			tmp1=-0.5diff(listTimes,1,1)parms.timeStep/86400/365.25; tmp1=[median(tmp1);tmp1];
173			listTimes2=listTimes2+tmp1;%this converts the enddate to the middate of pkg/diags
174			ii=find(listTimes2>=parms.yearInAve(1)&listTimes2<=parms.yearInAve(2)+1);
175			if parms.diagsAreMonthly;%then restrict to full years
176	gforget	1.6	ni=floor(length(ii)/12)*12;
177			if ni>0;
178			parms.recInAve=[ii(1) ii(floor(ni))];
179			else;
180			parms.recInAve=[ii(1) ii(end)];
181			end;
182	gforget	1.1	elseif ~isempty(ii);
183			parms.recInAve=[ii(1) ii(end)];
184			else;
185			parms.recInAve=[1 1];
186			end;
187	gforget	1.9	if doInteractive;
188			tmp1=sprintf('parms.recInAve = [%i %i] (time mean and variance records)',parms.recInAve); gcmfaces_msg(tmp1,'== ');
189			end;
190	gforget	1.1
191			fprintf('\n\n');
192