%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CMIP6 Grazing Rates %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Created November 2022
% Read Me: Computed Diagnostic Grazing Rates and Auxilliary Variable for CMIP6 models
% USed to create all data products in Comm Earth And Env paper

%% Startup
addpath(genpath('/Users/roh029/Desktop/Matlab/SOLACE'))
addpath(genpath('/Users/roh029/Desktop/Matlab/Utilities'))

%% Plot
%  Pick Model
cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv

%  plotname  = 'ACCESS-ESM-1-5';   % WOMBAT
%  plotname  = 'CanESM5-CMOC';     % CMOC
%  plotname  = 'NorESM2-ESM-LM';   % HAMOCC
%  plotname  = 'MIROC-ES2L';       % OECO
%  plotname  = 'CESM2';            % MARBL   
%  plotname  = 'CanESM5-CanOE';    % CanOE
%  plotname  = 'UKESM1-0-LL';      % Medusa2.1
%  plotname  = 'IPSL-CM6A-LR';     % PISCESv2
%  plotname  = 'CNRM-ESM2-1';      % PISCESv2
%  plotname  = 'CMCC-ESM2';        % BFM5.2
%  plotname   = 'GFDL-ESM4';       % Cobaltv2

load(plotname)
load coastlines

% Pick Vars:
% 1. Phytoplankton Biomass 2. Zooplankton Biomass  
% 3. NPP                   4. GSP                  5. EP
% 6. Spec. Division Rate   7. Spec. Grazing Rate
% 8. Z:P                   9. EP/NPP               10 GSP/NPP
CoreVar_Names(8:10) = {'Z:P','EP/NPP','GSP/NPP'}; 
CoreVar_Inv_Units(8:10) = {'(frac)','(frac)','(frac)'};
CoreVar_Int_Units(8:10) = {'(frac)','(frac)','(frac)'};

for makeplot = 1

seasnames = {'DJF','MAM','JJA','SON'};
seas(1,:)      = [12,1,2];     seas(2,:)      = [3:5];   
seas(3,:)      = [6:8];     seas(4,:)      = [9:11];

% Resize for mapping
[lat_size,lon_size] = size(Lon);
Lon_rz = squeeze([Lon,Lon,Lon]);
Lat_rz = squeeze([Lat,Lat,Lat]);
Lon_rz(:,1:lon_size) = Lon_rz(:,1:lon_size)-360;
Lon_rz(:,lon_size*2+1:end) = Lon_rz(:,lon_size*2+1:end)+360;

for p=3%1:3
  if p==1; vars = [3,4,5]; elseif p==2; vars = [1,2,8]; elseif p==3; vars = [6,7,9,10]; end


% Create Figure
  figure('Units','normalized','Position',[0 0 2 2]); clc
  set(gcf,'color','w');
  
  for v = 1:length(vars)
      if vars(v)==8 
      plotvar = (squeeze(CoreVar_Inv(2,:,:,:)))./(squeeze(CoreVar_Inv(1,:,:,:))); cRange = [0,2];%prctile(plotvar(:),[2,95]);
      elseif vars(v)==9
      plotvar = (squeeze(CoreVar_Inv(5,:,:,:)))./(squeeze(CoreVar_Inv(3,:,:,:)));  cRange = [0,.7]; %cRange = prctile(plotvar(:),[2,96]);
      elseif vars(v)==10 
      plotvar = (squeeze(CoreVar_Inv(4,:,:,:)))./(squeeze(CoreVar_Inv(3,:,:,:))); cRange = [0,1.5]; %cRange = prctile(plotvar(:),[2,98]);
      else
      plotvar = (squeeze(CoreVar_Inv(vars(v),:,:,:))); cRange = prctile(plotvar(:),[2,98]);
      end
      
      
      for s =1:4
          Plotvar = squeeze(nanmean(plotvar(seas(s,:),:,:)));
          Plotvar_rz = cat(2,Plotvar,Plotvar,Plotvar);
          
          h= subplot(4,4,(v-1)*4+s);
                % Adjust Positions [left bottom width height]
                posnew = get(h, 'Position'); 
                %posnew(4) = posnew(4) +.05;
                if     s==1; posnew(1) = posnew(1) - 0.05; 
                elseif s==2; posnew(1) = posnew(1) - 0.11; 
                elseif s==3; posnew(1) = posnew(1) - 0.17;
                elseif s==4; posnew(1) = posnew(1) - 0.222; 
                       PosPlot(v,:)=posnew;
                end
                set(h, 'Position', posnew)
                
          axesm('miller','maplatlim',[-75 75],'maplonlim',[-260, 100]); 
          pcolorm(Lat_rz,Lon_rz,Plotvar_rz); shading flat; hold on
          plotm(coastlat,coastlon,'k'); patchm(coastlat,coastlon,[114,114,114]/256);
          
          tightmap 
          %m_proj('miller','lon',[-260,100]);
          %m_pcolor(Lon_rz,Lat_rz,Plotvar_rz); shading flat; hold on
          %m_coast('patch',[.6 .6 .6]);
          %m_grid('xaxis','middle','yticklabels',[],'xticklabels',[]);
          
          % Colormap
          if      any(vars(v)==[1,3,6]) ; colors = cmocean('algae'); 
                        colors =[multigradient([ 1,1,1; colors(1,:)],'length',40); colors]; caxis(cRange)
          elseif  any(vars(v)==[2,4,7,10]); colors =  cmocean('matter'); 
                         colors =[multigradient([ 1,1,1; colors(1,:)],'length',40); colors]; caxis(cRange)
          elseif  any(vars(v)==[5,9]); colors=cmocean('dense'); 
                         colors =[multigradient([ 1,1,1; colors(1,:)],'length',40); colors]; caxis(cRange)
          elseif  any(vars(v)==[8]); colors =  cmocean('matter'); 
                         colors =[multigradient([ 1,1,1; colors(1,:)],'length',40); colors];  caxis(cRange)
                    
          end
          colormap(gca,colors)
         
          %set(gca,'ColorScale','log')
          if s==4; cb= colorbar; 
               set(get(cb,'title'),'string',CoreVar_Inv_Units{vars(v)},'interpreter','latex','fontsize',25);
               set(cb,'TickLabelInterpreter','latex','fontsize',25)
          end
          
          
          % Title; 
          if s==1; ylabel({CoreVar_Names{vars(v)};CoreVar_Inv_Units{vars(v)}}  ,'fontsize',25,'interpreter','latex'); end
          if v ==1; title(seasnames{s},'fontsize',25,'interpreter','latex'); end
      end  
  end   
  
  % Add global integrals
   for v = 1:length(vars)
       
      if vars(v)==8 
      plotvar = (squeeze(CoreVar_Int(:,2,:)))./(squeeze(CoreVar_Int(:,1,:)));
      elseif vars(v)==9
      plotvar = (squeeze(CoreVar_Int(:,5,:)))./(squeeze(CoreVar_Int(:,3,:)));
      elseif vars(v)==10 
      plotvar = (squeeze(CoreVar_Int(:,4,:)))./(squeeze(CoreVar_Int(:,3,:)));
      else
      plotvar = (squeeze(CoreVar_Int(:,vars(v),:)));
      end
       
    ax(v) =  axes;         % lef              bot           wid  h
    set(ax(v), 'Position',[ PosPlot(v,1)+PosPlot(v,3)+.045  ...
                            PosPlot(v,2), PosPlot(v,3)-.02, PosPlot(v,4) ]); 
    plot(1:12,squeeze(plotvar(1,:)),'ko-','linewidth',2);  hold on
      yline(mean(squeeze(plotvar(1,:))),'k')
    plot(1:12,squeeze(plotvar(2,:)),'bo-','linewidth',2);  hold on
    plot(1:12,squeeze(plotvar(3,:)),'go-','linewidth',2);  hold on
    ylabel({CoreVar_Names{vars(v)};CoreVar_Int_Units{vars(v)}}  ,'fontsize',18,'interpreter','latex'); 
    
    if any(vars(v)==[1:5]); title('Spatial Integral','fontsize',18,'interpreter','latex'); end
    if any(vars(v)==[6:10]); title('Spatial Mean','fontsize',18,'interpreter','latex'); end
    %if v==4; legend('Global','Northern Hemisphere','Southern Hemisphere', 'location','southoutside'); end
    xticks([1:12]); xlim([1,12])
    xticklabels({'J','F','M','A','M','J','J','A','S','O','N','D'})
    set(ax(v),'fontsize',18,'TickLabelInterpreter','latex')
   end
   
% Title
ax(20) =  axes;         % lef              bot           wid  h
set(ax(20), 'Position',[ 0,0,.06,1]); 
text(.5,.63 ,plotname, 'horizontalalignment','center','interpreter','latex','fontsize',80,'Rotation',90)
xticks([]); yticks([]); 
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create Data Products %%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%for RunAllDiagnostic=1
latmin =-60; latmax = 60;
 dth =200; % Depth threshold to integrated down the water columm ... Perhaps cause thats as deep as CESM2 is saved?


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for proc =1
%% NorESM2-ESM-LM: HAMOCCv2 (NOR) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_NorESM2-LM_historical_r2i1p1f1_gn
% Variables Used: phyc | zooc | intpp | epc100 |
%                 areacello | volcello 
% Caveats: Pthresh from Personal correspondence
% Notes:   Biomass is saved on regrided grid (gr) and fluxes/rates on
%          native grid (gn). Cell area is same on both but native grid has density
%          levels. Constant cell volume from gr is used here.
fname = 'NorESM2-ESM-LM';

for Compute = 1
fprintf('\n NorESM2-ESM-LM: HAMOCCv2L (NOR) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    % Volume
    volcello = ncread('volcello_Ofx_historical_NorESM2-LM_r1i1p1f1_gr.nc','volcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Volume = permute(volcello(:,:,1:70),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)

    % Area
    areacello = ncread('areacello_Ofx_historical_NorESM2-LM_r1i1p1f1_gn.nc','areacello');
    % float areacello(i (lon), j (lat)) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    Cell_Thick = Cell_Volume./permute(repmat(Cell_Area,[1,1,70]),[3,1,2]);
 
    % Depth Level
    depthlev = ncread('volcello_Ofx_historical_NorESM2-LM_r1i1p1f1_gr.nc','lev'); 
    % bin 13 is 100m
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = lon.var.longitude; Lat= lat.var.latitude;
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zooc    = zooc*1000;            % mol/m3   -> mmol/m3
    phyc    = phyc*1000;            % mol/m3   -> mmol/m3
    intpp   = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100  = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
     %zooc(zooc<1e-50)=nan;       
    zooc(zooc>10e15) = nan; 
    %phyc(phyc<1e-50)=nan;       
    phyc(phyc>10e15) = nan; 
    %intpp(intpp<1e-5)=nan;     
    intpp(intpp>10e15) = nan; 
   % epc100(epc100<1e-5)=nan;  
    epc100(epc100>10e15) = nan; 
    
    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    zooc    = circshift(zooc,-(lonoff-1),4);
    phyc = circshift(phyc,-(lonoff-1),4);
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'zooc', 'phyc', 'intpp', 'epc100', '-append')
    save(fname, 'Cell_Thick','Cell_Area','Cell_Volume','-append');
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
    % PRIMARY Source:  Tjiputra (2020); tjiputraOceanBiogeochemistryNorwegian2020       
                    % 1P1Z with single-prey grazing
                      % Phosphate Units
                      % P:N:C = 1:16:122 
    % Parameters in:  Tjiputra (2020); tjiputraOceanBiogeochemistryNorwegian2020    
                      % Unreported paremeter are the same as version 1
    % Equations in:   ilyinaGlobalOceanBiogeochemistry2013
                       % Maier-Reimer et al 2005; maier-reimerHAMburgOceanCarbon2005      (NORESM1)
                       % Via Schwinger 2016; schwingerEvaluationNorESMOCVersions2016 (HAMOCC5.1)
                       % "To avoid negative phytoplankton concentrations due to the constant forward time stepping,
                       %      zooplankton grazing is calculated semi-implicitely:"
   
    % Mortality:    Zloss = .06Z +.025Z2
                         % Losses to excretion, mortality  –  Tjiputra (2020)
                         % quadratic mortaltiy converted from 3e-6 1/((kmolPm3)d)
                         % NOTE: units of morality term suggest quadratic
                         % in shcwing and tjipura, but linear in MR2005--
                         % Assue it was changed

                       
    % Temp limiation and thresholds
    pthresh = 1e-11 *1e6 .*122;%             % Converting from 1e-11 kmol Pm-3 -> From Personal comunication with Jorg Scheinger
                                             % @ Surface Units of kmol P m−3. Value not given in Maier-Reimer et al 2005
                        
    % Zooplankton Group 1 (Generic Zooplankton)
    % MM-II w/ single prey repsonse 
    % Specific Grazing Rates on:
    % Phytoplankton --  K: 9.78 (mmol C/m3); g: 1.2; EQ: g = g*(P'/P'+ K)
    %    Conversion: K = 8 10^-8 kmol P m^3  ->> 8*10^-8 .*1000000 .* 122 mmol C m^3
        pstar     =  max(0,phyc-pthresh);
        g_Z_P = 1.2.* (pstar./(phyc + 9.78));
           
    % Total Grazing 
        GSP = g_Z_P.*zooc;  

    % Total Trophic Transfer 
        w =0.6; % Assimilation Efficiency
        e = .85; % ingestion(vs eggestion fraction)
        GSP_ht = GSP.*e.*w; 
        
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./zooc;
     
    % Bulk mortality 
      zmin = 0.0012;   % Minimum Z concentration. 1e-11 kmolP/m3 - > mmolC/m3
      zstar = max(0,zooc-zmin);
      Bz =0.06;        % Excretion Rate % ilyinaGlobalOceanBiogeochemistry2013
      mZ = .025;       % Mortality Rate 3e6 1/d/(kmolP/m3)  -> 1/d/(mmolC/m3) %ilyinaGlobalOceanBiogeochemistry2013
      z_loss       =  .06*zstar +.025.*zstar.^2;
      z_loss_spec  =  z_loss./zooc;
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec','GSP_ht', 'GSP','z_loss','z_loss_spec', '-append')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phyC_inv    = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv    = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
   
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv;  
    
    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
     
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'G',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP','Z Loss rate','Specific Z loss rate'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)','($d^{-1}$)'};
    
  % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phyC_inv',      'zooC_inv','gsp_inv', 'gsp_ht_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CANESM5 - CMOC (CAN) %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_CanESM5_historical_r1i1p2f1_gn
% DOI: https://www.wdc-climate.de/ui/cmip6?input=CMIP6.CMIP.CCCma.CanESM5.historical
%       doi:10.22033/ESGF/CMIP6.3610
% Variables Used: phyc | zooc | intpp | epc100 |
%                 areacello | thkcello 
% Caveats: 
% Notes:   

fname = 'CanESM5-CMOC';

for Compute = 1
fprintf('\n CanESM5-CMOC: CMOC (CAN) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    
    % Area
    areacello = ncread('areacello_Ofx_CanESM5_historical_r2i1p2f1_gn.nc','areacello');
    % float areacello(i (lon), j (lat)) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_CanESM5_historical_r1i1p1f1_gn.nc','thkcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Thick = permute(thkcello(:,:,1:45),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
  
    % Volume
    Cell_Volume = Cell_Thick.*permute(repmat(Cell_Area,[1,1,45]),[3,1,2]); % Unit: m^3, Dim: (level, lat, lon)
    
    % Depth Level
    depthlev = ncread('thkcello_Ofx_CanESM5_historical_r1i1p1f1_gn.nc','lev'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = lon.var.longitude; Lat= lat.var.latitude;
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    
    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zooc    = zooc*1000;            % mol/m3   -> mmol/m3
    phyc    = phyc*1000;            % mol/m3   -> mmol/m3
    intpp   = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100  = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %zooc(zooc<1e-50)=nan;       
    zooc(zooc>10e15) = nan; 
    %phyc(phyc<1e-50)=nan;       
    phyc(phyc>10e15) = nan; 
    %intpp(intpp<1e-5)=nan;     
    intpp(intpp>10e15) = nan; 
   % epc100(epc100<1e-5)=nan;  
    epc100(epc100>10e15) = nan; 
    
    %
    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    zooc    = circshift(zooc,-(lonoff-1),4);
    phyc    = circshift(phyc,-(lonoff-1),4);
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'zooc', 'phyc','intpp', 'epc100', '-append')
    save(fname, 'Cell_Thick','Cell_Area','Cell_Volume','-append');
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
    % PRIMARY Source: zahariev et al. (2008);  zaharievPreindustrialHistoricalFertilization2008
        % 1P1Z with single-prey response 
           % Phytoplankton
           % Zooplankton
           % Detrius       
        % N units
          % 106:16 C:N Ratio (Redfield) 
                                             
     % Parameters in:  zaharievPreindustrialHistoricalFertilization2008
     %                 Table 1
                   
     % Equations in:   zaharievPreindustrialHistoricalFertilization2008
     %                 Equation 10
 
 
     % Z1: "Zooplankton"  Grazing On
          % Prey: Phytoplankton                 
            % "Phytoplankton"
            %  Single-prey w/ switching  
            %  g   =  2.0 (1/d); 
            %  K   =  1.325 (mmolC/m3); Converted from 0.2 (mmolN/m3)
             g_Z_P = 2.0.*phyc.^2./(1.325.^2 + phyc.^2);

    % Total Grazing 
        GSP = g_Z_P.*zooc;  
        
    % Total Transfer to higher trophic levels
        g =0.7; % Assimilation efficeincy.
        GSP_ht = g.*GSP; 
        
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./zooc;
        
      % Bulk mortality 
      Bz = .2 +.05; % 1/d; Metabolic losses to nitrogen and detritus, repsectively
      Mz = 0.015; %  0.1   1/d/(mmolN/m3)  --> 1/d/(mmolC/m3) 
      
      z_loss       =  Bz.*zooc +Mz.*zooc.^2;
      z_loss_spec  =  z_loss./zooc;
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec','GSP_ht', 'GSP','z_loss','z_loss_spec', '-append')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

     % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phyC_inv    = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv    = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
   % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
   
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
     % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
    
     % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phyC_inv',      'zooC_inv','gsp_inv', 'gsp_ht_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ACCESS-WOMBAT (AUS) %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_ACCESS-ESM1-5_historical_r1i1p1f1_gn
% Variables Used: phyc | zooc | intpp | epc100 |
%                 areacello | volcello | thkcello
% Caveats: phyc | zooc | intpp were all pulled from gadi and processed
%          seperately as bottom on CMIP6_Graze_Clim script.
% Notes:   epc100 was hosted on the public server and already in units of molC/m2/s; however, phyc, zooc, 
%          and intpp were derived from raw output (phy, zoo, pprod_gross_2d) hosted on gadi in mmolN units.
%          These were all convereted to mmol Carbon in CMIP6_Graze_Clim on Cheyyene. Then saved as such in 
%          ACCESS-ESM-1-5_Clims          
          
fname = 'ACCESS-ESM-1-5';

for Compute = 1
fprintf('\n ACCESS-ESM1-5: WOMBAT (AUS) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    % Volume
    volcello = ncread('volcello_Ofx_ACCESS-ESM1-5_historical_r1i1p1f1_gn.nc','volcello');
    % float volcello(lev, i (lon), j (lat)) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Volume = permute(volcello(:,:,1:50),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)

    % Area
    areacello = ncread('areacello_Ofx_ACCESS-ESM1-5_historical_r1i1p1f1_gn.nc','areacello');
    % float areacello(i (lon), j (lat)) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_ACCESS-ESM1-5_historical_r1i1p1f1_gn.nc','thkcello');
    Cell_Thick = permute(thkcello(:,:,1:50),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
    % Cell_Thick_Cehck = Cell_Volume./permute(repmat(Cell_Area,[1,1,50]),[3,1,2]); 
 
    % Depth Level
    depthlev = ncread('volcello_Ofx_ACCESS-ESM1-5_historical_r1i1p1f1_gn.nc','lev'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = lon.var.longitude; Lat= lat.var.latitude;
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zooc    = zooc;                 % mmolC/m3    -> mmolC/m3       
    phyc    = phyc;                 % mmolC/m3    -> mmolC/m3       
    intpp   = intpp*60*60*24;       % mmolC/m2/s -> mmolC/m3/d     
    epc100  = epc100*1000*60*60*24; % molC/m2/s  -> mmolC/m3/d
    thetao  = thetao./6.625;        % Accidntently applied N->C Conversion while processing climatolgoy. Correct for that. 
    
    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    % remove anything with more export than pp
    %zooc(zooc<1e-10)=0;      
    zooc(zooc>10e15) = nan; 
    %phyc(phyc<1e-10)=0;      
    phyc(phyc>10e15) = nan; 
    %intpp(intpp<1e-10)=0;    
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-10)=0;   
    epc100(epc100>10e15) = nan; 
    thetao(thetao>10e15) = nan;  
    

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    zooc    = circshift(zooc,-(lonoff-1),4);
    phyc    = circshift(phyc,-(lonoff-1),4);
    thetao   = circshift(thetao,-(lonoff-1),4);
    
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
 
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'zooc', 'phyc', 'intpp', 'epc100','thetao', '-append')
    save(fname, 'Cell_Thick','Cell_Area','Cell_Volume','-append');
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
    % PRIMARY Source: Law et al (2017);  lawCarbonCycleAustralian2017
                    % 1P1Z with single-prey response
                    % P units
                      % But parameters reported (And presumbaly input) in N units. 
                      % Explicitly uses 106:16:1 P:N:C       
                    
    % Parameters in:  Law et al (2017);  lawCarbonCycleAustralian2017
                     % Operates in Phosphate but epsilon is given in Nitrate units
                     % "The equations of WOMBAT are given in Oke et al. (2013, Appendix B)" - pg 2570
    % Equations in:   Oke et al (2013);  okeEvaluationNearglobalEddyresolving2013
                     % "We assume a Redfield ratio of 1 : 16 : 106 : −172 (P : N : C : O2) with phosphate, nitrate, carbon and oxygen."
                     % Phosphate is macro-nutrient BUT eps is reported (and parameterized) in N.
                        %  Not sure where in source code it is converted thouhgh

                        
     % Z1: "Zooplankton"  Grazing On:
            % "Phytoplankton
            % D-III w/ Single prey repsonse 
            %  g= 1.575 (1/d); 
            %  eps = 1.6 (1/ (mmolN m^6)/d )  % 1/ (mmolN m^6)/d   Reported Units: 1/ (mmolN m^2)/d (here and in Oke et al) 
            %  K = 6.57  (sqrt(gmax_WOMBAT/ eps_WOMBAT).*(106/16)) % (mmolC m^3); converted from Type III Disk, and N to C (w/ 106/16)
               g_Z_P = (1.575.*phyc.^2)./(phyc.^2 + 6.57.^2);
     
    % Total Grazing 
        GSP = g_Z_P.*zooc;  
        
    % Total Transfer up the trophic ladder
      g = .925; % Assimilation efficeicncy
      GSP_ht = g.*GSP;

    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./zooc;
     
    % Bulk mortality 
      Bz = 0.01.*1.066.^thetao; % Excretion rate
      Mz = 0.05;                % Quadratic mortality Rate
      z_loss       =  Bz.*zooc + Mz.*zooc.^2;
      z_loss_spec  =  z_loss./zooc;
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP','z_loss','z_loss_spec', '-append')

    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP','GSP_ht', '-append')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

     % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phyC_inv    = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv    = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
   
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
    
   
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phyC_inv',      'zooC_inv','gsp_inv', 'gsp_ht_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MIROC-ES2L - OECO-v2 (JAP) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_MIROC-ES2L_historical_r1i1p1f2_gn
% Variables Used: phyc | zooc | intpp | epc100 | tos
%                 areacello | volcello | thkcello
% Caveats: OECO-v2 includes a Diazatrophs class but only a single phytoplankton variable was uploaded to the server (phyc), presumbaly because there was 
%          very little diazatroph biomass and they were combined. To compute specfiic grazing rates, we assume diazatroph biomass accounted for 2% of 
%          total phyc based on the mean  partitionioning in CESM2, the model with the most similiar BGC structure and diazatrophs.
% Notes:
fname = 'MIROC-ES2L';

for Compute = 1
fprintf('\n Starting CMIROC-ES2L - OECO-v2 (JAP) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids
    
    % Volume
    volcello = ncread('volcello_Ofx_MIROC-ES2H_historical_r3i1p4f2_gn.nc','volcello');
    % loat volcello(x, y, lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Volume = permute(volcello(:,:,1:63),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)

    % Area
    areacello = ncread('areacello_Ofx_MIROC-ES2H_historical_r2i1p4f2_gn.nc','areacello');
    % loat areacello(x, y) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_MIROC-ES2H_historical_r2i1p4f2_gn.nc','thkcello');
    Cell_Thick = permute(thkcello(:,:,1:63),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
    %Cell_Thick = Cell_Volume./permute(repmat(Cell_Area,[1,1,63]),[3,1,2]); 
    Cell_Thick(isnan(phyc(1,:,:,:)))=nan;

    % Depth Level
    depthlev = ncread('volcello_Ofx_MIROC-ES2H_historical_r3i1p4f2_gn.nc','lev'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = double(lon.var.longitude); Lat= double(lat.var.latitude);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff','-append')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zooc    = zooc*1000;            % mol/m3   -> mmol/m3
    phyc    = phyc*1000;            % mol/m3   -> mmol/m3
    intpp   = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100  = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos     = tos;                  % degC     -> degC

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %zooc(zooc<1e-50)=nan;       
    zooc(zooc>10e15) = nan; 
    %phyc(phyc<1e-50)=nan;       
    phyc(phyc>10e15) = nan; 
    %intpp(intpp<1e-5)=nan;     
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;  
    epc100(epc100>10e15) = nan; 
    tos(tos>10e15) = nan; 

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    zooc    = circshift(zooc,-(lonoff-1),4);
    phyc    = circshift(phyc,-(lonoff-1),4);
   
    % 3 dimensions
    intpp   = circshift(intpp,-(lonoff-1),3); 
    epc100  = circshift(epc100,-(lonoff-1),3);
    tos     = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'zooc', 'phyc', 'intpp', 'epc100', '-append')
    save(fname, 'Cell_Thick','Cell_Area','Cell_Volume','-append');
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
    % PRIMARY Source:  Hajima et al. (2020); hajimaDevelopmentMIROCES2LEarth2020      
                    % 2P1Z with single-prey grazing
                      % P Groups: Nondiazatrophs, Diazatrophs
                    % N Units
                      % C:N:P = 106:16:1 (takahashiRedfieldRatioBased1985)
                    % No Temperature limiation
   % Parameters in:  Schmittner et al 2005; schmittnerGlobalModelMarine2005
                     % Via Table B2 in Appendix B
   % Equations in:   Schmittner et al 2005; schmittnerGlobalModelMarine2005
                     % Via Table B2 in Appendix B 
                     % Eq 14 of schmittnerGlobalModelMarine2005
                     % Does not give seperate equestions for diaz and phy
                     %      but is provided as reference for both
                     
   % Zooplankton Group 1: "Zooplankton"  Grazing On:
   diazfrac = 2./100; % Asssumes that only 2% of all phytoplankton are diaztrophs based on CESM 2 mean output, which is the closest model in structure with a diazatrophs group. 
   
          % "Phytoplankton (nondiazatrophs)"
          % D-III w/ single prey repsonse   
          %  g   = 2 (1/d); 
          %  eps = 1 (1/(mmol N m3)^2/d)
          %  K   =  sqrt(2/1) * 106/16 = 9.3692 (mmol C/m3) 
             g_Z_nonDiaz = 2.* ((phyc.*(1-diazfrac)).^2 ./ ((phyc.*(1-diazfrac)).^2 + 9.3692.^2));

        % "Diazatrophs"
          % D-III w/ single prey repsonse   
          %  g   = 2 (1/d); 
          %  eps = 1 (1/(mmol m3)^2/d)
          %  K   =  sqrt(2/1) * 106/16 = 9.3692 (mmol C/m3) 
             g_Z_diaz = 2.* ((phyc.*(diazfrac)).^2 ./ ((phyc.*(diazfrac)).^2 + 9.3692.^2));

    % Total Grazing 
        GSP = g_Z_nonDiaz.*zooc +  g_Z_diaz.*zooc;  

    % Total Trophic Transfer  
        g = 0.75; % Assimiliation efficiecny 
        GSP_ht = g.*GSP;
        
    % Bul Zooplankton Specific Grazing Rates  
        graze_spec = GSP./zooc;
        
    % Bulk mortality 
      Bz = 0.03; % 1/d Excretion of zooplankton
      Mz = 0.03;  % Mortalit of zooplankton 
      z_loss       =  Bz*zooc + .03.*zooc.^2; % .03 converted from 0.2  1/(mmolN/m)3 *d
      z_loss_spec  =  z_loss./zooc;
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP','GSP_ht', 'z_loss','z_loss_spec', '-append')
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

 % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phyC_inv    = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
         phyC_inv(isnan(squeeze(phyc(:,1,:,:)))) =nan;
    zooC_inv    = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
         zooC_inv(isnan(squeeze(phyc(:,1,:,:)))) =nan;
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_inv(isnan(squeeze(phyc(:,1,:,:)))) =nan;
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv(isnan(squeeze(phyc(:,1,:,:)))) =nan;
      
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
      
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phyC_inv',      'zooC_inv','gsp_inv', 'gsp_ht_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CESM2 - MARBL (USA) %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_CESM2_historical_r1i1p1f1_gn
% Variables Used: phydiat |phydiaz | phypico | zooc | intpp | epc100 |
%                 tos | areacello | volcello
% Caveats: A depth and temperature dependent threshold was  used for
% computing available phytoplantkon biomass --- from PC

% Notes:  - Different phenology than all other models with much more spring
%           centri bloom. (no drift though)
%         - Doesnt appear to agree with Fig14 of long et 2021 but that is
%         chl not carbon or npp and i cant find any apparent problems in
%         the processing. Plus total NPP seems fine. Maybe just dont focus
%         on phenology?

fname = 'CESM2';

for Compute = 1
fprintf('\n Starting CESM2 - MARBL (USA) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load('CESM2_Clims.mat')  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    % Volume
    volcello = ncread('volcello_Ofx_CESM2_historical_r11i1p1f1_gn.nc','volcello');
    % loat volcello(nlon, nlat, lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Volume = permute(volcello(:,:,1:15),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)

    % Area
    areacello = ncread('areacello_Ofx_CESM2_historical_r11i1p1f1_gn.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    Cell_Thick = Cell_Volume./permute(repmat(Cell_Area,[1,1,15]),[3,1,2]);
 
    % Depth Level
    depthlev = ncread('volcello_Ofx_CESM2_historical_r11i1p1f1_gn.nc','lev'); 
    depthlev = depthlev(1:15)./100; % Convert from cm to m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = lon.var.lon; Lat= lat.var.lat;
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff', '-append')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zooc    = zooc*1000;            % mol/m3   -> mmol/m3
    phydiat = phydiat*1000;         % mol/m3   -> mmol/m3
    phydiaz = phydiaz*1000;         % mol/m3   -> mmol/m3
    phypico = phypico*1000;         % mol/m3   -> mmol/m3
    phyc    =   phydiat + phydiaz ...
              + phypico;
    intpp   = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m2/d
    epc100  = epc100*1000*60*60*24; % mol/m2/s -> mmol/m2/d
    tos     = tos;                  % degC     -> degC
    thetao  = thetao(:,1:15,:,:);    % Top 150 m of potential T

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %zooc(zooc<1e-15)=nan;       
    zooc(zooc>10e15) = nan; 
    %phyc(phyc<1e-5)=nan;        
    phyc(phyc>10e15) = nan; 
    %phydiat(phydiat<1e-5)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phydiaz(phydiaz<1e-5)=nan;  
    phydiaz(phydiaz>10e15) = nan;
    %phypico(phypico<1e-5)=nan;  
    phypico(phypico>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;   
    epc100(epc100>10e15) = nan; 
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan;  

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    zooc    = circshift(zooc,-(lonoff-1),4);
    phyc    = circshift(phyc,-(lonoff-1),4);
    phydiat = circshift(phydiat,-(lonoff-1),4);
    phydiaz = circshift(phydiaz,-(lonoff-1),4);
    phypico = circshift(phypico,-(lonoff-1),4);
    thetao = circshift(thetao,-(lonoff-1),4);
    
    % 3 dimensions
    intpp   = circshift(intpp,-(lonoff-1),3); 
    epc100  = circshift(epc100,-(lonoff-1),3);
    tos     = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname,  'depthlev','zooc', 'phyc', 'phydiat', 'phydiaz', 'phypico', 'intpp', 'epc100', 'tos','thetao', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
    % PRIMARY Source: Long et al (2021);  longSimulationsMarineBiogeochemistry2021
                    % 3P1Z with single-prey grazing
                      %  P: Diatoms, Diazatrophs, Small pico/nano, (implicit calficiers)
                    % Carbon units (for Zooplankton)
                      % C, P, Si, Fe, CaCO3  units (for Phytoplankton)
                    
    % Equations in:   Long et al (2021);  longSimulationsMarineBiogeochemistry2021
    % Parameters in:  Long et al (2021);  longSimulationsMarineBiogeochemistry2021

    % Temp limiation and thresholds
    Tf = 1.7.^((thetao - 30)/10);    
    %Tf = permute(repmat(Tf,[1,1,1,15]),[1,4,2,3] );  % Use surface temeprature for limitation 
    
    % Phytopankton loss threshols
           % Temperature and Depth Dependent threshold. Only keeps seed population
           % if warm/shallow. Otherwise can go to zero.
                % Sources:
                % Equations
                %        https://github.com/marbl-ecosys/MARBL/blob/6b2a5e41dc5b0faaf9b600c96432e2487515e878/src/marbl_interior_tendency_mod.F90#L1270
                % Temp threshold parameters
                %       https://github.com/marbl-ecosys/MARBL/blob/stable/src/marbl_settings_mod.F90#L2885
                % loss threshold parameters; depth threshold parameters
                %      https://marbl.readthedocs.io/_/downloads/en/latest/pdf/
                %   

           % Depth Functionality, 
           zt1  = 80;           % (m) Pthresh = C_loss thres for z shallwer than this (cm) 
           zt2  = 120;          % (m) pthresh = 0 for z deeper than this (cm)
           % Linear interpolatation inbetween
           depth_grid = permute(repmat(depthlev(depthlev<dth,:,:),[1,384, 320, 12]),[4,1,2,3]);
           f_loss_thresh =  min(max((zt2-depth_grid)./(zt2-zt1),0),1);  
           % 5         15        25        35        45        55        65        75        85        95        105       115            125      135   145
           % 1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    0.8750    0.6250    0.3750    0.1250         0         0    0     

           sp_lt2 = 0;          % (mmolC/m3) C_loss thres when Temp is below themp thresh
           sp_lt1 = 0.01;       % (mmolC/m3) C_loss thres when Temp is above themp thresh
           sp_temp_th = -10;    % No P_threshold if temp is beloe temp thersholf
           % Pico Phy Threshold below  Temp thresh
           pthresh_sp = ones(size(phypico)).*sp_lt2.*f_loss_thresh;
           % Pico Phy Thresholf above Temp thresh
           warm_thr  = sp_lt1.*f_loss_thresh;
           pthresh_sp(thetao>sp_temp_th) =  warm_thr(thetao>sp_temp_th);

           diat_lt2 = 0;        % (mmolC/m3) C_loss thres when Temp is below themp thresh
           diat_lt1 = 0.02;     % (mmolC/m3) C_loss thres when Temp is above themp thresh
           diat_temp_th = -10;  % No P_threshold if temp is beloe temp thersholf
           % Diat Phy Thresholf below Temp thresh
           pthresh_diat = ones(size(phydiat)).*diat_lt2.*f_loss_thresh;
           % Diat Phy Thresholf above Temp thresh
           warm_thr  = diat_lt1.*f_loss_thresh;
           pthresh_diat(thetao>diat_temp_th) =  warm_thr(thetao>diat_temp_th);

           diaz_lt2 = 0.001;    % (mmolC/m3) C_loss thres when Temp is below themp thresh
           diaz_lt1 = 0.01;     % (mmolC/m3) C_loss thres when Temp is above themp thresh
           diaz_temp_th = 15;   % No P_threshold if temp is beloe temp thersholf
           % Diaz Phy Thresholf below Temp thresh
           pthresh_diaz = ones(size(phydiaz)).*diaz_lt2.*f_loss_thresh;
           % Diaz Phy Thresholf above Temp thresh
           warm_thr  = diaz_lt1.*f_loss_thresh;
           pthresh_diaz(thetao>diaz_temp_th) =  warm_thr(thetao>diaz_temp_th);

    % Zooplankton loss thresholds
    % Depth Functionality, 
           zt1_z  = 110;           % (m) Pthresh = C_loss thres for z shallwer than this (cm) 
           zt2_z  = 150;           % (m) pthresh = 0 for z deeper than this (cm)
         
            z_loss_thresh =  min(max((zt2_z-depth_grid)./(zt2_z-zt1_z),0),1); 
            z_thresh = z_loss_thresh * .075;
            zstar     =  max(0, (zooc-z_thresh));

    % Zooplankton Group 1 (Generic Zooplankton)
    % Specific Grazing Rates on:
      % Diatoms --  K: 1.2 (mmol C/m3); g: 3.3 @ 30 Celcius (Tref); EQ: g = g*Tf*(P'/P'+ K)
        pstar     =  max(0, (phydiat-pthresh_diat));
        g_Z_diat  = (3.3.*Tf.*pstar) ./(pstar+1.2); % 1/d 

      % Diazatrophs -- K: 1.2 (mmol C/m3); g: 3.15 @ 30 Celcius (Tref); EQ: g = g*Tf*(P'/P'+ K)
        pstar     =  max(0, (phydiaz-pthresh_diaz));
        g_Z_diaz  = (3.15.*Tf.*pstar)./(pstar+1.2); % 1/d 

      % Small Phyto -- K: 1.2 (mmol C/m3); g: 3.3 @ 30 Celcius (Tref) EQ: g = g*Tf*(P'/P'+ K)
        pstar     =  max(0, (phypico-pthresh_sp));
        g_Z_sp    = (3.3.*Tf.*pstar)./(pstar+1.2); % 1/d 

    % Total Grazing 
        GSP = g_Z_diat.*zooc +  g_Z_diaz.*zooc + g_Z_sp.*zooc;  

    % Total Trophic Transfer 
        gdiat = 0.25; gdiaz = 0.3; gsp = 0.3;  % Gross growth efficiencies
        GSP_ht = gdiat.*g_Z_diat.*zooc +  gdiaz.*g_Z_diaz.*zooc + gsp.*g_Z_sp.*zooc;  

    % Bul Zooplankton Specific Grazing Rates  
        graze_spec = GSP./zooc;
        
        
   % Bulk mortality 
      Bz = .1.*Tf; % Linear mortality rate
      Mz = .4.*Tf; % Quadratic Mortality Rate
      z_loss       =  Bz.*zstar + Mz.*zstar.^(1.5); 
      z_loss_spec  =  z_loss./zooc;
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP','GSP_ht', 'z_loss','z_loss_spec', '-append')
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv = squeeze(sum(phydiat(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phydiaz_inv = squeeze(sum(phydiaz(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    physp_inv   = squeeze(sum(phypico(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phyC_inv    = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv    = squeeze(sum(zooc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
   
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
                     
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
      
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phydiaz_inv',   'physp_inv',  'phyC_inv',   'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv', 'gsp_inv', 'gsp_ht_inv',...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CANESM5 - CanOE (CAN) %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_CanESM5-CanOE_historical_r2i1p2f1_gn
% Variables Used: phydiat |phymisc | zoomicro | zoomacro | intpp | epc100 |
%                 tos | areacello | volcello
% Caveats: None

fname = 'CanESM5-CanOE';

for Compute = 1
fprintf('\n Starting CANESM5 - CanOE (CAN) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
 load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    % Area
    areacello = ncread('areacello_Ofx_CanESM5-CanOE_historical_r2i1p2f1_gn.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_CanESM5-CanOE_historical_r2i1p2f1_gn.nc','thkcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Thick = permute(thkcello(:,:,1:45),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
   
     % Volume
    Cell_Volume = Cell_Thick.*permute(repmat(Cell_Area,[1,1,45]),[3,1,2]); % Unit: m^3, Dim: (level, lat, lon)
    
    % Depth Level
    depthlev = ncread('thkcello_Ofx_CanESM5-CanOE_historical_r2i1p2f1_gn.nc','lev'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = lon.var.longitude; Lat= lat.var.latitude;
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso  = zmeso*1000;         % mol/m3   -> mmol/m3
    zooc   = zmicro + zmeso;
    phydiat  = phydiat*1000;         % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;         % mol/m3   -> mmol/m3
    phyc     = phydiat + phymisc;
    phyn     = phyn*1000;         % mol/m3   -> mmol/m3
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC
    thetao   = thetao +273.15;       % deg C-> Deg K
    
    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %zooc(zooc<1e-15)=nan;        
    zooc(zooc>10e15) = nan;
    %phyc(phyc<1e-15)=nan;        
    phyc(phyc>10e15) = nan;
    phyn(phyn>10e15) = nan;
    phyfe(phyfe>10e15) = nan;
    %zmicro(zmicro<1e-15)=nan;    
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-15)=nan;      
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-15)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-15)=nan;  
    phymisc(phymisc>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;    
    epc100(epc100>10e15) = nan; 
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan; 
    
    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    phyc = circshift(phyc,-(lonoff-1),4);
    phyn = circshift(phyn,-(lonoff-1),4);
    phyfe = circshift(phyfe,-(lonoff-1),4);
    zooc = circshift(zooc,-(lonoff-1),4);
    zmicro    = circshift(zmicro,-(lonoff-1),4);
    zmeso = circshift(zmeso,-(lonoff-1),4);
    phydiat = circshift(phydiat,-(lonoff-1),4);
    phymisc = circshift(phymisc,-(lonoff-1),4);
    thetao = circshift(thetao,-(lonoff-1),4);
    
    
    % 3 dimensions
    intpp   = circshift(intpp,-(lonoff-1),3); 
    epc100  = circshift(epc100,-(lonoff-1),3);
    tos     = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname,  'depthlev', 'zmicro', 'zmeso', 'zooc', 'phydiat', 'phymisc', 'phyc','phyn','phyfe', 'intpp', 'epc100', 'tos', 'thetao', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
% PRIMARY Source: Christian et al. (2021);  christianOceanBiogeochemistryCanadian2021
                    % 2P2Z with multi-prey response 
                       % Phytoplankton: Small, Large
                       % Zooplankton:   Small, Large
                       % Detrius:       Small, Large
                    % C units
                      % Flexible C:N:Fe:Chl ratios in Phyto but grazing depends on C
                      % Just C in Zoo       
             
                    
 % Parameters in:  christianOceanBiogeochemistryCanadian2021
                   
 % Equations in:   christianOceanBiogeochemistryCanadian2021
                    
 
 % NOTES: More complex version of ocean ecosystems, however many cmip6 experiments 
%         used previous BGcmodel: CMOC
  
  % Zoo. Grazing Eqautions
 

        % Z1: "Small Zooplankton" Grazing On
        % Prey: Small Phyto

            % "Small Phytoplankton"
            %  Single-prey   
            %  g   =  1.7 (1/d); 
            %  y   =  0.25 1/(mmolC/m3); Ivlev Parameter
            %  K   =  2.77 (mmolC/m3); Converted from -ln(.5)/.y
             g_Zmicro_misc = 1.7.*(1-exp(-0.25.*(phymisc)));

        % Z2: "Large Zooplankton" Grazing On
        % Prey: Large Phyto, Small Zooplantkon

            % "Large Phytoplankton"
            %  Single-prey w/ switching  
            %  g   =  0.85 (1/d); 
            %  y   =  0.25 1/(mmolC/m3); Ivlev Parameter
            %  K   =  2.77 (mmolC/m3); Converted from -ln(.5)/.y
             g_Zmeso_diat = 0.85.*(1-exp(-0.25.*(phydiat+zmicro))).*(phydiat./(phydiat+zmicro));

            % "Small Zooplankton"
            %  Single-prey w/ switching  
            %  g   =  0.85 (1/d); 
            %  y   =  0.25 1/(mmolC/m3); Ivlev Parameter
            %  K   =  2.77 (mmolC/m3); Converted from -ln(.5)/.y
             g_Zmeso_Zmicro = 0.85.*(1-exp(-0.25.*(phydiat+zmicro))).*(zmicro./(phydiat+zmicro));



    % Total Grazing 
        GSP = g_Zmicro_misc.*zmicro +  g_Zmeso_diat.*zmeso;  

    % Total Trophic Transfer 
        g = .8; % Assimiliation efficiency for all grazing terms
        GSP_ht = g.*GSP; 
        
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
    
   % Bulk mortality 
   % Respiration 
       rzs = .3; rzl = .1; % Spefic Raspiration rates
       RCN=106/16; 
       g_prime_micro = g_Zmicro_misc.*zmicro .* (phyn./phyc).*RCN;
       g_prime_meso  = g_Zmeso_diat.*zmeso   .* (phyn./phyc).*RCN;
       CXS_micro     =  g_prime_micro .* ((phyc./phyn).*1/RCN -1); 
       CXS_meso      =  g_prime_meso .* ((phyc./phyn).*1/RCN -1); 
       
    
       a = -37.4e3./8.314 ;%(j/mol/K)
       tf =  exp(a.*(1./thetao -1./298.15) ); % T in Kelvin
       rs = max(rzs.*tf.*zmicro-CXS_micro,0);
       rl = max(rzl.*tf.*zmeso-CXS_meso,0);
      
      z_loss       =  .05.*zmicro + .06.*zmicro.^2 +rs ...
                     + .1.*zmeso+ .06.*zmeso.^2 +rl ...
                     + (1-g) .* g_Zmeso_Zmicro.*zmeso;      % Biomass lost from total zooplankton due to assimiliation effciency during carnivory
      % Note, both Z pools also include explict respiration term
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP_ht', 'GSP', 'z_loss','z_loss_spec', '-append')
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    phyC_inv      = squeeze(sum(phyc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)   .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv      = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
   
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
       
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
  
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',    'graze_inv','gsp_ht_inv','gsp_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% UKESM1 - MEDUSA2.1 (UK) %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_UKESM1-0-LL_historical_r9i1p1f2_gn
% Variables Used: phydiat |phymisc | zoomicro | zoomacro | intpp | epc100 | detoc
%                 tos | areacello | volcello
% Caveats: - % indexed in chuks of 12 for climatology instead of using datetime stamp which appears worng... % Previously was Only 10 year climatology to avoid sesonal drift
%          -Cell Thickness and Volume are prognostic. We use the 50 year mean field
%               -Does not appear to have much seaonal variability but could
%                 make climatology instead
% Note:    - Depth integrated and Spatially integrated values for detritus
%            are not saved
%          - GSP is only Zooplankton consumption of primary producers (not
%            other Z, bac,  or detritus)
%          - Explicitly saved graze seems way too high, much higher than
%          npp... and well different from diagnostic one...

fname = 'UKESM1-0-LL';

for Compute = 1
fprintf('\n Starting UKESM1 - MEDUSA2.1 (UK) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

    % Volume
%     volcello = ncread('volcello_Omon_UKESM1-0-LL_historical_r9i1p1f2_gn_185001-189912.nc','volcello');
%     % load volcello(nlon, nlat, lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
%     volcello = squeeze(nanmean(volcello(:,:,:,1:120),4)); % only average first 10 years if using 10 year climatology
%     Cell_Volume = permute(volcello(:,:,1:75),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)
    load UKESM_clim_CellVol_10yr Cell_Volume % avoid load massive time varying file!
     
    % Area
    %areacello = ncread('areacello_Ofx_UKESM1-0-LL_piControl_r1i1p1f2_gn.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    %Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)
    load UKESM_clim_CellVol_10yr Cell_Area % avoid load massive time varying file!
  
    % Thickness
    % thkcello = ncread('thkcello_Omon_UKESM1-0-LL_historical_r7i1p1f3_gn_185001-189912.nc','thkcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    % Cell_Thick = permute(thkcello(:,:,1:75),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
    Cell_Thick = Cell_Volume./permute(repmat(Cell_Area,[1,1,75]),[3,1,2]); 
    
    % Depth Level
    depthlev = cumsum(nanmean(nanmean(Cell_Thick,2),3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = double(lon.var.longitude); Lat= double(lat.var.latitude);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro   = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso    = zmeso*1000;         % mol/m3   -> mmol/m3
    phydiat  = phydiat*1000;       % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;       % mol/m3   -> mmol/m3
    phyc     = phydiat+phymisc;
    zooc     = zmicro+zmeso;
    detoc    = detoc*1000;         % mol/m3   -> mmol/m3
    
    graz     = graz*1000*60*60*24;  % mol/m3/s -> mmol/m3/d
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
     %zooc(zooc<1e-15)=nan;        
    zooc(zooc>10e15) = nan;
    %phyc(phyc<1e-15)=nan;        
    phyc(phyc>10e15) = nan;
    %zmicro(zmicro<1e-15)=nan;    
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-15)=nan;      
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-15)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-15)=nan;  
    phymisc(phymisc>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;    
    epc100(epc100>10e15) = nan; 
    %detoc(detoc<1e-5)=nan;         
    detoc(detoc>10e15) = nan;
    %graz(graz<1e-5)=nan;           
    graz(graz>10e15) = nan;
    tos(tos>10e15) = nan; 

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    phyc = circshift(phyc,-(lonoff-1),4);
    zooc = circshift(zooc,-(lonoff-1),4);
    zmicro = circshift(zmicro,-(lonoff-1),4);
    zmeso  = circshift(zmeso,-(lonoff-1),4);
    phydiat  = circshift(phydiat,-(lonoff-1),4);
    phymisc  = circshift(phymisc,-(lonoff-1),4);
    detoc    = circshift(detoc,-(lonoff-1),4);
    graz     = circshift(graz,-(lonoff-1),4); 
    
    % 3 dimensions
    intpp   = circshift(intpp,-(lonoff-1),3); 
    epc100  = circshift(epc100,-(lonoff-1),3);
    tos     = circshift(tos,-(lonoff-1),3);
    %Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);   % Grid has already been shifted!
    %Cell_Thick  = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    %Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'phyc','zooc', 'zmicro', 'zmeso', 'phydiat', 'phymisc', 'detoc', 'intpp', 'epc100','graz', 'tos', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
% PRIMARY Source: Yool et al (2021);  yoolEvaluatingPhysicalBiogeochemical2021
                    % 2P2Z with multi-prey response w/ switching (i.e. Fasham et al, 1990)
                       % Phytoplankton: Non-diatoms (nanophytoplankton), Diatoms (microphytoplankton)
                       % Zooplankton:   Microzooplankton,                Macrozooplankton
                       % Detrius:       Small (explicit),                Large (implict)
                    % N units
                      % Fixed C:N ratios in Phyto, Zoo
                      % Variable C:N Det as it is fed by phyto and zoo and zoo have lower C:N
                      % 6.625 = 106:16 C:N (phyto), 5.625 (zoo)       
                    
 % Parameters in:  Yool et al (2013);  yoolMEDUSA2IntermediateComplexity2013
                   
 % Equations in:   Yool et al (2013);  yoolMEDUSA2IntermediateComplexity2013
                    
 
 % NOTES: Medusa 2.1 (Yool et al., 2021) evolved from Medusa 2.1 (Yool et al., 2013)
 %        HadGEM2-ES2 (w/ bgc from diat-HadOCC) is insitiutional
 %        predecessor to the UKESM1 (w/ bgc from MEDUSA. 
  

                       
        % Z1: "microooplankton"  Grazing On
        % Prey: nondiatoms, Detritus (slow)

        % Preferences
            % Base Preferencea
            prf_det      = 0.25;
            prf_misc     = 0.75;

            % Total apparent food
            food = prf_misc.*phymisc.^2 + prf_det.*detoc.^2;

            %  phymisc: "Nondiatoms (nanophytoplankton)"
            %  Multi-prey w/ switching  
            %  g   =  2.0 (1/d); 
            %  K   =  5.3 (mmolC/m3); Converted from 0.8 mmol N/m3
             g_Zmicro_phymisc = (2.*prf_misc.*phymisc.^2) ./ (5.3.^2 + food);

                                  
        % Z2: "mespooplankton"  Grazing On
        % Prey: nondiatoms, Diatoms, microzooplankton, Detritus (slow)
        % Preferences

            % Base Preferences
            prf_det      = 0.15;
            prf_misc     = 0.15;
            prf_diat     = 0.35;
            prf_zmicro   = 0.35;

            % Total apparent food
            food = prf_misc.*phymisc.^2 + prf_det.*detoc.^2 + prf_diat.*phydiat.^2 + prf_zmicro.*zmicro.^2;

            % "Nondiatoms (nanophytoplankton)"
            %  Multi-prey w/ switching  
            %  g   =  0.5 (1/d); 
            %  K   =  1.9875 (mmolC/m3); Converted from 0.3 mmol N/m3
             g_Zmeso_phymisc = (0.5.*prf_misc.*phymisc.^2) ./ (1.9875.^2 + food);

            % "Diatoms (micro)"
            %  Multi-prey w/ switching  
            %  g   =  0.5 (1/d); 
            %  K   =  1.9875 (mmolC/m3); Converted from 0.3 mmol N/m3
             g_Zmeso_phydiat = (0.5.*prf_diat.*phydiat.^2) ./ (1.9875.^2 + food);

             % "MicroZ "
            %  Multi-prey w/ switching  
            %  g   =  0.5 (1/d); 
            %  K   =  1.9875 (mmolC/m3); Converted from 0.3 mmol N/m3
             g_Zmeso_Zmicro = (0.5.*prf_zmicro.*zmicro.^2) ./ (1.9875.^2 + food);
             
    % Total Grazing 
        GSP = g_Zmicro_phymisc.*zmicro +  g_Zmeso_phymisc.*zmeso +  g_Zmeso_phydiat.*zmeso;  

    % Total Trophic Transfer
      g= .8.*.8.*.2; 
      GSP_ht = g.*GSP; 
      
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
        
    % Total Grazing -- Explicit Calculation    
        GSP_exp = graz;
        
    % Bulk Zooplankton Specific Grazing Rates -- Explict Calculation
        graze_spec_exp = GSP_exp./(zmicro+zmeso);
     
   % Bulk mortality 
      z_loss       =  .02.*zmicro + .1.*(zmicro./(2.81  + zmicro)).*zmicro   ...
                   +  .02.*zmeso  + .2.*(zmeso./(4.22    + zmeso)).*zmeso ...
                   +   (1-g).* g_Zmeso_Zmicro.*zmeso;            % K parameter converted from N using zooplankton c:N of 5.625
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP_ht', 'GSP','graze_spec_exp', 'GSP_exp', 'z_loss','z_loss_spec',   '-append')
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    phyC_inv      = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)   .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv      = squeeze(sum(zooc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    detoc_inv     = squeeze(sum(detoc(:,depthlev<dth,:,:)   .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
     
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_inv_exp = squeeze(sum(GSP_exp(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv   = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
       
    % Specific Rates (1/d)
    div_inv     = npp_inv./phyC_inv;
    graze_inv   = gsp_inv./zooC_inv; 
    graze_inv_exp   = gsp_inv_exp./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
    
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',     'graze_inv', 'detoc_inv', 'graze_inv_exp','gsp_inv_exp','gsp_inv', 'gsp_ht_inv',...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IPSL-CM6A PISCESv2 (FRA) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_IPSL-CM6A-LR_historical_r10i1p1f1_gn
% Variables Used: phydiat |phymisc | zoomicro | zoomacro | intpp | epc100 | detoc
%                 tos | areacello | volcello
% Caveats: 
% Note:    Depth integrated and Spatially integrated values for detritus
%          are not saved
%          GSP is only Zooplankton consumption of primary producers (not
%          other Z, bac,  or detritus)

fname = 'IPSL-CM6A-LR';

for Compute = 1
fprintf('\n Starting IPSL-CM6A-LR - PISCESv2 (FRA) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
 load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

   
    % Area
    areacello = ncread('areacello_Ofx_IPSL-CM6A-LR_historical_r32i1p1f1_gn.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Odec_IPSL-CM6A-LR_historical_r24i1p1f1_gn_1855-2005.nc','thkcello');
    thkcello_m =nanmean(thkcello,4);
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Thick = permute(thkcello_m(:,:,1:75),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
   
    % Volume
    Cell_Volume = Cell_Thick.*permute(repmat(Cell_Area,[1,1,75]),[3,1,2]); % Unit: m^3, Dim: (level, lat, lon)
    
     % Depth Level
    depthlev = ncread('thkcello_Odec_IPSL-CM6A-LR_historical_r24i1p1f1_gn_1855-2005.nc','olevel'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = double(lon.var.nav_lon); Lat= double(lat.var.nav_lat);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro   = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso    = zmeso*1000;         % mol/m3   -> mmol/m3
    phydiat  = phydiat*1000;       % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;       % mol/m3   -> mmol/m3
    phyc     = phydiat+phymisc;
    zooc     = zmicro+zmeso;
    detoc    = detoc*1000;         % mol/m3   -> mmol/m3
    graz     = graz*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC
    o2       = o2.*1000;             % mol/m3 ->mmol/m3
   
     
    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
     %zooc(zooc<1e-15)=nan;        
    zooc(zooc>10e15) = nan;
    %phyc(phyc<1e-15)=nan;        
    phyc(phyc>10e15) = nan;
    %zmicro(zmicro<1e-15)=nan;    
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-15)=nan;      
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-15)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-15)=nan;  
    phymisc(phymisc>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;    
    epc100(epc100>10e15) = nan; 
    %detoc(detoc<1e-5)=nan;         
    detoc(detoc>10e15) = nan;
    %graz(graz<1e-5)=nan;           
    graz(graz>10e15) = nan;
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan;   
    o2(o2>10e15) = nan; 
    
    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    phyc     = circshift(phyc,-(lonoff-1),4);
    zooc     = circshift(zooc,-(lonoff-1),4);
    zmicro   = circshift(zmicro,-(lonoff-1),4);
    zmeso    = circshift(zmeso,-(lonoff-1),4);
    phydiat  = circshift(phydiat,-(lonoff-1),4);
    phymisc  = circshift(phymisc,-(lonoff-1),4);
    detoc    = circshift(detoc,-(lonoff-1),4);
    graz     = circshift(graz,-(lonoff-1),4); 
    thetao   = circshift(thetao,-(lonoff-1),4);
    o2   = circshift(o2,-(lonoff-1),4);
    
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    tos          = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick   = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname,  'depthlev','phyc','zooc', 'zmicro', 'zmeso', 'phydiat', 'phymisc', 'detoc', 'intpp', 'epc100', 'graz', 'tos','thetao', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume','o2', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
 % PRIMARY Source: Aumont et al. (2015);  aumontPISCESv2OceanBiogeochemical2015
                    % 2P2Z with multi-prey response 
                       % Phytoplankton: Nanophytoplankton,  Diatoms
                       % Zooplankton:   Micro,              Meso
                       % Detrius:       POC
                    % C units
                      % Just C in Zoo       
             
                    
 % Parameters in:  aumontPISCESv2OceanBiogeochemical2015
                   
 % Equations in:   aumontPISCESv2OceanBiogeochemical2015
                    
% Zooplantkton Grazing Formulaton
         
% Auxilliary Variables
 
    % Class dependent T-lim terms scale the max class
    %       specific grazing rate at 0 C
     b_micro = 1.079;  b_meso = 1.079; b_zoo =1.079;
     ft_micro = b_micro.^thetao;
     ft_meso = b_meso.^thetao;
     ft_z = b_zoo.^thetao;
     
    % Phyotplankton Threshold 
     th_prey     = 0.001;
     th_glob     = 0.3;


% Z1: "MicroZooplankton"  Grazing On:
%     Prey: Nanophytoplankton, Diatoms, POC

    % Base Preferencea
    prf_misc     = 1;
    prf_diat     = 0.5;
    prf_poc      = 0.1;


    % Total apparent food
    food_th = prf_misc.* max(0,phymisc-th_prey) + ...
              prf_diat.* max(0,phydiat-th_prey) + ...
              prf_poc.* max(0,detoc-th_prey);

    food   = prf_misc.* phymisc + ...
             prf_diat.* phydiat + ...
             prf_poc.* detoc;

    food_lim = max(0, food_th-min(.5*food_th,th_glob)); 

    % "NanoPhytoplankton"
    %  Single-prey w/ switching  
    %  g   =  3.0 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmicro_phymisc = (3.0 .* ft_micro .* (food_lim./food) .* prf_misc .* max(0,phymisc-th_prey)) ...
                       ./(20 + food);

    % "Diatoms"
    %  Single-prey w/ switching  
    %  g   =  3.0 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmicro_phydiat = (3.0 .* ft_micro .* (food_lim./food) .* prf_diat .* max(0,phydiat-th_prey)) ...
                       ./(20 + food);


% Z2: "Mesozooplankton"  Grazing On:
%      Prey: Nanophytoplankton, Diatoms, POC, and microzooplanton

    % Base Preferencea
    prf_misc     = 0.3;
    prf_diat     = 1;
    prf_poc      = 0.3;
    prf_zmicro      = 1.0;


    % Total apparent food
    food_th = prf_misc.* max(0,phymisc-th_prey) + ...
              prf_diat.* max(0,phydiat-th_prey) + ...
              prf_poc.* max(0,detoc-th_prey) + ...
              prf_zmicro.* max(0,zmicro-th_prey);

    food   = prf_misc.* phymisc + ...
             prf_diat.* phydiat + ...
             prf_poc.* detoc + ...
             prf_zmicro.* zmicro;

    food_lim = max(0, food_th-min(.5*food_th,th_glob)); 

    % "NanoPhytoplankton"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_phymisc = (0.75 .* ft_meso .* (food_lim./food) .* prf_misc .* max(0,phymisc-th_prey)) ...
                       ./(20 + food);

    % "Diatoms"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_phydiat = (0.75 .* ft_meso .* (food_lim./food) .* prf_diat .* max(0,phydiat-th_prey)) ...
                       ./(20 + food);
                   
    % "Micro Zoo"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_Zmicro = (0.75 .* ft_meso .* (food_lim./food) .* prf_zmicro .* max(0,zmicro-th_prey)) ...
                       ./(20 + food);
         
    % Total Grazing 
        GSP = g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro + ...
               g_Zmeso_phymisc.*zmeso  +   g_Zmeso_phydiat.*zmeso;  

    % Total Trophic Transfer (ingestion...but not respiration) 
      gmicro = .3; gmeso = .35;
      GSP_ht = gmicro.*(g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro )+ ...
               gmeso.*(g_Zmeso_phymisc.*zmeso  +   g_Zmeso_phydiat.*zmeso);  
      
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
        
    % Total Grazing -- Explicit Calculation    
        GSP_exp = graz;
        
    % Bulk Zooplankton Specific Grazing Rates -- Explict Calculation
        graze_spec_exp = GSP_exp./(zmicro+zmeso);
     
   % Bulk mortality 
    % Oxygen Degredation Factor -- Mortality is enhanced when oxygen is depleted. In other words, microzooplankton (but also mesozooplankton,see below) are treated as being unable to cope with anoxic 
                                % waters. This increased mortality also avoids respiration in waters devoid of oxygen.
     
     DO2 = 3.*min(1,max(0,0.4.*(1-o2)./(6+o2))); % 1 and 6 are the vlaues of o2min,1 and o2 min2 in units of umolO2/L
   
      z_loss       =  ft_z.*.03.* (zmicro./(.2+zmicro) + DO2).* zmicro + ft_z.*.004.*zmicro.^2  ...
                     + (1-gmeso).* g_Zmeso_Zmicro.*zmeso ...
                     + ft_z.*.005.*(zmeso./(.2+zmeso)  + DO2).* zmeso+ ft_z.*.03.*zmeso.^2; 
           
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP_ht', 'GSP','graze_spec_exp', 'GSP_exp',  'z_loss','z_loss_spec',  '-append')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    phyC_inv      = squeeze(sum(phyc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)   .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv      = squeeze(sum(zooc(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    detoc_inv   = squeeze(sum(detoc(:,depthlev<dth,:,:) .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
     
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;   
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv  = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_inv_exp = squeeze(sum(GSP_exp(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    
    % Specific Rates (1/d)
    div_inv       = npp_inv./phyC_inv;
    graze_inv     = gsp_inv./zooC_inv; 
    graze_inv_exp = gsp_inv_exp./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate', 'GSP transfer'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)'};
       
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv','gsp_ht_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',     'graze_inv', 'detoc_inv', 'graze_inv_exp','gsp_inv_exp','gsp_inv',...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CNRM-ESM2.1 PISCESv2 (FRA) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_CNRM-ESM2-1_historical_r1i1p1f2_gn
% Variables Used: phydiat |phymisc | zoomicro | zoomacro | intpp | epc100 | detoc
%                 tos | areacello | volcello
% Caveats: -% Core inv variables are capped at the 99.9 percentile to remove strange
%             export/P biomass event off eat coast of russia which is an outlier biassing seasonal cycle
% Note:    % New version of IPSL i believe!
%          Depth integrated and Spatially integrated values for detritus
%          are not saved
%          GSP is only Zooplankton consumption of primary producers (not
%          other Z, bac,  or detritus)

fname = 'CNRM-ESM2-1';

for Compute = 1
fprintf('\n Starting CNMR-ESM2.1 - PISCESv2 (FRA) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
 load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

   
    % Area
    areacello = ncread('areacello_Ofx_CNRM-ESM2-1_historical_r2i1p1f2_gn.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_CNRM-CM6-1_dcppC-ipv-pos_r1i1p1f2_gn.nc','thkcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Thick = permute(thkcello(:,:,1:75),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
   
    % Volume
    Cell_Volume = Cell_Thick.*permute(repmat(Cell_Area,[1,1,75]),[3,1,2]); % Unit: m^3, Dim: (level, lat, lon)
    
      % Depth Level
    depthlev = ncread('thkcello_Ofx_CNRM-CM6-1_dcppC-ipv-pos_r1i1p1f2_gn.nc','lev'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = double(lon.var.lon); Lat= double(lat.var.lat);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro   = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso    = zmeso*1000;         % mol/m3   -> mmol/m3
    phydiat  = phydiat*1000;       % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;       % mol/m3   -> mmol/m3
    phyc     = phydiat+phymisc;
    zooc     = zmicro+zmeso;
    detoc    = detoc*1000;         % mol/m3   -> mmol/m3
    graz     = graz*1000*60*60*24;  % mol/m3/s -> mmol/m3/d
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC
    o2       = o2.*1000;             % mol/m3 ->mmol/m3
    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %zooc(zooc<1e-15)=nan;        
    zooc(zooc>10e15) = nan;
    %phyc(phyc<1e-15)=nan;        
    phyc(phyc>10e15) = nan;
    %zmicro(zmicro<1e-15)=nan;    
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-15)=nan;      
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-15)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-15)=nan;  
    phymisc(phymisc>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;    
    epc100(epc100>10e15) = nan; 
    %detoc(detoc<1e-5)=nan;         
    detoc(detoc>10e15) = nan;
    %graz(graz<1e-5)=nan;           
    graz(graz>10e15) = nan;
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan;  
    o2(o2>10e15) = nan;  
    
    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    phyc     = circshift(phyc,-(lonoff-1),4);
    zooc     = circshift(zooc,-(lonoff-1),4);
    zmicro   = circshift(zmicro,-(lonoff-1),4);
    zmeso    = circshift(zmeso,-(lonoff-1),4);
    phydiat  = circshift(phydiat,-(lonoff-1),4);
    phymisc  = circshift(phymisc,-(lonoff-1),4);
    detoc    = circshift(detoc,-(lonoff-1),4);
    graz     = circshift(graz,-(lonoff-1),4); 
    thetao   = circshift(thetao,-(lonoff-1),4);
    o2       = circshift(o2,-(lonoff-1),4);
   
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    tos          = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick   = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    

    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname,  'depthlev','phyc','zooc', 'zmicro', 'zmeso', 'phydiat', 'phymisc', 'detoc', 'intpp', 'epc100', 'tos', 'graz','thetao','o2', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
 % PRIMARY Source: Aumont et al. (2015);  aumontPISCESv2OceanBiogeochemical2015
                    % 2P2Z with multi-prey response 
                       % Phytoplankton: Nanophytoplankton,  Diatoms
                       % Zooplankton:   Micro,              Meso
                       % Detrius:       POC
                    % C units
                      % Just C in Zoo       
             
                    
 % Parameters in:  aumontPISCESv2OceanBiogeochemical2015
                   
 % Equations in:   aumontPISCESv2OceanBiogeochemical2015
                    
 
 % NOTES: More complex version of ocean ecosystems, however many cmip6 experiments 
%         used previous BGcmodel: CMOC
  
% Zooplantkton Grazing Formulaton
         
    % Class dependent T-lim terms scale the max class
    %       specific grazing rate at 0 C
     b_micro = 1.079;  b_meso = 1.079; b_zoo =1.079;
     ft_micro = b_micro.^thetao;
     ft_meso = b_meso.^thetao;
     ft_z = b_zoo.^thetao;
     
    % Phyotplankton Threshold 
     th_prey     = 0.001;
     th_glob     = 0.3;


% Z1: "MicroZooplankton"  Grazing On:
%     Prey: Nanophytoplankton, Diatoms, POC

    % Base Preferencea
    prf_misc     = 1;
    prf_diat     = 0.5;
    prf_poc      = 0.1;


    % Total apparent food
    food_th = prf_misc.* max(0,phymisc-th_prey) + ...
              prf_diat.* max(0,phydiat-th_prey) + ...
              prf_poc.* max(0,detoc-th_prey);

    food   = prf_misc.* phymisc + ...
             prf_diat.* phydiat + ...
             prf_poc.* detoc;

    food_lim = max(0, food_th-min(.5*food_th,th_glob)); 

    % "NanoPhytoplankton"
    %  Single-prey w/ switching  
    %  g   =  3.0 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmicro_phymisc = (3.0 .* ft_micro .* (food_lim./food) .* prf_misc .* max(0,phymisc-th_prey)) ...
                       ./(20 + food);

    % "Diatoms"
    %  Single-prey w/ switching  
    %  g   =  3.0 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmicro_phydiat = (3.0 .* ft_micro .* (food_lim./food) .* prf_diat .* max(0,phydiat-th_prey)) ...
                       ./(20 + food);


% Z2: "Mesozooplankton"  Grazing On:
%      Prey: Nanophytoplankton, Diatoms, POC, and microzooplanton

    % Base Preferencea
    prf_misc     = 0.3;
    prf_diat     = 1;
    prf_poc      = 0.3;
    prf_zmicro      = 1.0;


    % Total apparent food
    food_th = prf_misc.* max(0,phymisc-th_prey) + ...
              prf_diat.* max(0,phydiat-th_prey) + ...
              prf_poc.* max(0,detoc-th_prey) + ...
              prf_zmicro.* max(0,zmicro-th_prey);

    food   = prf_misc.* phymisc + ...
             prf_diat.* phydiat + ...
             prf_poc.* detoc + ...
             prf_zmicro.* zmicro;

    food_lim = max(0, food_th-min(.5*food_th,th_glob)); 

    % "NanoPhytoplankton"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_phymisc = (0.75 .* ft_meso .* (food_lim./food) .* prf_misc .* max(0,phymisc-th_prey)) ...
                       ./(20 + food);

    % "Diatoms"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_phydiat = (0.75 .* ft_meso .* (food_lim./food) .* prf_diat .* max(0,phydiat-th_prey)) ...
                       ./(20 + food);
                   
    % "Micro Zoo"
    %  Single-prey w/ switching  
    %  g   =  0.75 (1/d); 
    %  K   =  20 (mmolC/m3); Converted from 20 umol C/L
     g_Zmeso_Zmicro = (0.75 .* ft_meso .* (food_lim./food) .* prf_zmicro .* max(0,zmicro-th_prey)) ...
                       ./(20 + food);
         
    % Total Grazing 
        GSP = g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro + ...
               g_Zmeso_phymisc.*zmeso  +   g_Zmeso_phydiat.*zmeso;  

    % Total Trophic Transfer (ingestion...but not respiration) 
      gmicro = .3; gmeso = .35;
      GSP_ht = gmicro.*(g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro )+ ...
               gmeso.*(g_Zmeso_phymisc.*zmeso  +   g_Zmeso_phydiat.*zmeso);  
      
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
        
    % Total Grazing -- Explicit Calculation    
        GSP_exp = graz;
        
    % Bulk Zooplankton Specific Grazing Rates -- Explict Calculation
        graze_spec_exp = GSP_exp./(zmicro+zmeso);
     
   % Bulk mortality 
    % Oxygen Degredation Factor -- Mortality is enhanced when oxygen is depleted. In other words, microzooplankton (but also mesozooplankton,see below) are treated as being unable to cope with anoxic 
                                % waters. This increased mortality also avoids respiration in waters devoid of oxygen.
     
     DO2 = 3.*min(1,max(0,0.4.*(1-o2)./(6+o2))); % 1 and 6 are the vlaues of o2min,1 and o2 min2 in units of umolO2/L
   
      z_loss       =  ft_z.*.03.* (zmicro./(.2+zmicro) + DO2).* zmicro + ft_z.*.004.*zmicro.^2  ...
                     + (1-gmeso).* g_Zmeso_Zmicro.*zmeso ...
                     + ft_z.*.005.*(zmeso./(.2+zmeso)  + DO2).* zmeso+ ft_z.*.03.*zmeso.^2; 
           
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP_ht', 'GSP','graze_spec_exp', 'GSP_exp',  'z_loss','z_loss_spec',  '-append')

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
        phydiat_inv(phydiat_inv>prctile(phydiat_inv(:),99.9))=prctile(phydiat_inv(:),99.9);
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
        phymisc_inv(phymisc_inv>prctile(phymisc_inv(:),99.9))=prctile(phymisc_inv(:),99.9);
    phyC_inv      = phydiat_inv + phymisc_inv;
                  %phyC_inv(phyC_inv>598)=598; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
                  zmicro_inv(zmicro_inv>prctile(zmicro_inv(:),99.9))=prctile(zmicro_inv(:),99.9);
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
                  zmeso_inv(zmeso_inv>prctile(zmeso_inv(:),99.9))=prctile(zmeso_inv(:),99.9);
    zooC_inv      = zmicro_inv + zmeso_inv;
                 % zooC_inv(zooC_inv>653)=653; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    detoc_inv   = squeeze(sum(detoc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
                  detoc_inv(detoc_inv>prctile(detoc_inv(:),99.9))=prctile(detoc_inv(:),99.9);
                  
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
                npp_inv(npp_inv>216)=216; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    exp_inv     = epc100;   
                exp_inv(exp_inv>38.8)=38.8; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
                gsp_inv(gsp_inv>98.8)=98.8; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    gsp_ht_inv  = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
                gsp_ht_inv(gsp_ht_inv>98.8)=98.8; % cap exp at 99.9 percentile to remove strange export event off eat coast of russia
    gsp_inv_exp = squeeze(sum(GSP_exp(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    
    % Specific Rates (1/d)
    div_inv       = npp_inv./phyC_inv;
    graze_inv     = gsp_inv./zooC_inv; 
    graze_inv_exp = gsp_inv_exp./zooC_inv; 
    
    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv; 
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate',''};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)'};
    
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',     'graze_inv','detoc_inv','graze_inv_exp','gsp_inv_exp','gsp_inv',  ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end

 


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CMCC-ESM2  BFM5.2 (ITA) %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run: Omon_CMCC-ESM2_historical_r1i1p1f1_gn
% Variables Used: phydiat |phymisc | zoomicro | zoomacro | intpp | epc100 | bacc
%                 tos | areacello | volcello
% Caveats: 
% Note:    Depth integrated and Spatially integrated values for detritus
%          are not saved
%          GSP is only Zooplankton consumption of primary producers (not
%          other Z, bac,  or detritus)

fname = 'CMCC-ESM2';

for Compute = 1
fprintf('\n Starting CMCC-ESM2  BFM5.2 (ITA) \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
 load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

      % Volume
    volcello = ncread('volcello_Ofx_CMCC-ESM2_historical_r1i1p1f1_gn.nc','volcello');
    % load volcello(nlon, nlat, lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Volume = permute(volcello(:,:,1:50),[3,2,1]); % Unit: m^3, Dim: (level, lat, lon)
    
    % Area
    areacello = ncread('areacello_Ofx_CMCC-ESM2_historical_r1i1p1f1_gn.nc','areacello');   
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    Cell_Thick = Cell_Volume./permute(repmat(Cell_Area,[1,1,50]),[3,1,2]); 
     
   % Depth Level
    depthlev = ncread('volcello_Ofx_CMCC-ESM2_historical_r1i1p1f1_gn.nc','lev'); 
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    Lon = double(lon.var.longitude); Lat= double(lat.var.latitude);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro   = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso    = zmeso*1000;         % mol/m3   -> mmol/m3
    phydiat  = phydiat*1000;       % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;       % mol/m3   -> mmol/m3
    phyc     = phydiat+phymisc;
    zooc     = zmicro+zmeso;
    bacc     = bacc*1000;         % mol/m3   -> mmol/m3
    
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC
    o2       = o2.*1000;             % mol/m3->mmol/m3 
    
    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
      %zooc(zooc<1e-15)=nan;        
    zooc(zooc>10e15) = nan;
    %phyc(phyc<1e-15)=nan;        
    phyc(phyc>10e15) = nan;
    %zmicro(zmicro<1e-15)=nan;    
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-15)=nan;      
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-15)=nan;  
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-15)=nan;  
    phymisc(phymisc>10e15) = nan;
    %intpp(intpp<1e-5)=nan;      
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;    
    epc100(epc100>10e15) = nan; 
    %bacc(bacc<1e-5)=nan;       
    bacc(bacc>10e15) = nan;
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan;  
    o2(o2>10e15) = nan; 

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    thetao   = circshift(thetao,-(lonoff-1),4);
    phyc     = circshift(phyc,-(lonoff-1),4);
    zooc     = circshift(zooc,-(lonoff-1),4);
    zmicro   = circshift(zmicro,-(lonoff-1),4);
    zmeso    = circshift(zmeso,-(lonoff-1),4);
    phydiat  = circshift(phydiat,-(lonoff-1),4);
    phymisc  = circshift(phymisc,-(lonoff-1),4);
    bacc     = circshift(bacc,-(lonoff-1),4);
    o2       = circshift(o2,-(lonoff-1),4);
    
    
    % 3 dimensions
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    tos          = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick   = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname,  'depthlev','phyc','zooc', 'zmicro', 'zmeso', 'phydiat', 'phymisc', 'bacc', 'intpp', 'epc100', 'tos','thetao','o2', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
 % PRIMARY Source: Lovato et al. (2022);  lovatoCMIP6SimulationsCMCC2022
                    % 2P2Z with multi-prey response w/ 
                       % Phytoplankton: diatoms, nano-flagellates
                       % Zooplankton:   Micro,   Meso
                       % Detrius:       
                       % Bacteria: Heterotrophic bacteria
                    % C units for Z
                      % Fixed stoichiometry
                      % Hpyotoplantk ohave C, N, P and Fe
               
 % Parameters in:  5.2: Based on Vichi et al 2020; vichiBiogeochemicalFluxModel2020
 %                 5.1: Or Vichi et al 2015: vichiBiogeochemicalFluxModel2015 
 %                 Parameters from personal correspondence
                   
 % Equations in:   Based on Vichi et al 2020; vichiBiogeochemicalFluxModel2020
 
% Zoo. Grazing Eqautions

 
                    % Auxilliary Variables
                      % Temperature Limitation 
                        % Class dependent T-lim terms scale the max class
                        %       specific grazing rate at 0 C
                         q10 = 2.0;  
                         ft = q10.^((thetao-10)/10);
          
                    % Z1: "MicroZooplankton"  Grazing On
                    % Prey: Nanoflaggelates, Diatoms, Bacteria, MicroZ
                         
                        % Base Preferencea
                        dprf_misc     = 1;
                        dprf_diat     = 0.01;
                        dprf_bacc     = 0.2;
                        dprf_zmicro   = 0.2;
                    
                        % Apparent Preferences
                        Kpref = 1.665; % mmolC/m3 converted From 20 mgC/m3
                        prf_misc = dprf_misc .* (phymisc./(phymisc + Kpref)); 
                        prf_diat = dprf_diat .* (phydiat./(phydiat + Kpref)); 
                        prf_bacc = dprf_bacc  .* (bacc./(bacc + Kpref));
                        prf_zmicro  = dprf_zmicro  .* (zmicro ./(zmicro + Kpref)); 
                        
                        % Total apparent food
                        food    = prf_misc.* phymisc + ...
                                  prf_diat.* phydiat + ...
                                  prf_bacc.* bacc + ... 
                                  prf_zmicro .* zmicro;
  
                        % "NanoPhytoplankton"
                        %  Single-prey w/ switching  
                        %  g   =  3.0 (1/d); 
                        %  K   =  1.665 (mmolC/m3); Converted from 20 mg/m3
                         g_Zmicro_phymisc = (3 .* ft .* prf_misc .* phymisc) ...
                                           ./(1.665 + food);
                   
                        % "Diatoms"
                        %  Single-prey w/ switching  
                        %  g   =  3.0 (1/d); 
                        %  K   =  1.665 (mmolC/m3); Converted from mg/m3
                         g_Zmicro_phydiat = (3 .* ft .* prf_diat .* phydiat) ...
                                           ./(1.665 + food);
                                       
                                       
                       % "Microzoo"
                        %  Single-prey w/ switching  
                        %  g   =  3.0 (1/d); 
                        %  K   =  1.665 (mmolC/m3); Converted from mg/m3
                         g_Zmicro_Zmicro = (3 .* ft .* prf_zmicro .* zmicro) ...
                                           ./(1.665 + food);
   
                                       
                    % Z2: "Mesozooplankton"  Grazing On
                    % Prey: Diatoms, Microzooplanton, Mesozooplankton
                      

                        % Base (and apparent) Preferencea
                        prf_diat     = 1.0;
                        prf_zmicro   = 0.5;
                        prf_zmeso    = 1.0;
                    
                        % Food  
                        food   = prf_diat.* phydiat + ...
                                 prf_zmicro.* zmicro + ...
                                 prf_zmeso.* zmeso;
                                     
                        % "Diatoms"
                        %  Single-prey w/ switching  
                        %  g   =  2.0 (1/d); 
                        %  e   =  0.025 (m3/mgC/d) 
                        %  K   =  6.66 (mmolC/m3); (g/e)/12 Converted from 0.025 (m3/mgC/d) 
                        %                        and g; using 12.011 = mgC/mmolC
                           g_Zmeso_phydiat = (2.0 .* ft .* prf_diat .* phydiat) ...
                                           ./(6.66 + food);
  
                                       
                        % "MicroZoo"
                        %  Single-prey w/ switching  
                        %  g   =  2.0 (1/d); 
                        %  e   =  0.025 (m3/mgC/d) 
                        %  K   =  6.66 (mmolC/m3); (g/e)/12 Converted from 0.025 (m3/mgC/d) 
                        %                        and g; using 12.011 = mgC/mmolC
                           g_Zmeso_Zmicro = (2.0 .* ft .* prf_zmicro .* zmicro) ...
                                           ./(6.66 + food);
                                       
                        % "MesoZoo"
                        %  Single-prey w/ switching  
                        %  g   =  2.0 (1/d); 
                        %  e   =  0.025 (m3/mgC/d) 
                        %  K   =  6.66 (mmolC/m3); (g/e)/12 Converted from 0.025 (m3/mgC/d) 
                        %                        and g; using 12.011 = mgC/mmolC
                           g_Zmeso_Zmeso = (2.0 .* ft .* prf_zmeso .* zmeso) ...
                                           ./(6.66 + food);
  
         
    % Total Grazing 
        GSP = g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro + ...
                                       +   g_Zmeso_phydiat.*zmeso;  

    % Total Trophic transfer Grazing 
        g_micro = .5; g_meso = .6; % Assimilation efficienies
        % Note assimilate efficeinc if Growth/ingestion; and some fraciton
        % of that immediatly excreted. But this is inclued in assimiliation
        % efficiency and only det3erines what pool it is routed to
        % B_micro = 0; B_meso = 0; % Excreted fraction of uptake (faces production)
      
        GSP_ht  =    g_micro.*(g_Zmicro_phymisc.*zmicro +  g_Zmicro_phydiat.*zmicro) + ...
                                       +   g_meso.*(g_Zmeso_phydiat.*zmeso);      
        
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
        
        
   % Bulk mortality 
      b_micro = 0.01; b_meso = 0.02; % basal specific Respiration rate
     
      fzo          =   (o2./(30+o2));                           % asssumed type 2 based on units of p_chro which ar mmolO2/m3 
      z_loss       =   (.01+.25.*(1-fzo)).*ft.*zmicro + b_micro.*ft.*zmicro  ...
                    +  (1-g_meso).*g_Zmeso_Zmicro.*zmeso + (1-g_micro).*g_Zmicro_Zmicro.*zmicro ...
                    +  0.01.*(1-fzo).*ft.*zmeso + .01.*zmeso.^2 + b_meso.*ft.*zmeso...
                    +  (1-g_meso).*g_Zmeso_Zmeso.*zmeso; 
                 
                  % Note: dont include canabalism on own class.
                  % Not:L Assume same Q10 for Z
                  % Note doesnt include explicit respiration
                  % Note doesnt include excreted from of uptake (Bz) based
                  % on elemental ratios
                  
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP_ht','GSP','z_loss','z_loss_spec', '-append')
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    phyC_inv      = squeeze(sum(phyc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)   .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)    .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv      = squeeze(sum(zooc(:,depthlev<dth,:,:)     .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    bacc_inv     = squeeze(sum(bacc(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    
    % Specific Rates (1/d)
    div_inv       = npp_inv./phyC_inv;
    graze_inv     = gsp_inv./zooC_inv; 
  
    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate',''};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)',''};
    
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',     'graze_inv', 'bacc_inv', 'gsp_inv','gsp_ht_inv', ...
                     'Lon',           'Lat',           'lonoff', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end





%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GFDL-ESM4.1 - COBALTv2 (USA) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Run:  Omon_GFDL-ESM4_historical_r1i1p1f1_gr
% Variables Used: phydiat | phymisc | phydiaz | phypico | zoomicro |
%                zoomacro | intpp | epc100 | baccos | 
%                 tos | areacello | volcello
% Caveats: DO diagnostically compute grazing rates and GSP. But also look at 
%          'graz' which is GSP (zoo grazing on phy). 
%          assum zmacro is 7% of meso (macro _meso) based on maredat median
%          obs
%          Bacteria is assume to covary with the the P profile relatvie to
%          surfce value
%          Total Z is zmeso (large + meso) + z micro 
%          Total P is diat + diaz + pico + misc;  
% Note:    zmeso is large AND medium Z
%          Not sure what phymisc is?
%          Graze (which is GSP) is saved explicitly
%          Only surface bacc (baccos) was saved
%          GSP is only Zooplankton consumption of primary producers (not
%          other Z, bac,  or detritus)

fname = 'GFDL-ESM4';

for Compute = 1
fprintf('\n Starting GFDL-ESM4.1 - COBALTv2 (USA)  \n')
cd /Users/roh029/Desktop/Data.nosync/CMIP6
 load([fname,'_Clims.mat'])  % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/Graze_Clim
cd /Users/roh029/Desktop/Data.nosync/CMIP6/coords
load(fname)        % Created in /glade/u/home/trohr/AAPP/Data/Analysis/CMIP6/CMIP6_Lat_Lon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Model Grid Dimensions')

    cd /Users/roh029/Desktop/Data.nosync/CMIP6/ModelGrids

     % Area
    areacello = ncread('areacello_Ofx_GFDL-ESM4_historical_r1i1p1f1_gr.nc','areacello');
    % loat areacello(nlon, nlat) | volcello:_FillValue = 1.e+20f | volcello:units = "m2" ;
    Cell_Area = permute(areacello(:,:),[2,1]); % Unit: m^2, Dim: (lat, lon)

    % Thickness
    thkcello = ncread('thkcello_Ofx_GFDL-ESM4_historical_r1i1p1f1_gr.nc','thkcello');
    % float volcello(i (lon), j (lat), lev) | volcello:_FillValue = 1.e+20f | volcello:units = "m3" ;
    Cell_Thick = permute(thkcello(:,:,1:35),[3,2,1]); % Unit: m, Dim: (level, lat, lon)
   
    % Volume
    Cell_Volume = Cell_Thick.*permute(repmat(Cell_Area,[1,1,35]),[3,1,2]); % Unit: m^3, Dim: (level, lat, lon)
    
         
   % Depth Level
    depthlev = ncread('thkcello_Ofx_GFDL-ESM4_historical_r1i1p1f1_gr.nc','lev'); 
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Regridding and Converting Units')
    
    %%%%%%%%  Regrid Cordinates to all start on 1 Degree East
    lon = double(lon.var.lon); lat= double(lat.var.lat);
    [Lon, Lat] = meshgrid(lon,lat);
    lon_dum = Lon(1,:);
    [a,lonoff]=min(abs(lon_dum-0)); % Find 1 degree east
    Lon = circshift(Lon,-(lonoff-1),2);
    Lat = circshift(Lat,-(lonoff-1),2);
        % check that longitude starts at 1E
        fprintf(['\n      Shifted Lon is: ' ,num2str(Lon(1,1))])
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'Lon', 'Lat', 'lonoff')

    %%%%%%%%% Check and Convert Units %%%%%%%%%%%
    zmicro   = zmicro*1000;        % mol/m3   -> mmol/m3
    zmeso    = zmeso*1000;         % mol/m3   -> mmol/m3
    phydiat  = phydiat*1000;       % mol/m3   -> mmol/m3
    phymisc  = phymisc*1000;       % mol/m3   -> mmol/m3
    phypico  = phypico*1000;       % mol/m3   -> mmol/m3
    phydiaz  = phydiaz*1000;       % mol/m3   -> mmol/m3
    phyc     = phydiat + phymisc + ...
               phydiaz + phypico;
    zooc     = zmicro  + zmeso;
    baccos   = baccos*1000;          % mol/m3   -> mmol/m3
    
    intpp    = intpp*1000*60*60*24;  % mol/m2/s -> mmol/m3/d
    epc100   = epc100*1000*60*60*24; % mol/m2/s -> mmol/m3/d
    graz     = graz*1000*60*60*24;   % mol/m3/s -> mmol/m3/d
    tos      = tos;                  % degC     -> degC

    %%%%%%%%%% Check and Remove Nan/Fill/Land/Neg values %%%%%%%
    % FillValue = 1.e+20f; Land values appear to vary though
    %phyc(phyc<1e-5)=nan;           
    phyc(phyc>10e15) = nan; 
    %zooc(zooc<1e-5)=nan;           
    zooc(zooc>10e15) = nan; 
    %zmicro(zmicro<1e-5)=nan;       
    zmicro(zmicro>10e15) = nan; 
    %zmeso(zmeso<1e-5)=nan;         
    zmeso(zmeso>10e15) = nan; 
    %phydiat(phydiat<1e-5)=nan;     
    phydiat(phydiat>10e15) = nan; 
    %phymisc(phymisc<1e-5)=nan;     
    phymisc(phymisc>10e15) = nan;
    %phydiaz(phydiaz<1e-5)=nan;     
    phydiaz(phydiaz>10e15) = nan; 
    %phypico(phypico<1e-5)=nan;     
    phypico(phypico>10e15) = nan;
    %baccos(baccos<1e-5)=nan;       
    baccos(baccos>10e15) = nan;
    %graz(graz<1e-5)=nan;           
    graz(graz>10e15) = nan;
    %intpp(intpp<1e-5)=nan;         
    intpp(intpp>10e15) = nan; 
    %epc100(epc100<1e-5)=nan;       
    epc100(epc100>10e15) = nan; 
    tos(tos>10e15) = nan; 
    thetao(thetao>10e15) = nan; 

    %%%%%%%%%% Phase shift with lon %%%%%%%%%%%%%
    % 4 dimesnsions
    thetao   = circshift(thetao,-(lonoff-1),4);
    phyc     = circshift(phyc,-(lonoff-1),4);
    zooc     = circshift(zooc,-(lonoff-1),4);
    zmicro   = circshift(zmicro,-(lonoff-1),4);
    zmeso    = circshift(zmeso,-(lonoff-1),4);
    phydiat  = circshift(phydiat,-(lonoff-1),4);
    phymisc  = circshift(phymisc,-(lonoff-1),4);
    phydiaz  = circshift(phydiaz,-(lonoff-1),4);
    phypico  = circshift(phypico,-(lonoff-1),4);
    graz     = circshift(graz,-(lonoff-1),4);
    
    % 3 dimensions
    baccos       = circshift(baccos,-(lonoff-1),3);
    intpp        = circshift(intpp,-(lonoff-1),3); 
    epc100       = circshift(epc100,-(lonoff-1),3);
    tos          = circshift(tos,-(lonoff-1),3);
    Cell_Volume  = circshift(Cell_Volume,-(lonoff-1),3);
    Cell_Thick   = circshift(Cell_Thick,-(lonoff-1),3);
    
    % 2 Dimensions
    Cell_Area  = circshift(Cell_Area,-(lonoff-1),2);
    
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims
    save(fname, 'depthlev', 'phyc','zooc', 'zmicro', 'zmeso', 'phydiat', 'phymisc', 'phydiaz', 'phypico', 'graz', 'baccos', 'intpp', 'epc100', 'tos','thetao', '-append')
    save(fname, 'Cell_Thick', 'Cell_Area', 'Cell_Volume', '-append')
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Diagnostic Grazing Rates')
 % PRIMARY Source: Stock et al. (2020);  stockOceanBiogeochemistryGFDL2020
                    % 3P3Z with multi-prey response w/ mild density dependent swithcing
                       % Phytoplankton: Small, Large, Diaz 
                       % Zooplankton:   Micro, Meso, Macro
                       % Detrius:       POC
                       % Bacteria
                    % N units
                      % 106:16 C:N for All Phytoplankton
             
                    
 % Parameters in:  stockOceanBiogeochemistryGFDL2020  | Everything but Prefences
 %                 stockBottomTopForcing2008          |  Preferences
                   
 % Equations in:   stockOceanBiogeochemistryGFDL2020
                    
 
 % NOTES: Same Parameter used in COBALT (stock et al, 2014)
  
 % Determine missing variables

 

 for d=1:35; psd(:,d,:,:,:)=squeeze(phyc(:,d,:,:))./squeeze(phyc(:,1,:,:));end
 bacc=permute(repmat(baccos,[1,1,1,35]),[1,4,2,3]).*psd; % Assume homogenous profile of bacteria 
 
 % Meso macro Z partiitoning - PEr cOmms woth Charlie Stock
 an_npp = squeeze(nanmean(intpp)); %  Off hand, I would say that large is likely a bit more than medium in highly productive regions, comparable or a bit lower in mesotrophic regions, and much lower in the heart of the oligotrophic gyre.  
 biomes = prctile(an_npp(:), [33,66]);
 b_map = an_npp.*0+3;
 b_map(an_npp<biomes(2))=2;
 b_map(an_npp<biomes(1))=1;
  
 zmes = zmeso; zmac = zmeso;            % Assume 93% of "mesozooplankton" are medium vs large Z based on 
 % Highly productive region             % large is likely a bit more than medium in highly productive regions
 zmes(b_map==3) = zmes(b_map==3).*0.4; 
 zmac(b_map==3) = zmac(b_map==3).*0.6;
 
 % mestropic regions                     % large comparable or a bit lower in mesotrophic regions
 zmes(b_map==2) = zmes(b_map==2).*0.6; 
 zmac(b_map==2) = zmac(b_map==2).*0.4;
 
% oligotropic regions                     % large much lower in the heart of the oligotrophic gyre
 zmes(b_map==1) = zmes(b_map==1).*0.9; 
 zmac(b_map==1) = zmac(b_map==1).*0.1;
 
  % Zoo. Grazing Eqautions .  --- Just for reference as using explicit
  % 'graz'

                    % Auxilliary Variables
                      % Temperature Limitation 
                        % Class dependent T-lim terms scale the max class
                        %       specific grazing rate at 0 C
                         k_zoo = 0.063; 
                         ft = exp(k_zoo.*thetao);
                      % Phyotplankton Threshold 
                        p_thresh = 0.001;
                       
                        
                        
                    % Z1: "MicroZooplankton"  Grazing On
                    % Prey: Small phyto, Bacteria 
                        
                        % Base Preferencea
                        prf_sp     = 1;
                        prf_B      = 0.25;
                    
                        % Apparent Prerence 
                        Prf_sp     = (prf_sp.*phypico) ./ sqrt( (prf_sp.*phypico).^2 + (prf_B.*bacc).^2 ) ;
                        Prf_B      = (prf_B.*bacc)   ./ sqrt( (prf_sp.*phypico).^2 + (prf_B.*bacc).^2 );
                        
                        % Total apparent food
                        food = Prf_sp.*phypico  + ...
                                  Prf_B.*bacc ;
                        
                        % "Small Phytoplankton"
                        %  Single-prey w/ switching  
                          g   =  1.28; %(1/d); 
                          K   =  8.28; %(mmolC/m3);
                          g_COBALT_Z1sp = (g .* ft .* Prf_sp .* phypico) ...
                                           ./(K + food);
                   
              % Z2: "Medium Zooplankton"  Grazing On
                    % Prey: Large phyto, Diazatrophs, Small Zooplankton
                   
                        
                        % Base Preferencea
                        prf_lp     = 1;
                        prf_diaz   = 1;
                        prf_sz     = 1;
                    
                        % Apparent Prerence 
                        Prf_lp     = (prf_lp.*(phydiat+phymisc))       ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_sz.*zmicro).^2) ;
                        Prf_diaz   = (prf_diaz.*phydiaz)   ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_sz.*zmicro).^2) ;
                        Prf_sz     = (prf_sz.*zmicro)   ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_sz.*zmicro).^2) ;
                        
                        % Total apparent food
                        food = Prf_lp.*(phydiat+phymisc)  + Prf_diaz.*phydiaz + Prf_sz.*zmicro ;
                        
                        % "Large Phytoplankton"
                        %  Single-prey w/ switching  
                          g   =  0.57; %(1/d); 
                          K   =  8.28; %(mmolC/m3); 
                         g_COBALT_Z2lp = (g .* ft .* Prf_lp .* (phydiat+phymisc)) ...
                                           ./(K + food);
                                       
                        % "Diazatrophs"
                        %  Single-prey w/ switching  
                          g   =  0.57; %(1/d); 
                          K   =  8.28; %(mmolC/m3);
                         g_COBALT_Z2diaz = (g .* ft .* Prf_diaz .* phydiaz) ...
                                           ./(K + food);
                   
                       % "small Z"
                        %  Single-prey w/ switching  
                          g   =  0.57; %(1/d); 
                          K   =  8.28; %(mmolC/m3);
                         g_COBALT_Z2_Z1 = (g .* ft .* Prf_sz .* zmicro) ...
                                           ./(K + food);
                                       
                  % Z3: "Large Zooplankton"  Grazing On
                    % Prey: Large phyto, Diazatrophs, Medium Zooplankton
                    
                        
                        % Base Preferencea
                        prf_lp     = 1;
                        prf_diaz   = 1;
                        prf_mz     = 1;
                    
                        % Apparent Prerence 
                        Prf_lp     = (prf_lp.*(phydiat+phymisc))       ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_mz.*zmes).^2) ;
                        Prf_diaz   = (prf_diaz.*phydiaz)   ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_mz.*zmes).^2) ;
                        Prf_mz     = (prf_mz.*zmes)   ...
                                     ./ sqrt( (prf_lp.*(phydiat+phymisc)).^2 + (prf_diaz.*phydiaz).^2 + (prf_mz.*zmes).^2) ;
                        
                        % Total apparent food
                        food = Prf_lp.*(phydiat+phymisc)  + Prf_diaz.*phydiaz + Prf_mz.*zmes ;
                        
                        % "Large Phytoplankton"
                        %  Single-prey w/ switching  
                          g   =  0.23; %(1/d); 
                          K   =  8.28; %(mmolC/m3); Converted from 20 umol C/L
                         g_COBALT_Z3lp = (g .* ft .* Prf_lp .* (phydiat+phymisc)) ...
                                           ./(K + food);
                                       
                 
                        % "Diazatrophs"
                        %  Single-prey w/ switching  
                          g   =  0.23; %(1/d); 
                          K   =  8.28; %(mmolC/m3); Converted from 20 umol C/L
                         g_COBALT_Z3diaz = (g .* ft .* Prf_diaz .* phydiaz) ...
                                           ./(K + food);
                                       
                        % "Medium Zooplankton"
                        %  Single-prey w/ switching  
                          g   =  0.23; %(1/d); 
                          K   =  8.28; %(mmolC/m3); Converted from 20 umol C/L
                          g_COBALT_Z3_Z2= (g .* ft .* Prf_mz .* zmes) ...
                                           ./(K + food);
                                       
                                       
                  % Z4: "Higher Predators"  Grazing On
                    % Prey: Medium Zooplankton, Large Zooplankton
                    
                        
                        % Base Preferencea
                        prf_lz     = 1;
                        prf_mz     = 1;
                    
                        % Apparent Prerence 
                        Prf_lz   = (prf_lz.*zmac)   ...
                                     ./ sqrt( (prf_lp.*(zmes)).^2 + (prf_mz.*zmac).^2) ;
                        Prf_mz     = (prf_mz.*zmes)   ...
                                      ./ sqrt( (prf_lp.*(zmes)).^2 + (prf_mz.*zmac).^2);
                        
                        % Total apparent food
                        food = Prf_lz.*(zmac)  + Prf_diaz.*zmes  ;
                        
                        % "Large zooplankton"
                        %  Single-prey w/ switching  
                          g   =  0.09; %(1/d); 
                          K   =  8.28; %(mmolC/m3); Converted from 20 umol C/L
                         g_COBALT_HPlz = (g .* ft .* Prf_lz .* zmac) ...
                                           ./(K + food);
                                       
                 
                        % "Medium Zooplankton"
                        %  Single-prey w/ switching  
                          g   =  0.09; %(1/d); 
                          K   =  8.28; %(mmolC/m3); Converted from 20 umol C/L
                         g_COBALT_HPmz = (g .* ft .* Prf_mz .* zmes) ...
                                           ./(K + food);
                                       
                                       
                                       
                   
               
    % Total Grazing 
        GSP = g_COBALT_Z1sp.*zmicro + ...
              g_COBALT_Z2lp.*zmes   + g_COBALT_Z2diaz.*zmes + ...
              g_COBALT_Z3lp.*zmac   + g_COBALT_Z3diaz.*zmac;                         ;  

    % Total Trophic Transfer
    % note .7 referes to 70% that is ingested and not routed throgh slopy
    % feeding to detritus.  .4 is ussed for assimilation effifency which
    % also considers metabolic processes (ths not gross).  
     g = .4;
     GSP_ht = g.*GSP;
     
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec = GSP./(zmicro+zmeso);
        
                                       
    % Total Grazing 
        GSP_exp = graz;
        
    % Bulk Zooplankton Specific Grazing Rates  
        graze_spec_exp = GSP_exp./(zmicro+zmeso);
    
    
     % Bulk mortality 
     % HP biomassadjusts to M+L Z pool
      Bz1 =.018;  Bz2 =.008;   Bz3 =.0032; g = .4;
      tfr = 2.^(thetao./10); % T limitaiton from stockControlsRatioMesozooplankton2010
      z_loss       =  + Bz1.*tfr.*zmicro + Bz2.*tfr.*zmes + Bz3.*tfr.*zmac ...
                      + g_COBALT_HPlz .* (zmeso)  ...
                      + g_COBALT_HPmz .* (zmeso) ...   
                      + (1-g).* g_COBALT_Z2_Z1.*zmes + (1-g).*g_COBALT_Z3_Z2.*zmac;  % innefficeinciers
      z_loss_spec  =  z_loss./(zmicro+zmeso);
      
        
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/Processed_Clims   
    save(fname, 'graze_spec', 'GSP','GSP_ht', 'graze_spec_exp', 'GSP_exp','z_loss','z_loss_spec',  '-append')
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Depth Integrated Tracers')
clearvars CoreVar_Inv

    % Tracers (mmol/m2)  --  Multiply concentration by gridcell height then multiply
    phydiat_inv   = squeeze(sum(phydiat(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phymisc_inv   = squeeze(sum(phymisc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    phydiaz_inv   = squeeze(sum(phydiaz(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phypico_inv   = squeeze(sum(phypico(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    phyC_inv      = squeeze(sum(phyc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmicro_inv    = squeeze(sum(zmicro(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zmeso_inv     = squeeze(sum(zmeso(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    zooC_inv      = squeeze(sum(zooc(:,depthlev<dth,:,:)  .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan'));
    
    % Bulk Rates (mmol/m2/day)
    npp_inv     = intpp;
    exp_inv     = epc100;
    gsp_inv     = squeeze(sum(GSP(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_ht_inv     = squeeze(sum(GSP_ht(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    gsp_inv_exp     = squeeze(sum(GSP_exp(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    
    % Specific Rates (1/d)
    div_inv         = npp_inv./phyC_inv;
    graze_inv       = gsp_inv./zooC_inv; 
    graze_inv_exp   = gsp_inv_exp./zooC_inv; 

    % Combine Standard Variables for Easy Computations
    CoreVar_Inv(1,:,:,:) =  phyC_inv;
    CoreVar_Inv(2,:,:,:) =  zooC_inv;
    CoreVar_Inv(3,:,:,:) =  npp_inv;
    CoreVar_Inv(4,:,:,:) =  gsp_inv;
    CoreVar_Inv(5,:,:,:) =  exp_inv;
    CoreVar_Inv(8,:,:,:) =  gsp_ht_inv;
    % Rates
    CoreVar_Inv(6,:,:,:) =  div_inv;
    CoreVar_Inv(7,:,:,:) =  graze_inv;
    
    % Add Loss terms
    Z_Loss_inv     = squeeze(sum(z_loss(:,depthlev<dth,:,:)      .* permute(repmat(Cell_Thick(depthlev<dth,:,:),[1,1,1,12]),[4,1,2,3]), 2,'omitnan')); 
    Z_loss_inv     = Z_Loss_inv./zooC_inv;
    CoreVar_Inv(9,:,:,:) =  Z_Loss_inv;
    CoreVar_Inv(10,:,:,:) = Z_loss_inv;
    
    
    CoreVar_Names = {'Phytoplankton Biomass',    'Zooplankton Biomass',...
                         'NPP',                 'GSP',            'EP',...
                         'Spec. Division Rate', 'Spec. Grazing Rate'};
    CoreVar_Inv_Units = {'(mmol $m^{-2}$)',          '(mmol $m^{-2}$)',...
                         '(mmol $m^{-2}$ $d^{-1}$)', '(mmol $m^{-2}$ $d^{-1}$)',  '(mmol $m^{-2}$ $d^{-1}$)'...
                         '($d^{-1}$)',               '($d^{-1}$)'};
    
    % Save
    cd  /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname,      'phydiat_inv',   'phymisc_inv',   'phydiaz_inv',   'phypico_inv', 'phyC_inv',   'zmicro_inv', 'zmeso_inv', 'zooC_inv', ...
                     'exp_inv',       'npp_inv',       'div_inv',     'graze_inv', 'gsp_inv', 'graze_inv_exp', 'gsp_inv_exp',...
                     'Lon',           'Lat',           'lonoff','gsp_ht_inv', ...
                     'CoreVar_Inv',   'CoreVar_Names', 'CoreVar_Inv_Units');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n    Computing Spatial Integrals')

    mmolC_to_PG = 12.*1e-18;
    G_lats = find(Lat(:,1)>latmin & Lat(:,1)<latmax);
    N_lats = find(Lat(:,1)>0 & Lat(:,1)<latmax);
    S_lats = find(Lat(:,1)<0 & Lat(:,1)>latmin);
    
    % Global 
        % Integrals
        CoreVar_Int(1,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,G_lats,:) .*...
                                         permute(repmat(Cell_Area(G_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(1,6,:)   = CoreVar_Int(1,3,:)./CoreVar_Int(1,1,:);
        CoreVar_Int(1,7,:)   = CoreVar_Int(1,4,:)./CoreVar_Int(1,2,:);
        CoreVar_Int(1,10,:)   = CoreVar_Int(1,9,:)./CoreVar_Int(1,2,:);

    % Nothern Hemisphre
        % Integrals
        CoreVar_Int(2,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,N_lats,:) .*...
                                         permute(repmat(Cell_Area(N_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(2,6,:)   = CoreVar_Int(2,3,:)./CoreVar_Int(2,1,:);
        CoreVar_Int(2,7,:)   = CoreVar_Int(2,4,:)./CoreVar_Int(2,2,:);
        CoreVar_Int(2,10,:)   = CoreVar_Int(2,9,:)./CoreVar_Int(2,2,:);

    % Southern Hemisphere
        % Integrals
        CoreVar_Int(3,[1:5,8,9],:) =  nansum(nansum(CoreVar_Inv([1:5,8,9],:,S_lats,:) .*...
                                         permute(repmat(Cell_Area(S_lats,:),[1,1,7,12]),[3,4,1,2]) ...
                                  ,3),4) .* mmolC_to_PG;
        % Averages                      
        CoreVar_Int(3,6,:)   = CoreVar_Int(3,3,:)./CoreVar_Int(3,1,:);
        CoreVar_Int(3,7,:)   = CoreVar_Int(3,4,:)./CoreVar_Int(3,2,:);
        CoreVar_Int(3,10,:)   = CoreVar_Int(3,9,:)./CoreVar_Int(3,2,:);
        
    % Convert Productivity to per Year (from per day)
      CoreVar_Int(:,[3:5,8,9],:) =  CoreVar_Int(:,[3:5,8,9],:).*365;
      
    % Units
    CoreVar_Int_Units = {'PgC',          'PgC',...
                         'PgC/yr',       'PgC/yr',    'PgC/yr'...
                         '($d^{-1}$)',   '($d^{-1}$)', 'PgC/yr','PgC/yr','($d^{-1}$)'};
  
    % Save
    cd /Users/roh029/Desktop/Data.nosync/CMIP6/DepthInv
    save(fname, 'CoreVar_Int','CoreVar_Int_Units','-append');

end

end