% Matlab file that runs the GISS 1-dimensional model. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Clears all possible leftovers from previous use % of Matlab so as to avoid potential problems. clear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Various parameters : % You can choose to change them as you like, % except when noted otherwise. % However, for your convienience, remember to % add the semi-columns in the end of each line. % Depth of the mixed layer (in meters). Hml = 100 ; % Depth of the thermocline (in meters). H = 900 ; % Number of layers in the thermocline. Nz = 30 ; 'Pressing the return key will keep the default values' % Heat diffusivity coefficient. % In the first line units are cm^2/sec. % You are allowed to change only the first line, not the second. % In the second line I multiply it by % 10^-4 to convert the units to m^2/sec. k = input('What is the value of the heat diffusivity coefficient (in cm^2/sec) ?'); if isempty(k) k = 1.5 ; end k = k*1.0e-4 ; % Equilibrium temperature change for double CO2 (in degrees Kelvin). Te = input('What is the value of the equilibrium temperature change for doubled CO2 (in C) ?'); if isempty(Te) Te = 3.0 ; end % Data collection interval (in years). % Should be chosen according to the smoothness % of the changes of our forcing function. % No need to make it very small though. % I have initially set it equal to 1 year. T_interval = 1 ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % From now on you are not supposed to change anything, % unless you are looking for trouble! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Give filenames for the input and the output files. in_data = input('What is the name of the input file ? : ','s') ; if isempty(in_data) ' You did not give the name of the input file ' ' This forces me to suspend the program. ' ' You will have to start all over again. Sorry! ' break end %out_data = input('What is the name of the output file ? : ','s') ; %if isempty(out_data) %' You did not give the name of the output file ' %' This forces me to suspend the program. ' %' You will have to start all over again. Sorry! ' % break %end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Read the various variables that enter in the forcing function. fid=fopen(in_data,'rt') ; fgetl(fid) ; inputs = fscanf(fid,'%f %f %f %f %f %f %f %f',[8 inf]) ; inputs = inputs' ; clear fid %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Other parameters: % Those you ABSOLUTELY cannot change. % Depth of each layer in the thermocline (in meters). dz = H/Nz ; % Time step for forward integration (in sec). % I set an explicit upper limit for dt of 0.1 years. if 1/2*dz^2/(2*k)/(24*3600*365) > 0.1 dt=0.1*(24*3600*365) else dt = 24*3600*365*1/round((24*3600*365)/(1/2*dz^2/(2*k))) ; end % Time interval of integration (end date, begin date). % Taken from the first and last date that appear in the input file. Date_begin = inputs(1,1) ; Date_end = inputs(length(inputs(:,1)),1) ; % Weird date variable (in sec). Don't worry about it. time = [Date_begin*24*3600*365 : dt : Date_end*24*3600*365] ; % Number of time steps. Ntime = (Date_end-Date_begin)*24*3600*365/dt ; % Number of time steps between data collection. N_interval = T_interval*24*3600*365/dt ; % Heat capacity of water. % In the first line units are cal/(cm^3 degrees Kelvin). % In the second line I multiply by the appropriate % factor to convert the units to Joules/(m^3 degrees Kelvin) C = 1 ; C = C*4.2e6 ; % Heat conductivity coefficient (in Watts/(m degrees Kelvin). % If I hadn't changed the units of "C" and "k", then the units % would be cal/(m sec degrees Kelvin), and since "C" would be % equal to unity, "lamda" would be numerically equal to "k" lamda = C*k ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Interpolate the input variables onto the % time grid that has just been created. DCO2 = interp1(inputs(:,1).*24*3600*365,inputs(:,2),time,'linear') ; CO2 = DCO2 + 293 ; DCH4 = interp1(inputs(:,1).*24*3600*365,inputs(:,3),time,'linear') ; DN2O = interp1(inputs(:,1).*24*3600*365,inputs(:,4),time,'linear') ; DCC2 = interp1(inputs(:,1).*24*3600*365,inputs(:,5),time,'linear') ; DCC3 = interp1(inputs(:,1).*24*3600*365,inputs(:,6),time,'linear') ; DV = interp1(inputs(:,1).*24*3600*365,inputs(:,7),time,'linear') ; DS = interp1(inputs(:,1).*24*3600*365,inputs(:,8),time,'linear') ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize the run and take first measurements. % You will find the definitions of the various variables % in later parts of the program. F(1) = 0 ; Fd(1) = 0 ; DTml(1) = 0 ; j=[1:Nz+1] ; DTtherm(j) = 0.*ones(size(j)) ; count = 1 ; data_collect(count,1) = time(1)./(24*3600*365) ; data_collect(count,2) = DTml(1) ; data_collect(count,3) = F(1) ; data_collect(count,4) = Fd(1) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Beginning of time step loop. for i=2:Ntime+1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Empirical formula giving heat flux from % the atmosphere into the earth's surface. % In the first formula the units are cal/(min cm^2). % I then multiply by the appropriate factor to % convert the units to Watts/m^2. F(i) = 2.6e-5.*DCO2(i-1)./(1+0.0022*CO2(i-1)).^0.6 ... - 5.88e-3/Te.*DTml(i-1) ... + 0*3.685e-4/Te^2.*DTml(i-1).^2 ... - 0*4.172e-7/Te.*DCO2(i-1).*DTml(i-1) ... + 1.197e-3.*DCH4(i-1).^0.5 ... + 5.88e-3.*DN2O(i-1).^0.6 ... + 3.15e-4.*DCC3(i-1) ... + 3.78e-4.*DCC2(i-1) ... - 1.197e-4.*DCH4(i-1).*DN2O(i-1) ... - 2.40e-2.*DV(i-1) ... + 2.10e-3.*DV(i-1).^2 ... - 0*1.17e-3/Te.*DTml(i-1).*DV(i-1) ... + 3.184e-1.*DS(i-1) ; F(i) = F(i).*4.2*1.0e4*1/60 ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Heat flux from the thermocline into the mixed layer (in Watts/m^2). Fd(i) = lamda.*(DTtherm(2)-DTml(i-1))./dz ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Forward integration to get the next time step temperatures % at the mixed layer and the various layers of the thermocline. DTml(i) = DTml(i-1) + (F(i) + Fd(i)).*dt/(C*Hml) ; DTtherm(1) = DTml(i) ; j=[2:Nz] ; DTtherm(j)=DTtherm(j) + ... k.*(DTtherm(j+1)-2.0.*DTtherm(j)+DTtherm(j-1))./dz^2*dt ; DTtherm(Nz+1) = DTtherm(Nz) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Data collection. if round(i/N_interval) == i/N_interval count=count+1 data_collect(count,1) = time(i)./(24*3600*365) ; data_collect(count,2) = DTml(i) ; data_collect(count,3) = F(i) ; data_collect(count,4) = Fd(i) ; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of time step loop. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Save data in appropriate file, called "data.dat". %save out_data data_collect -ascii % Plot change in temperature subplot(3,1,1) plot(data_collect(:,1), data_collect(:,2)) title('\Delta T - Mixed Layer') ylabel('^\circ C') subplot(3,1,2) plot(data_collect(:,1), data_collect(:,3)) title('Heat Flux from Atmosphere \rightarrow Ocean') ylabel('W/m^2') subplot(3,1,3) plot(data_collect(:,1), data_collect(:,4)) title('Heat Flux from Thermocline \rightarrow Mixed Layer') ylabel('W/m^2') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program.