[Dart-dev] [7901] DART/trunk/models/cam: Documentation for the latest version of the CAM
nancy at ucar.edu
nancy at ucar.edu
Fri Apr 17 14:10:18 MDT 2015
Revision: 7901
Author: nancy
Date: 2015-04-17 14:10:18 -0600 (Fri, 17 Apr 2015)
Log Message:
-----------
Documentation for the latest version of the CAM
model_mod, plus documentation (and matlab scripts
to generate the figures) that explain how the
assimilation avoids problems at the top of the model
by not updating the highest levels of the model state
since it is so highly damped.
Modified Paths:
--------------
DART/trunk/models/cam/model_mod.html
Added Paths:
-----------
DART/trunk/models/cam/doc/highest_state_p_Pa.pdf
DART/trunk/models/cam/matlab/highest_state_p_Pa.m
DART/trunk/models/cam/matlab/highest_state_scale_h.m
Property Changed:
----------------
DART/trunk/models/cam/doc/
-------------- next part --------------
Property changes on: DART/trunk/models/cam/doc
___________________________________________________________________
Added: svn:mergeinfo
+ /DART/branches/cam/models/cam/doc:6639-7899
/DART/branches/cam-update/doc:4903-4923
/DART/branches/development/models/cam/doc:4680-6255
/DART/branches/gitm/models/cam/doc:5143-6215
/DART/branches/gitm_lanai/models/cam/doc:6571-6652
/DART/branches/helen/models/cam/doc:5995-6161
Copied: DART/trunk/models/cam/doc/highest_state_p_Pa.pdf (from rev 7899, DART/branches/cam/models/cam/doc/top_of_model_damping.pdf)
===================================================================
(Binary files differ)
Added: DART/trunk/models/cam/matlab/highest_state_p_Pa.m
===================================================================
--- DART/trunk/models/cam/matlab/highest_state_p_Pa.m (rev 0)
+++ DART/trunk/models/cam/matlab/highest_state_p_Pa.m 2015-04-17 20:10:18 UTC (rev 7901)
@@ -0,0 +1,283 @@
+function highest_state_p_Pa(fname, div24del2flag, cutoff, vert_norm, varargin)
+% plots the vertical damping profile for CAM, WACCM
+%
+% highest_state_p_Pa(fname, div24del2flag, cutoff, vert_norm)
+% highest_state_p_Pa(fname, div24del2flag, cutoff, vert_norm, 'nu_top', nu_top)
+%
+% fname = pathname of initial files containing vertical grid information.
+%
+% div24del2flag = CAM namelist div24del2flag value;
+% 0 for CAM-SE,
+% 2 for CAM-FV div2 damping
+% 4 for CAM-FV div4 damping
+% 24 for CAM-FV div4+del2 damping
+%
+% cutoff = DART namelist variable 'cutoff'
+%
+% vert_norm = DART namelist variable 'vert_normalization_pressure'
+%
+% Optional arguments:
+% nu_top = needed only for CAM-SE (div24del2flag = 0).
+% Get from CAM's atm_in namelist.
+%
+% Plot CAM's (WACCM's) extra diffusion in the top layers
+% These are functions of CAM's ptop and -FV's div24del2flag.
+% Plot cam/model_mod.f90's highest_state_pressure_Pa and distance increment added to
+% the nominal distances between an ob and a state variable location,
+% which are functions of DART's cutoff and vert_normalization_pressure,
+% as well as CAM's ptop and div24del2flag.
+%
+% EXAMPLE:
+% fname = 'caminput.nc';
+% div24del2flag = 2;
+% cutoff = 0.2; % DART input.nml cutoff
+% vert_norm = 20000.0; % DART input.nml vert_normalization_pressure
+% highest_state_p_Pa(fname, div24del2flag, cutoff, vert_norm)
+%
+% EXAMPLE:
+% fname = 'caminput.nc';
+% div24del2flag = 0;
+% cutoff = 0.2; % DART input.nml cutoff
+% vert_norm = 20000.0; % DART input.nml vert_normalization_pressure
+% highest_state_p_Pa(fname, div24del2flag, cutoff, vert_norm,'nu_top',0.00005)
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+default_nu_top = 0.00001;
+args = inputParser;
+
+addRequired(args,'fname', at ischar);
+addRequired(args,'div24del2flag', at isnumeric);
+addRequired(args,'cutoff', at isnumeric);
+addRequired(args,'vert_norm', at isnumeric);
+addParamValue(args,'nu_top',default_nu_top, at isnumeric);
+parse(args,fname,div24del2flag,cutoff,vert_norm,varargin{:});
+
+if ~isempty(fieldnames(args.Unmatched))
+ disp('Extra inputs:')
+ disp(args.Unmatched)
+end
+
+% if you want to echo the input
+fprintf('fname : %s\n', args.Results.fname)
+fprintf('div24del2flag : %i\n', args.Results.div24del2flag)
+fprintf('cutoff : %f radians\n', args.Results.cutoff)
+fprintf('vert_normalization_pressure : %f Pa/radian\n', args.Results.vert_norm)
+if (div24del2flag == 0)
+ fprintf('nu_top : %f\n', args.Results.nu_top)
+end
+
+% Read in nlevs, hybrid coords (As and Bs) from the initial file.
+if (exist(fname,'file') ~= 2)
+ error('file/fname <%s> does not exist',fname)
+end
+
+if (div24del2flag == 4)
+ error('diffusion flavor 4 ("div4") is independent of height. Can set highest_state_pressure_Pa = 0')
+end
+
+As = local_nc_varget(fname, 'hyam');
+Bs = local_nc_varget(fname, 'hybm');
+fulls = local_nc_varget(fname, 'hyai');
+
+ncid = netcdf.open(fname,'NOWRITE');
+dimid = netcdf.inqDimID(ncid,'lev');
+[~, nlevs] = netcdf.inqDim(ncid, dimid);
+netcdf.close(ncid);
+
+% Convert cutoff from radians to Pa.
+% For illustration; DART does all of this in radian space.
+
+c_v = cutoff*vert_norm;
+fprintf('cutoff*vert_norm is %f\n',c_v)
+
+% Calculate p(k=1,nlevs) from As and Bs,
+
+ptop = 1.0e5 * fulls(1);
+p = 1.0e5 * (As + Bs);
+
+
+if (div24del2flag == 0)
+ % Calculate the "regular diffusion" in the top 3 layers for CAM-SE
+ diff_bot_k = 3;
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = zeros(nlevs,1);
+ tau(1) = args.Results.nu_top * 4;
+ tau(2) = args.Results.nu_top * 2;
+ tau(3) = args.Results.nu_top;
+
+elseif (div24del2flag == 2)
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = 8.0*(1.0 + tanh(log(ptop./p)));
+ % Calculate the div2 diffusion coefficient from p
+ % and the pressure level above which distance increments will be added.
+ tau_div2 = ones(nlevs,1);
+ for k = 1:nlevs
+ if (tau(k) < 1.0)
+ diff_bot_k = k;
+ break
+ else
+ tau_div2(k) = tau(k);
+ end
+ end
+elseif (div24del2flag == 24)
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = 8.0*(1.0 + tanh(log(ptop./p)));
+ % Calculate the del2 diffusion coefficient from p.
+ tau_del2 = zeros(nlevs,1);
+ for k = 1:nlevs
+ if (tau(k) < 0.3)
+ diff_bot_k = k;
+ break
+ else
+ tau_del2(k) = tau(k);
+ end
+ end
+else
+ error('div24del2flag not recognized. Choose from {0,2,24,4}')
+end
+
+highest_spPa = p(diff_bot_k) + c_v*(1.0 + sqrt(1.0 + 2*(p(diff_bot_k) - ptop)./(c_v)));
+
+% Distance increment in units of Pa (in model_mod it's radians)
+delta_d = zeros(nlevs,1);
+for k = 1:nlevs
+ if (p(k) < highest_spPa)
+ delta_d(k) = (highest_spPa - p(k))^2 ./ (highest_spPa - ptop);
+ else
+ % Save a k-range for plotting only the interesting pressure range.
+ %highest_spPa_k = k+3;
+ highest_spPa_k = nlevs;
+ break
+ end
+end
+
+clf; orient tall; wysiwyg
+
+% little summary of the input at the bottom
+axes('position',[0.1, 0.03, 0.8, 0.1125]);
+axis off
+axis([0 1 0 1.5])
+h1 = text(0.1, 0.00, sprintf('file providing model levels = %s',args.Results.fname));
+h2 = text(0.1, 0.25, sprintf('div24del2flag = %i',args.Results.div24del2flag));
+h3 = text(0.1, 0.50, sprintf('cutoff = %f radians',args.Results.cutoff));
+h4 = text(0.1, 0.75, sprintf('vert_normalization_pressure = %f Pa/radian',args.Results.vert_norm));
+if (div24del2flag == 0)
+ h5 = text(0.1, 1.00, sprintf('nu_top = %f',args.Results.nu_top));
+ set(h5,'Interpreter','none','FontName','Courier New','FontSize',12)
+end
+set(h1,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h2,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h3,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h4,'Interpreter','none','FontName','Courier New','FontSize',12)
+
+% Plot on a reversed log10 scale
+ax1 = axes('position',[0.125, 0.2, 0.75, 0.7]);
+set(ax1, 'Ydir','reverse', ...
+ 'Yscale','log', ...
+ 'Layer','top', ...
+ 'FontSize',20)
+
+hold on;
+ylabel('Pressure (Pa)')
+
+if (div24del2flag == 0)
+ % Plot (2) (div2_tau vs p(k)
+ h = plot(tau(1:highest_spPa_k), p(1:highest_spPa_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau)*1.05];
+elseif (div24del2flag == 2)
+ % Plot (2) (div2_tau vs p(k)
+ h = plot(tau_div2(1:highest_spPa_k), p(1:highest_spPa_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau_div2)*1.05];
+elseif (div24del2flag == 24)
+ h = plot(tau_del2(1:highest_spPa_k), p(1:highest_spPa_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau_del2)*1.05];
+end
+
+% A line at p(top)
+
+ptop_y = [ptop, ptop];
+ptop_line = line(ptop_x,ptop_y,'Color','k');
+set(ptop_line,'LineStyle','--','Marker','.', 'LineWidth',2.0);
+
+% A line at highest_state_pressure_Pa for default cutoff*vert_norm
+% Other cutoff*vert_norm in a lighter shade? NO; handle in plot of H_Pa vs c*v, below
+
+highest_spPa_y = [highest_spPa, highest_spPa];
+highest_spPa_line = line(ptop_x,highest_spPa_y,'Color','r');
+set(highest_spPa_line,'LineStyle',':','Marker','.', 'LineWidth',2.0);
+
+xlabel('diffusion coefficient \tau')
+
+% Second x-axis (on top) for plotting distance increments.
+
+ax2_XLim = [-.05*max(delta_d), 1.05*max(delta_d)];
+ax2 = axes('position',get(ax1,'Position'), ...
+ 'FontSize',20, ...
+ 'XAxisLocation','top', ...
+ 'YAxisLocation','right',...
+ 'Color','none',...
+ 'XColor','r','YColor','k',...
+ 'XLim',ax2_XLim, ...
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'), ...
+ 'Layer','top', ...
+ 'Yscale','log');
+% set(ax2, 'Ydir','reverse', 'Yscale','log', 'Layer','top')
+
+ylabel('Pressure (Pa)')
+xlabel('Extra distance used in cutoff calculation')
+
+% A curve of delta-d above highest_state_pressure_Pa for default cutoff*vert_norm (only)
+d_line = line(delta_d(1:highest_spPa_k),p(1:highest_spPa_k),'Color','r');
+set(d_line,'Marker','x', 'LineWidth',2.0);
+
+% A line at 2*cutoff
+cutoff_x = [2*c_v, 2*c_v];
+y_axis_height = get(ax2,'YLim');
+cutoff_y = [y_axis_height(1), p(nlevs-5)];
+cutoff_line = line(cutoff_x, cutoff_y, 'Color','r');
+set(cutoff_line,'LineStyle','--','Marker','.', 'LineWidth',2.0);
+
+% create a png file for import to powerpoint
+
+outfname = sprintf('tau_dist_c%3.2f_v%i.png', cutoff,vert_norm);
+print(gcf,'-dpng',outfname);
+
+% Plot (min) highest_state_pressure_Pa vs cutoff*vert_norm.
+
+legend_handles = [cutoff_line, d_line, highest_spPa_line, ptop_line, h];
+lh = legend(legend_handles,'obs impact limit','extra "distance"', ...
+ 'highest\_state\_pressure\_Pa','pressure at model top','CAM \tau');
+set(lh,'Interpreter','tex','Location','SouthEast','box','off')
+
+%=====================================================================
+
+
+function value = local_nc_varget(fname,varname)
+%% If the variable exists in the file, return the contents of the variable.
+% if the variable does not exist, return empty value instead of error-ing
+% out.
+
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname,'NOWRITE');
+ value = netcdf.getVar(ncid, varid);
+ netcdf.close(ncid)
+else
+ value = [];
+end
+
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Property changes on: DART/trunk/models/cam/matlab/highest_state_p_Pa.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/cam/matlab/highest_state_scale_h.m
===================================================================
--- DART/trunk/models/cam/matlab/highest_state_scale_h.m (rev 0)
+++ DART/trunk/models/cam/matlab/highest_state_scale_h.m 2015-04-17 20:10:18 UTC (rev 7901)
@@ -0,0 +1,281 @@
+function highest_state_scale_h(fname, div24del2flag, cutoff, vert_norm, varargin)
+% plots the vertical damping profile for CAM, WACCM
+%
+% highest_state_scale_h(fname, div24del2flag, cutoff, vert_norm)
+% highest_state_scale_h(fname, div24del2flag, cutoff, vert_norm, 'nu_top', nu_top)
+%
+% fname = pathname of initial files containing vertical grid information.
+%
+% div24del2flag = CAM namelist div24del2flag value;
+% 0 for CAM-SE,
+% 2 for CAM-FV div2 damping
+% 4 for CAM-FV div4 damping
+% 24 for CAM-FV div4+del2 damping
+%
+% cutoff = DART namelist variable 'cutoff'
+%
+% vert_norm = DART namelist variable 'vert_normalization_scale_height'
+%
+% Optional arguments:
+% nu_top = needed only for CAM-SE (div24del2flag = 0).
+% Get from CAM's atm_in namelist.
+%
+% Plot CAM's (WACCM's) extra diffusion in the top layers
+% These are functions of CAM's ptop and -FV's div24del2flag.
+% Plot cam/model_mod.f90's highest_state_pressure_Pa and distance increment added to
+% the nominal distances between an ob and a state variable location,
+% which are functions of DART's cutoff and vert_normalization_scale_height,
+% as well as CAM's ptop and div24del2flag.
+%
+% EXAMPLE:
+% fname = 'caminput.nc';
+% div24del2flag = 2;
+% cutoff = 0.2; % DART input.nml cutoff
+% vert_norm = 2.5; % DART input.nml vert_normalization_scale_height
+% highest_state_scale_h(fname, div24del2flag, cutoff, vert_norm)
+%
+% EXAMPLE:
+% fname = 'caminput.nc';
+% div24del2flag = 0;
+% cutoff = 0.2; % DART input.nml cutoff
+% vert_norm = 2.5; % DART input.nml vert_normalization_scale_height
+% highest_state_scale_h(fname, div24del2flag, cutoff, vert_norm,'nu_top',0.00005)
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+default_nu_top = 0.00001;
+args = inputParser;
+
+addRequired(args,'fname', at ischar);
+addRequired(args,'div24del2flag', at isnumeric);
+addRequired(args,'cutoff', at isnumeric);
+addRequired(args,'vert_norm', at isnumeric);
+addParamValue(args,'nu_top',default_nu_top, at isnumeric);
+parse(args,fname,div24del2flag,cutoff,vert_norm,varargin{:});
+
+if ~isempty(fieldnames(args.Unmatched))
+ disp('Extra inputs:')
+ disp(args.Unmatched)
+end
+
+% if you want to echo the input
+fprintf('fname : %s\n', args.Results.fname)
+fprintf('div24del2flag : %i\n', args.Results.div24del2flag)
+fprintf('cutoff : %f radians\n', args.Results.cutoff)
+fprintf('vert_normalization_scale_height : %f Pa/radian\n', args.Results.vert_norm)
+if (div24del2flag == 0)
+ fprintf('nu_top : %f\n', args.Results.nu_top)
+end
+
+% Read in nlevs, hybrid coords (As and Bs) from the initial file.
+if (exist(fname,'file') ~= 2)
+ error('file/fname <%s> does not exist',fname)
+end
+
+if (div24del2flag == 4)
+ error('diffusion flavor 4 ("div4") is independent of height. Can set highest_state_scale_height = 0')
+end
+
+As = local_nc_varget(fname, 'hyam');
+Bs = local_nc_varget(fname, 'hybm');
+fulls = local_nc_varget(fname, 'hyai');
+
+ncid = netcdf.open(fname,'NOWRITE');
+dimid = netcdf.inqDimID(ncid,'lev');
+[~, nlevs] = netcdf.inqDim(ncid, dimid);
+netcdf.close(ncid);
+
+% Convert cutoff from radians to Pa.
+% For illustration; DART does all of this in radian space.
+
+c_v = cutoff*vert_norm;
+fprintf('cutoff*vert_norm is %f\n',c_v)
+
+% Calculate scale heights(k=1,nlevs) from As and Bs,
+% Note that fulls(1) needs *10^5 to be a pressure,
+% but then that pressure is in the denominator under p_surf(=10^5)
+
+ptop = 1.0e5 * fulls(1);
+p = 1.0e5 * (As + Bs);
+ptop_sh = log(1.0 ./fulls(1) );
+height = log(1.0 ./(As + Bs));
+
+if (div24del2flag == 0)
+ % Calculate the "regular diffusion" in the top 3 layers for CAM-SE
+ diff_bot_k = 3;
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = zeros(nlevs,1);
+ tau(1) = args.Results.nu_top * 4;
+ tau(2) = args.Results.nu_top * 2;
+ tau(3) = args.Results.nu_top;
+
+elseif (div24del2flag == 2)
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = 8.0*(1.0 + tanh(log(ptop./p)));
+ % Calculate the div2 diffusion coefficient from p
+ % and the pressure level above which distance increments will be added.
+ tau_div2 = ones(nlevs,1);
+ for k = 1:nlevs
+ if (tau(k) < 1.0)
+ diff_bot_k = k;
+ break
+ else
+ tau_div2(k) = tau(k);
+ end
+ end
+elseif (div24del2flag == 24)
+ % Calculate the diffusion coefficient without height restrictions.
+ tau = 8.0*(1.0 + tanh(log(ptop./p)));
+ % Calculate the del2 diffusion coefficient from p.
+ tau_del2 = zeros(nlevs,1);
+ for k = 1:nlevs
+ if (tau(k) < 0.3)
+ diff_bot_k = k;
+ break
+ else
+ tau_del2(k) = tau(k);
+ end
+ end
+else
+ error('div24del2flag not recognized. Choose from {0,2,24,4}')
+end
+
+highest_sh = height(diff_bot_k) - c_v*(1.0 + sqrt(1.0 + 2*(ptop_sh - height(diff_bot_k))./(c_v)));
+
+% Distance increment in units of Pa (in model_mod it's radians)
+delta_d = zeros(nlevs,1);
+for k = 1:nlevs
+ if (height(k) >= highest_sh)
+ delta_d(k) = (highest_sh - height(k))^2 ./ (ptop_sh - highest_sh);
+ else
+ % Save a k-range for plotting only the interesting pressure range.
+ %highest_sh_k = k+3;
+ highest_sh_k = nlevs;
+ break
+ end
+end
+
+clf; orient tall; wysiwyg
+
+% little summary of the input at the bottom
+axes('position',[0.1, 0.03, 0.8, 0.1125]);
+axis off
+axis([0 1 0 1.5])
+h1 = text(0.1, 0.00, sprintf('file providing model levels = %s',args.Results.fname));
+h2 = text(0.1, 0.25, sprintf('div24del2flag = %i',args.Results.div24del2flag));
+h3 = text(0.1, 0.50, sprintf('cutoff = %f radians',args.Results.cutoff));
+h4 = text(0.1, 0.75, sprintf('vert_normalization_scale_height = %f scale_height/radian',args.Results.vert_norm));
+if (div24del2flag == 0)
+ h5 = text(0.1, 1.00, sprintf('nu_top = %f',args.Results.nu_top));
+ set(h5,'Interpreter','none','FontName','Courier New','FontSize',12)
+end
+set(h1,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h2,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h3,'Interpreter','none','FontName','Courier New','FontSize',12)
+set(h4,'Interpreter','none','FontName','Courier New','FontSize',12)
+
+% Plot on a linear scale.
+ax1 = axes('position',[0.125, 0.2, 0.75, 0.7],'FontSize',20);
+
+hold on;
+ylabel('Scale Height')
+
+if (div24del2flag == 0)
+ % Plot (2) (div2_tau vs p(k)
+ h = plot(tau(1:highest_sh_k), height(1:highest_sh_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau)*1.05];
+elseif (div24del2flag == 2)
+ % Plot (2) (div2_tau vs p(k)
+ h = plot(tau_div2(1:highest_sh_k), height(1:highest_sh_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau_div2)*1.05];
+elseif (div24del2flag == 24)
+ h = plot(tau_del2(1:highest_sh_k), height(1:highest_sh_k), 'k.-', ...
+ 'LineWidth',2.0, 'MarkerSize',20);
+ ptop_x = [0.0, max(tau_del2)*1.05];
+end
+
+% A line at p(top)
+
+ptop_y = [ptop_sh, ptop_sh];
+ptop_line = line(ptop_x,ptop_y,'Color','k');
+set(ptop_line,'LineStyle','--','Marker','.', 'LineWidth',2.0);
+
+% A line at highest_state_scale_height for default cutoff*vert_norm
+% Other cutoff*vert_norm in a lighter shade? NO; handle in plot of H_Pa vs c*v, below
+
+highest_sh_y = [highest_sh, highest_sh];
+highest_sh_line = line(ptop_x,highest_sh_y,'Color','r');
+set(highest_sh_line,'LineStyle',':','Marker','.', 'LineWidth',2.0);
+
+xlabel('diffusion coefficient \tau')
+
+% Second x-axis (on top) for plotting distance increments.
+
+ax2_XLim = [-.05*max(delta_d), 1.05*max(delta_d)];
+ax2 = axes('position',get(ax1,'Position'), ...
+ 'FontSize',20, ...
+ 'XAxisLocation','top', ...
+ 'YAxisLocation','right',...
+ 'Color','none',...
+ 'XColor','r','YColor','k',...
+ 'XLim',ax2_XLim, ...
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'), ...
+ 'Layer','top');
+% set(ax2, 'Ydir','reverse', 'Yscale','log', 'Layer','top')
+
+ylabel('Scale Height')
+xlabel('Extra distance used in cutoff calculation')
+
+% A curve of delta-d above highest_state_pressure_Pa for default cutoff*vert_norm (only)
+d_line = line(delta_d(1:highest_sh_k),height(1:highest_sh_k),'Color','r');
+set(d_line,'Marker','x', 'LineWidth',2.0);
+
+% A line at 2*cutoff
+cutoff_x = [2*c_v, 2*c_v];
+y_axis_height = get(ax2,'YLim');
+cutoff_y = [y_axis_height(2),height(nlevs-10)];
+cutoff_line = line(cutoff_x, cutoff_y, 'Color','r');
+set(cutoff_line,'LineStyle','--','Marker','.', 'LineWidth',2.0);
+
+% create a png file for import to powerpoint
+
+outfname = sprintf('tau_dist_c%3.2f_v%i.png', cutoff, vert_norm);
+print(gcf,'-dpng',outfname);
+
+% Plot (min) highest_state_scale_height vs cutoff*vert_norm.
+
+legend_handles = [cutoff_line, d_line, highest_sh_line, ptop_line, h];
+lh = legend(legend_handles,'obs impact limit','extra "distance"', ...
+ 'highest\_state\_scale\_height','scale height at model top','CAM \tau');
+set(lh,'Interpreter','tex','Location','SouthEast','box','off')
+
+%=====================================================================
+
+
+function value = local_nc_varget(fname,varname)
+%% If the variable exists in the file, return the contents of the variable.
+% if the variable does not exist, return empty value instead of error-ing
+% out.
+
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname,'NOWRITE');
+ value = netcdf.getVar(ncid, varid);
+ netcdf.close(ncid)
+else
+ value = [];
+end
+
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Property changes on: DART/trunk/models/cam/matlab/highest_state_scale_h.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/models/cam/model_mod.html
===================================================================
--- DART/trunk/models/cam/model_mod.html 2015-04-17 19:59:35 UTC (rev 7900)
+++ DART/trunk/models/cam/model_mod.html 2015-04-17 20:10:18 UTC (rev 7901)
@@ -27,11 +27,13 @@
<A HREF="#Namelist">NAMELIST</A> /
+<A HREF="#SetupVariations">SETUP VARIATIONS</A> /
<A HREF="#Interface">INTERFACES</A> /
<A HREF="#PrivateComponents">PRIVATE COMPONENTS</A> /
-<A HREF="#Discussion">Discussion</A> /
+<A HREF="#Discussion">DISCUSSION</A> /
<A HREF="#FilesUsed">FILES</A> /
<A HREF="#References">REFERENCES</A> /
+<A HREF="#Contributors">CONTRIBUTORS</A> /
<A HREF="#Errors">ERRORS</A> /
<A HREF="#KnownBugs">BUGS</A> /
<A HREF="#FuturePlans">PLANS</A> /
@@ -45,20 +47,44 @@
<P>
The DART system supports data assimilation into the Community Atmosphere Model, CAM,
which is the atmospheric component of the Community Earth System Model
-(<a href="http://www2.cesm.ucar.edu/models">CESM</a>) .
-The Whole Atmosphere Community Climate Model, WACCM, is a variant of CAM
-with a higher model top and more complete handling of chemistry.
-WACCM can also be used in DART using this model_mod.
-It is being used by graduate students,
-post-graduates, and scientists at universities and research labs.
-Others are using atmospheric analyses for their time period and resolution of interest,
-which were produced here at NCAR using CESM+DART.
+(<a href="http://www2.cesm.ucar.edu/models">CESM</a>).
+This DART interface is being used by graduate students,
+post-graduates, and scientists at universities and research labs
+to conduct data assimilation reseearch.
+Others are using the products of data assimilation (analyses),
+which were produced here at NCAR using CESM+DART,
+to conduct related research.
The variety of research can be sampled on the DART
-<a href:"http://www.image.ucar.edu/DAReS/Publications/index.php">Publications</a> page.
-In addition to the standard DART features described elsewhere,
-current capabilities include the ability to:
+<a href="http://www.image.ucar.edu/DAReS/Publications/index.php">Publications</a> page.
</P>
+<P>
+"CAM" refers to a family of related atmospheric components,
+which can be built with 2 independent main characteristics. CESM labels these as:
+<pre>
+ 'resolution' = horizontal (not vertical) grid AND dynamical core (fluid dynamics equations on that grid)
+ 3 supported dycores; eulerian, FV, and SE
+ SE refined grids will have some support, but by their nature invite the use
+ of user defined grid and map files.
+ 'compset' = vertical grid AND parameterizations (aka physics):
+ Parameterization is the equations describing a physical process like
+ convection, radiation, chemistry, ...
+ Vertical grid is determined by the needs of the chosen parameterizations.
+ Spacing and height of the top (ptop) vary.
+ The combinations of parameterizations and vertical grids are named:
+ CAM3.5, CAM5, CAM#, ...
+ WACCM, WACCM#, WACCM-X,
+ CAM-Chem,
+ ...
+</pre>
+There are minor characteristics choices within each of these,
+but only chemistry choices in WACCM and CAM-Chem have an impact on DART.
+As of April, 2015, all of these variants are handled by the same model_mod.f90,
+namelist, and build scripts, with differences in the assimilation set up
+described <a href="#SetupVariations">here</a>.
+</P>
+
+This DART+CAM interface has the following features.
<UL>
<LI>Assimilate within the CESM software framework by using the multi-instance
capability of CESM1.1.1 (and later).
@@ -67,7 +93,7 @@
when needed, is not being actively supported for these CESM versions.
<LI>Use either the eulerian, finite-volume (FV), or spectral-element (SE) dynamical core.</LI>
<LI>Use any resolution of CAM, including refined mesh grids in CAM-SE.
-As of April, 2014 this is limited by ability of the memory of a node of your hardware to contain
+As of April, 2015 this is limited by the ability of the memory of a node of your hardware to contain
the state vector of a single ensemble member. Work is under way to relax this restriction.
</LI>
<LI>Assimilate a variety of observations; to date the observations
@@ -82,7 +108,7 @@
This allows users to change the model state without recompiling
(but other restrictions remain).</LI>
<LI>Generate analyses on the CAM grid which have only CAM model error in them,
-rather than a combination of another model's error and observation error.</LI>
+rather than another model's.</LI>
<LI>Generate such analyses with as few as 20 ensemble members.</LI>
</UL>
@@ -92,13 +118,13 @@
<a href="http://www.image.ucar.edu/pub/DART/CAM/">
http://www.image.ucar.edu/pub/DART/CAM/</a>
that are helpful for interfacing CAM with DART.
-In the current (2014) mode, CESM+DART can easily be started from a single
+In the current (2015) mode, CESM+DART can easily be started from a single
model state, which is perturbed to create an ensemble of the desired size.
A spin-up period is then required to allow the ensemble members to diverge.
</P>
<P>
-Sample sets of observations, which can be used with DART-CAM assimilations,
+Sample sets of observations, which can be used with CESM+DART assimilations,
can be found at
<a href="http://www.image.ucar.edu/pub/DART/Obs_sets/">
http://www.image.ucar.edu/pub/DART/Obs_sets/</a>
@@ -128,7 +154,7 @@
Some information can be tracked down in the atm.log file, but not all of it.
<LI> The restart files (for non-chemistry model versions) are much larger than the
initial files (and we need to deal with an ensemble of them).
-<LI> The temperature on the restart files is virtual equivalent potential temperature (?),
+<LI> The temperature on the restart files is virtual equivalent potential temperature,
which requires (at least) surface pressure, specific humidity, and sensible temperature
to calculate.
<LI> CAM does not call the initialization routines when restart files are used,
@@ -141,8 +167,9 @@
<P>
The DART interfaces to CAM and many of the other CESM components
have been integrated with the CESM set-up and run scripts.
+<A NAME="SetupScripts"></A>
Unlike previous versions of DART-CAM, CESM runs using its normal scripts,
-then stops and calls a DART script which runs a single assimilation step,
+then stops and calls a DART script, which runs a single assimilation step,
then returns to the CESM run script to continue the model advances.
See the <a href="../CESM/model_mod.html">CESM interface
documentation</a> for more information on running DART with CESM.
@@ -151,6 +178,21 @@
Each supported CESM version has similar, but unique,
sets of set-up scripts and CESM SourceMods.
Those generally do not affect the cam/model_mod.f90 interface.
+Current (April, 2015) set-up scripts are:
+<UL>
+<LI> CESM1_2_1_setup_pmo: sets up a perfect_model_mod experiment,
+which creates synthetic observations from a free model run,
+based on the user's somewhat restricted choice of model, dates, etc.
+The restrictions are made in order to
+streamline the script, which will shorten the learning curve for new users.
+<LI> CESM1_2_1_setup_pmo_advanced: same as CESM1_2_1_setup_pmo,
+but can handle more advanced set-ups: recent dates (non-default forcing files),
+refined-grid CAM-SE, etc.
+<LI> CESM1_2_1_setup_hybrid: streamlined script (see CESM1_2_1_setup_pmo)
+which sets up an ensemble assimilation using CESM's multi-instance capability.
+<LI> CESM1_2_1_setup_advanced: like CESM1_2_1_setup_pmo_advanced,
+but for setting up an assimilation.
+</UL>
</P>
<P>
@@ -215,13 +257,13 @@
Character strings that contain a '/' must be
enclosed in quotes to prevent them from
prematurely terminating the namelist.
+The values shown here are the default values.
</P>
<div class=namelist>
<pre>
&model_nml
output_state_vector = .false.,
- model_version = '4.0',
model_config_file = 'caminput.nc',
cam_phis = 'cam_phis.nc',
cs_grid_file = 'SEMapping_cs_grid.nc'
@@ -233,10 +275,10 @@
state_names_0d = '',
state_names_1d = '',
state_names_2d = 'PS'
- state_names_3d = 'T', 'US', 'VS', 'Q', 'CLDLIQ','CLDICE',
+ state_names_3d = 'T', 'U', 'V', 'Q', 'CLDLIQ','CLDICE',
which_vert_1d = -2 ,
which_vert_2d = -1 ,
- which_vert_3d = 6*1 ,
+ which_vert_3d = 6*1 ,
pert_names = '',
pert_sd = -888888.0,
pert_base_vals = -888888.0,
@@ -246,8 +288,9 @@
max_obs_lat_degree = 90.0,
Time_step_seconds = 21600,
Time_step_days = 0,
+ print_details = .false.
+ model_version = '4.0',
impact_only_same_kind = '',
- print_details = .false.
/
</pre>
</div>
@@ -256,31 +299,27 @@
<br />
<P>
+The names of the fields to put into the state vector come from the CAM initial
+file field names.
The specification of lists of names and numbers for the various dimensions
enables the very flexible definition of the state vector. It can be
done via the namelist, instead of recompiling DART for each different set.
-One hurdle that still remains is that distinct filter initial condition files are
-necessary for each distinct set of fields which compose the state vector.
+In the CESM+DART framework, filter initial condition files are
+created based on the state vector defined in the namelist.
</P>
<P>
The dimension of these lists is currently hardwired to size 100.
If more fields need to be assimilated (e.g. many chemical species),
look for the integer parameter MAX_STATE_NAMES in the source code
-and change it to a long enough value and recompile DART.
-Longer term we intend to investigate using 2 different namelists
-inside model_mod; one for setting the length of the lists and another
-to actually read in the data which fills the lists.
+and change it to a large enough value and recompile DART.
</P>
<P>
-The values for which_vert_#d is described in
-DART/location/threed_sphere/location_mod.html.
-</P>
+The values for which_vert_#d are described in
-<P>
-The names of the fields to put into the state vector come from the CAM initial
-file field names.
+<em class=file>location_mod</em>:<A HREF="../../location/threed_sphere/location_mod.html#get_close_obs">
+location_mod.html.</a>
</P>
<div>
@@ -300,12 +339,6 @@
prognostic flavor (gridded data) for easier plotting (recommended).
</TD></TR>
-<TR><TD> model_version </TD>
- <TD> character(len=128) </TD>
- <TD>The number of the CAM version being used, i.e. '3.0.7'.
-(no letters allowed, so rename cam3_0_p1)
-</TD></TR>
-
<TR><TD> model_config_file </TD>
<TD> character(len=128) </TD>
<TD>CAM initial file used to provide configuration information, like
@@ -315,19 +348,19 @@
<TR><TD> cam_phis </TD>
<TD> character(len=128) </TD>
- <TD>CAM topography file for the Eularian and Finite Volume versions.
+ <TD>CAM topography file for the Eulerian and Finite Volume versions.
Not used in the Spectral Element version.
</TD></TR>
<TR><TD> cs_grid_file </TD>
<TD> character(len=128) </TD>
- <TD>CAM Spectral Element grid file. Not used in the Eularian or
+ <TD>CAM Spectral Element grid file. Not used in the Eulerian or
Finite Volume versions.
</TD></TR>
<TR><TD> homme_map_file </TD>
<TD> character(len=128) </TD>
- <TD>CAM Spectral Element mapping file. Not used in the Eularian or
+ <TD>CAM Spectral Element mapping file written by CAM-SE. Not used in the Eulerian or
Finite Volume versions.
</TD></TR>
@@ -348,37 +381,32 @@
<TR><TD>which_vert_#d,<br /> #=1,2,3 </TD>
<TD>integer, dimension(100) </TD>
<TD>Vertical location types of fields in state_names_#d.
-See the <a href="../../location/sphere_threed/location_mod.html">3D sphere location</a>
+See the <a href="../../location/threed_sphere/location_mod.html">3D sphere location</a>
documentation for the mapping of integer values to vertical location types.
</TD></TR>
<TR><TD> pert_names </TD>
<TD>character(len=8), dimension(100) </TD>
<TD>To make filter generate an ensemble from a single model state by
-randomly perturbing it, list the field(s) to be perturbed here. To make
-trans_pv_sv_pert0 reset a single field to a constant value, list that field
-here. Trans_pv_sv_pert0 would be run as many times as there are ensemble
-members (see pert_base_vals), in order to provide spread to that variable.
+randomly perturbing it, list the field(s) to be perturbed here
+(and see the DART namelist settings for a
+<A HREF="#PerturbedEnsemble">perturbed initial ensemble)</A>.
</TD></TR>
<TR><TD> pert_sd </TD>
<TD>real(r8), dimension(100) </TD>
<TD>If positive, it's the standard deviation of the perturbation for each
-field in the pert_names list (filter). If negative, then pert_names can
-contain only one entry, and that field will be set to a different constant
-value for each ensemble member (trans_pv_sv_pert0). Those values come from
-pert_base_vals. Defaults to a MISSING real value and unused unless
-pert_names is set.
+field in the pert_names list (filter). Unused unless pert_names is set.
</TD></TR>
<TR><TD> pert_base_vals </TD>
<TD>real(r8), dimension(100) </TD>
- <TD>If pert_sd is negative, this is the list of values to use for
+ <TD>If pert_sd is positive, this the list of values to which the field(s)
+listed in pert_names will be reset if filter is told to create an ensemble
+from a single state vector. Otherwise, it's is the list of values to use for
each ensemble member when perturbing the single field named in pert_names.
-Otherwise, it's the list of values to which the field(s) listed in pert_names
-will be reset if filter is told to create an ensemble from a single state
-vector. Defaults to a MISSING real value and unused unless pert_names is set
-and pert_base_vals is /= this missing value (-888888.0d0).
+Unused unless pert_names is set
+and pert_base_vals is not the DART missing value.
</TD></TR>
<TR><TD> max_obs_lat_degree </TD>
@@ -389,7 +417,7 @@
<TR><TD> vert_coord </TD>
<TD> character(len=8) </TD>
<TD>The vertical coordinate to which all vertical locations are converted
-in model_mod. "log_invP" ("scale height" for WACCM) or "pressure" (default).
+in model_mod. "log_invP" ("scale height" especially for WACCM) or "pressure".
</TD></TR>
<TR><TD> highest_obs_pressure_Pa </TD>
@@ -402,15 +430,18 @@
<TR><TD> highest_state_pressure_Pa </TD>
<TD> real(r8) </TD>
- <TD>Influence of all obs on model points higher than this is reduced.
-NOTE that this
-has a non-backwards incompatible change from previous versions. It is now
-specified in Pascals, not millibars. Divide by 100 to convert the units.
+ <TD>Influence of all observations on model points higher than this is reduced.
+NOTE that this has a non-backwards incompatible change from previous versions.
+It is now specified in Pascals, not millibars. Divide by 100 to convert the units to mb.
+Details of calculating the minimum value recommended for a given vertical grid
+and set of DART namelist variables can be found in
+<A HREF="./doc/highest_state_p_Pa.pdf">highest_state_p_Pa.pdf</A>,
</TD></TR>
<TR><TD> Time_step_seconds </TD>
<TD> real(r8) </TD>
- <TD>Minimum forecast duration (the part < 1 day)
+ <TD>Minimum forecast duration (the part < 1 day).
+CESM1_2_1_setup_{hybrid,advanced} assume that this is 6 hours (21600 s).
</TD></TR>
<TR><TD> Time_step_days </TD>
@@ -418,19 +449,24 @@
<TD>Minimum forecast duration (the part > 24*3600 sec)
</TD></TR>
+<TR><TD> print_details </TD>
+ <TD> logical </TD>
+ <TD>If true, print out detailed information about the sizes, shapes, offsets,
+etc. of items in the CAM state vector. If false, just print out the names of
+the fields as they are read into the state vector.
+</TD></TR>
+
+<TR><TD> model_version </TD>
+ <TD> character(len=128) </TD>
+ <TD> Not used in CESM+DART; deprecated. The number of the CAM version being used, i.e. '3.0.7'.
+</TD></TR>
+
<TR><TD> impact_only_same_kind </TD>
<TD> character(len=32) </TD>
<TD>Name of one observation kind which can only affect
state variables of the same kind.
</TD></TR>
-<TR><TD> print_details </TD>
- <TD> logical </TD>
- <TD>If true, print out detailed information about the sizes, shapes, offsets,
-etc of items in the CAM state vector. If false, just print out the names of
-the fields as they are read into the state vector.
-</TD></TR>
-
</TBODY>
</TABLE>
</div>
@@ -440,21 +476,244 @@
<!--==================================================================-->
+<A NAME="SetupVariations"></A>
+<H2> Setup Variations </H2>
+<P>
+The variants of CAM require slight changes to the setup scripts
+(in $DART/models/cam/shell_scripts) and in the namelists
+(in $DART/models/cam/work/input.nml).
+From the DART side, assimilations can be started from a pre-existing ensemble,
+or an ensemble can be created from a single initial file before the first assimilation.
+In addition, there are setup differences between 'perfect model' runs,
+which are used to generate synthetic observations, and assimilation runs.
+Those differences are extensive enough that they've been coded into separate
+<A HREF="#SetupScripts">setup scripts</A>:
+</P>
+
+<P>
+Since the CESM compset and resolution, and the initial ensemble source
+are essentially independent of each other, changes for each of those
+may need to be combined to perform the desired setup.
+</P>
+
+<A NAME="#PerturbedEnsemble"></A>:
+<P>
+The default values in work/input.nml and shell_scripts/CESM1_2_1_setup_{pmo,hybrid}
+are set up for a CAM-FV, single assimilation cycle
+using the default values as found in model_mod.f90 and
+starting from a single model state, which must be perturbed into an ensemble.
+The following are suggestions for setting it up for other assimilations.
+Namelist variables listed here might be in any namelist within input.nml.
+</P>
+
+
+<H3> CAM-FV </H3>
+<P>
+<a name="CAM-FV"></a>
+If built with the FV dy-core, the number of model top levels with extra diffusion
+in CAM is controlled by div24del2flag. The recommended minium values of
+highest_state_pressure_Pa come from that variable, and cutoff*vert_normalization_X:
+<pre>
+ 2 ("div2") -> 2 levels -> highest_state_pressure_Pa = 9400. Pa
+ 4,24 ("del2") -> 3 levels -> highest_state_pressure_Pa = 10500. Pa
+</pre>
+<pre>
+ vert_coord = 'pressure'
+ state_num_1d = 0,
+ state_num_2d = 1,
+ state_num_3d = 6,
+ state_names_1d = ''
+ state_names_2d = 'PS'
+ state_names_3d = 'T', 'US', 'VS', 'Q', 'CLDLIQ', 'CLDICE'
+ which_vert_1d = 0,
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list