% PLOTALL
%%%%%%%%%%%%%%%
% plots horizontal maps, zonal transects, streamfunction, overturning
% etc.
% cleaned, updated: 9 Aug 01 swd
% change from plotresults: 8 May 01 swd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up for several grids...(via tags: bathynum and ocmipornot)
% USE AT YOUR OWN RISK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
 clear all
 path(path,'/data25/stephd/Global/Functions');
%
%%%%%%%%%% CHOICES %%%%%%%%%%%%%%%%
%%%%%%%%%%
  titstr=['0204.03'];
  yeard=365;
  run1str=['/data25/stephd/c46c3/newerice/0204.03'];
  timestr=['0003139000'];   dt=43200;
  bathynum=10;    % bathymetric set up
  ocmipornot=0;  %1=ocmip setup
%%%%%%%%%
  kn=2;       % level to plot
  knsp=3;    % level to plot for subplt
  choice=2;  % 1=snapshot; 2=mean
%%%%%%%%%
  if (ocmipornot==1), bathynum=20; end
  atlanornot=0;  %1=turn fields to atlantic centric
%
  if (bathynum==3), direc=['/data25/stephd/Newdata/input_thirdtry']; end
  if (bathynum==9), 
       direc=['/data25/stephd/Newdata/input_thirdtry_nomed']; end
  if (bathynum==10),
       direc=['/data25/stephd/Nesdata/input_tenthtry_ts']; end
  if (bathynum==11),
       direc=['/data25/stephd/Nesdata/input_tenthtry_noarctic']; end
  if (bathynum==12), 
       direc=['/data25/stephd/Nesdata/input_tenthtry_ts_ideal_sill2']; end
  if (bathynum==13), 
       direc=['/data25/stephd/Nesdata/input_tenthtry_ts_ideal_sill2_var2']; end
  if (bathynum==14), 
       direc=['/data25/stephd/Nesdata/input_tenthtry_ts_newideal_new']; end
  if (ocmipornot==1), direc=['/data25/stephd/c44/ocmip/data_ocmip']; end
  if (bathynum==12 | bathynum==13 | bathynum==14),
     dzfile=['dz_ideal']; 
  else
     dzfile=['dz_vero'];
  end % if
%%%%%%%%%%%%%%%%%%%   END CHOICES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% some parameters
  dtday=dt/(60*60*24); degm=112*1e3; degrad=pi/180; kgm=1;
%%%%%%%%%%%%%%%%%%%%%%%% GRID %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% horizontal grid
 if ((bathynum>0 & bathynum<13) | bathynum==14), 
  dx=4; dy=4;
  latg=[-90,-84:dy:84,90]; longg=[0:dx:360];
  latc=[-87,-82:dy:82,87]; longc=2:dx:358;
 end 
 if (bathynum==13),
   dx=4; dy=4;
   latg=[-90,-84:dy:84,90]; latc=[-87,-82:dy:82,87]; 
   i1=24; i2=136; i3=300;
   longg=[0:dx:i1 i1+3.6  i1+6.8  i1+9.6 ...
                 i1+12 i1+14 i1+16 ...
                 i1+18 i1+20 i1+22.4 ...
                 i1+25.2 i1+28.4 i1+32 ...
          i1+36:dx:i2 i2+3.6  i2+6.8  i2+9.6 ...
                 i2+12 i2+14 i2+16 ...
                 i2+18 i2+20 i2+22.4 ...
                 i2+25.2 i2+28.4 i2+32 ...
          i2+36:dx:i3 i3+3.6  i3+6.8  i3+9.6 ...
                 i3+12 i3+14 i3+16 ...
                 i3+18 i3+20 i3+22.4 ...
                 i3+25.2 i3+28.4 i3+32 ...
          i3+36:4:360];
    for i=1:length(longg)-1
      longc(i)=longg(i)+( longg(i+1)-longg(i) ) /2;
    end % for i
 end % if
 if (bathynum==20),
  dx=2.8125; dy=2.8125;
  latc=[-88.59375:2.8125:88.59375]; longc=1.40625:2.8125:358.59375; 
  latg=[-90:2.8125:90]; longg=0:2.8125:360;
 end % if
 latu=latc;     longu=longg; latv=latg;     longv=longc;
 nx=length(longc); ny=length(latc); nx1=nx-1; ny1=ny-1;
%
% find length between grid points and center points
 for i=1:nx, for j=1:ny
    dyg(i,j)=(latg(j+1)-latg(j))*degm;
    if (j<ny), dyc(i,j)=(latc(j+1)-latc(j))*degm; end
 end, end % for i,j
 dyc(:,ny)=dyg(:,ny)/2+dyg(:,1)/2;
 for i=1:nx, for j=1:ny
     dxg(i,j)=(longg(i+1)-longg(i))*cos(latg(j)*degrad)*degm;
     if (i<nx) dxc(i,j)=(longc(i+1)-longc(i))*cos(latc(j)*degrad)*degm; end
 end, end % for i,j
 dxc(nx,:)=dxg(nx1,:)/2+dxg(1,:)/2;
% dxcg at center of grid
 for i=1:nx, for j=1:ny
     dxcg(i,j)=(longg(i+1)-longg(i))*cos(latc(j)*degrad)*degm;
 end, end % for i j
% vertical resolution
  load ([dzfile,'.dat']); eval(['dz=',dzfile,';']); dz=dz';
  if (bathynum==14), clear tmp, tmp=dz(1:length(dz)-1); clear dz, dz=tmp; end
  depth=cumsum(dz); nz=length(dz); depthm(1)=0+dz(1)/2;
  for k=2:nz, depthm(k)=depthm(k-1)+dz(k-1)/2+dz(k)/2; end
%%%%%%
% face areas
  for k=1:nz,
    area_xyc(:,:,k)=dxcg(:,:).*dyg(:,:);
    area_xzg(:,:,k)=dxg(:,:)*dz(k);
    area_yzg(:,:,k)=dyg(:,:)*dz(k);
    volc(:,:,k)=area_xyc(:,:,k)*dz(k);
  end % for k
% now that areas are found, convert all nx+1,ny+1 arrays to nx,ny
  clear tmp, tmp=longg(1:nx); clear longg, longg=tmp;
  clear tmp, tmp=longu(1:nx); clear longu, longu=tmp;
  clear tmp, tmp=dxg(1:nx,:); clear dxg, dxg=tmp;
  clear tmp, tmp=latg(1:ny); clear latg, latg=tmp;
  clear tmp, tmp=latv(1:ny); clear latv, latv=tmp;
  clear tmp, tmp=dyg(:,1:ny); clear dyg, dyg=tmp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
  if (atlanornot==1)
    olongg=longg; olongc=longc; olongu=longu; olongv=longv;
    longg=turnlong(olongg,longg); longc=turnlong(olongc,longc); 
    longu=turnlong(olongu,longu); longv=turnlong(olongv,longv);
    dyg=turn2(olongg,dyg); dyc=turn2(olongc,dyc);
    dxg=turn2(olongg,dxg); dxc=turn2(olongc,dxc);
    area_xyc=turn3(olongg,area_xyc); area_xzg=turn3(olongg,area_xzg);
    area_yzg=turn3(olongg,area_yzg); volc=turn3(olongg,volc);
  end % if atlanornot 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read bathymetry and masks
  if (ocmipornot==1),
   fid=fopen([direc,'/bathymetry.bin'],'r','b');
  else
   fid=fopen([direc,'/bathy.bin'],'r','b');
  end % if
  bathy=fread(fid,'float32'); bathy=reshape(bathy,nx,ny);
% masks
  fid=fopen([direc,'/hFacC.bin'],'r','b');
  if (bathynum==11),  fid=fopen([direc,'/pmask.bin'],'r','b'); end
  pmask=fread(fid,'float32'); pmask=reshape(pmask,nx,ny,nz);
% find u,v mask
  for k=1:nz, for i=1:nx, for j=2:ny,
      vmask(i,j,k)=min(pmask(i,j,k),pmask(i,j-1,k));
  end, end, end % for k,i,j
  for k=1:nz, for i=2:nx, for j=1:ny,
      umask(i,j,k)=min(pmask(i,j,k),pmask(i-1,j,k));
      if (i==2), umask(1,j,k)=min(pmask(1,j,k),pmask(nx,j,k)); end
  end, end, end % for k,i,j
% masks for each ocean (down to antarctica)
  fid=fopen([direc,'/tocean.bin'],'r','b');
  tocean=fread(fid,'float32'); tocean=reshape(tocean,nx,ny);
% masks for each ocean, including southern ocean
  fid=fopen([direc,'/cocean.bin'],'r','b');
  cocean=fread(fid,'float32'); cocean=reshape(cocean,nx,ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% original T and S fields (ie Levitus)
  if (ocmipornot==1), 
    fid=fopen([direc,'/lev_clim_temp.bin'],'r','b'); monnum=1;
    theta=fread(fid,'float32');
    theta=reshape(theta,nx,ny,nz); theta_clim=theta;
  else
    fid=fopen([direc,'/theta.bin'],'r','b'); monnum=12;
    theta=fread(fid,'float32'); 
    theta=reshape(theta,nx,ny,nz,monnum); theta_clim=mean(theta,4);
  end
% clear fi, fi=find(isnan(theta_clim)==1); theta_clim(fi)=900;
% clear fi, fi=find(theta_clim==0); theta_clim(fi)=900;
%
  if (ocmipornot==1), 
    fid=fopen([direc,'/lev_clim_salt.bin'],'r','b'); monnum=1;
    salt=fread(fid,'float32');
    salt=reshape(salt,nx,ny,nz); salt_clim=salt;
  else
    fid=fopen([direc,'/sal.bin'],'r','b'); monnum=12;
    salt=fread(fid,'float32');
    salt=reshape(salt,nx,ny,nz,monnum); salt_clim=mean(salt,4);
  end
%%%%%
  if (atlanornot==1),
     bathy=turn2(olongc,bathy); pmask=turn3(olongc,pmask);
     umask=turn3(olongu,umask); vmask=turn3(olongv,vmask);
     tocean=turn2(olongc,tocean); cocean=turn2(olongc,cocean);
     theta_clim=turn3(olongc,theta_clim);
     salt_clim=turn3(olongc,salt_clim);
  end % if atlanornot
%%%%%
% land 
  filand=find(pmask==0); filand1=find(pmask(:,:,1)==0);
  if (kn>1); depthkn1=depth(kn-1); else; depthkn1=0; end
% salt_clim(filand)=NaN; theta_clim(filand)=NaN;
% time
  time=str2num(timestr)*dt; timeyear=time/(yeard*24*60*60);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
%  SETUP PLOTS    %
%%%%%%%%%%%%%%%%%%%
 figure(1), hold off, subplot(1,1,1), pcolor(longc,latc,bathy');
 colorbar, title([titstr,': bathymetry']);
 figure(2), hold off, subplot(1,1,1), pcolor(longc,latc,tocean');
 title([titstr,': Ocean Masks']);
 figure(3), hold off, subplot(1,1,1), subplot(3,1,1)
 pcolor(longc,latc,theta_clim(:,:,kn)'); shading flat, colorbar,
 title([titstr,': Levitus Theta, ',num2str(depthm(kn)),'m'])
 figure(4), hold off, subplot(1,1,1),  subplot(3,1,1)
 pcolor(longc,latc,salt_clim(:,:,kn)'); shading flat, colorbar,
 title([titstr,': Levitus Salinity, ',num2str(depthm(kn)),'m'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize summary plots %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fgn=115;
 figure(fgn), hold off, subplot(1,1,1);
 figure(fgn+1), hold off, subplot(1,1,1);
 figure(fgn+2), hold off, subplot(1,1,1);
 figure(fgn+3), hold off, subplot(1,1,1);
% global
  ax1=1; ax2=360; ax3=-90; ax4=90; inton= 20; vmult=150;
  if (atlanornot==1), ax1=-180; ax2=180; end
% north atlantic for subplt (note could also be another level);
  ax1_sp=280; ax2_sp=360; ax3_sp=-10; ax4_sp=90; inton_sp=5; vmult_sp=200;
  if (atlanornot==1), ax1_sp=-100; ax2_sp=20; end
%%%%%%%%%%%%
%
 ['finished set up, now starting reading fields'],
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     MAIN FIELDS                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (choice==1), maxfld=6; else, maxfld=7; end
if (choice>0),
 figure(10), hold off, subplot(1,1,1), 
 for ifld=1:maxfld,
  if (ifld==1), if (choice==1), str=['T']; else, str=['Ttave']; end
                lat=latc; long=longc; cax1=-2; cax2=28; 
                cint1=-20; cint2=50; cint=5;  end
  if (ifld==2), if (choice==1), str=['S']; else, str=['Stave']; end
                lat=latc; long=longc;  cax1=28; cax2=40;
                cint1=20; cint2=50; cint=.5; end
  if (ifld==3), if (choice==1), str=['U']; else, str=['uVeltave']; end
                lat=latu; long=longu; end
  if (ifld==4), if (choice==1), str=['V']; else, str=['vVeltave']; end
                lat=latv; long=longv; end
  if (ifld==5), if (choice==1), str=['W']; else, str=['wVeltave']; end
                lat=latc; long=longc; cax1=-9e-6; cax2=6e-6; cint=3e-6; end
  if (ifld==6), if (choice==1), str=['Eta']; else, str=['ETAtave']; end
                lat=latc; long=longc; cax1=-2; cax2=1.2; 
                cint1=-5; cint2=5; cint=.25; end
  if (ifld==7 & choice>1 ), str=['Convtave']; lat=latc; long=longc; 
                                      cax1=0; cax2=1; 
                                      cint1=0; cint2=1; cint=.2; end
% read fields
   fid=fopen([run1str,'/',str,'.',timestr,'.data'],'r','b');
   clear field, field=fread(fid,'float32');
   if (ifld==6), 
        field=reshape(field,nx,ny);
        if (atlanornot==1), field=turn2(olongc,field); end
   else
        field=reshape(field,nx,ny,nz); 
        if (atlanornot==1), 
            if (ifld==1 |ifld==2| ifld==5 | ifld==6); olong=olongc; end
            if (ifld==3), olong=olongu; end, if (ifld==4), olong=olongv; end
            field=turn3(olong,field); end
   end  % ifld
%%%
   if (ifld==1), temp=field;  clear filand_temp,
                 filand_temp=find(temp==0); 
                 filand1_temp=find(temp(:,:,1)==0); 
                 temp_pmask=ones(size(temp)); temp_pmask(filand_temp)=0; end
   if (ifld==2), salt=field; end
   if (ifld==3), uvel=field; end
   if (ifld==4), vvel=field; end
   if (ifld==5), wvel=field; end
   if (ifld==6), eta=field; end
   if (ifld==7), conv=field; end
%%%%%%%%%%%%%%%%%%%%
% calculate max/min of each field
   clear tmp, tmp=field;
   if (ifld~=6), tmp(filand_temp)=NaN; else, tmp(filand1_temp)=NaN; end
   totmax(ifld)=max(max(max(tmp))); 
   totmin(ifld)=min(min(min(tmp))); 
% calculate mean
   if (ifld~=6), volume=volc; volume(filand_temp)=0;
   else, volume=volc(:,:,1); volume(filand1_temp)=0; end
   tmp=field.*volume; totmean(ifld)=sum(sum(sum(tmp)))/(sum(sum(sum(tmp))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate zonal transports
   if (ifld==3), 
% for ACC
%  j1=2; j2=10; i1=73; if (atlanornot==1), i1=i1-45; end
      lat1=-83; lat2=-50; if (atlanornot==0), long1=290; else; long1=-70; end
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      transp_accd=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
%
%  j1=2; j2=15; i1=6;  if (atlanornot==1), i1=i1+45; end
      lat1=-83; lat2=-30; long1=22; 
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      transp_accg=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
%
% j1=2; j2=13; i1=37; if (atlanornot==1), i1=i1+45; end
      lat1=-83; lat2=-38; long1=146; 
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      transp_accf=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
% for mediterranean
% j1=31; j2=33; i1=nx; if (atlanornot==1), i1=i1-45; end
      lat1= 34; lat2= 42; long1=178; 
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear tmp2 fi, tmp2=field(i1,j1:j2,:).*area_xzg(i1,j1:j2,:);
      fi=find(tmp2<0); transp_med=sum(tmp2(fi));
% indonesian throughflow
     if (bathynum==12 | bathynum==13), 
      lat1=-6; lat2=6; long1=142; 
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear tmp2 fi, tmp2=field(i1,j1:j2,:).*area_xzg(i1,j1:j2,:);
      fi=find(tmp2<0); transp_indtf=sum(tmp2(fi));
     end % if bathynum 12
   end % if
% meridional transports
   if (ifld==4),
% equatorial atlantic
%  j1=22; i1=80; i2=85; if (atlanornot==1), i1=i1-45; i2=i2-45; end
      lat1=-2; if (atlanornot==0), long1=318; long2=338; else
                                      long1=-42; long2=-22; end
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
      transp_eqatl=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
% gulf stream
%  j1=30; i1=70; i2=75; if (atlanornot==1),i1=i1-45; i2=i2-45; end
      lat1=30; if (atlanornot==0), long1=278; long2=298; else
                                      long1=-82; long2=-62; end
      if (bathynum==12 | bathynum==13), if (atlanornot==0), long1=310; long2=318;
                         else, long1=-50; long2=-32; end, end
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
      transp_gstre=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
% kuroshio
%  j1=29; i1=31; i2=36; if (atlanornot==1), i1=i1+45; i2=i2+45; end
      lat1=26; long1=122; long2=142; 
       if (bathynum==12 | bathynum==13), long1=142; long2=162; end
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
      transp_kuros=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
% indonesian throughflow
%  j1=22; i1=30; i2=35; if (atlanornot==1), i1=i1+45; i2=i2+45; end
     if (bathynum~=12 & bathynum~=13)
      lat1=-2; long1=118; long2=138;                        
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
      transp_indtf=sum(sum(field(i1:i2,j1,:).*area_xzg(i1:i2,j1,:) ));
     end % if bathynum 12
% arctic exchange
%  j1=39; i1=1; i2=7; if (atlanornot==1), i1=83-45; i2=7+45; end
      lat1=66; long2=26; if (atlanornot==0), long1=2; else, long1=-30; end
      if (bathynum==12 | bathynum==13), 
        if (atlanornot==0), long1=2; else, long1=-50; end, end
      clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
      clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
      clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
      clear tmp, tmp=field.*area_xzg;
      clear tmp2 fi, tmp2=tmp(i1:i2,j1,:); fi=find(tmp2>0); 
      transp_arcex=sum(tmp2(fi));
      if (atlanornot==0), long1=330; long2=360;
         if (bathynum==12 | bathynum==13), long1=310; end
         clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
         clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); 
           if (isempty(fi)==1), i2=nx; else, i2=fi(1); end
         clear tmp2 fi, tmp2=tmp(i1:i2,j1,:); fi=find(tmp2>0);
         transp_arcex= transp_arcex+sum(tmp2(fi));
      end % if
   end % if
%%%%%%%%%%%%%%%%%%%%
%     XY FIGURES   %
%%%%%%%%%%%%%%%%%%%%
   clear tmp, if (ifld==6); tmp=field; else, tmp(:,:)=field(:,:,kn); end
   clear fi, fi=find(pmask(:,:,kn)~=1); tmp(fi)=NaN;
   if (ifld<7), 
     figure(10), subplot(3,2,ifld), pcolor(long,lat,tmp'); 
   else
     figure(11), hold off, subplot(1,1,1), pcolor(long,lat,tmp');
   end
   shading flat, colorbar
   hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
   set(han,'LineWidth',1.5);
   if (ifld==6), title([direc,': ',str]); 
   else, title([str,' at ',num2str(depthm(kn)),'m']); end
   if (ifld==5), title([str,' at ',num2str(depthkn1),'m']); end
%
% plot difference from Levitus (for salt and temp)
   if (ifld==1 | ifld==2), clear tmp2
    if (ifld==1), tmp2=tmp-theta_clim(:,:,kn); end
    if (ifld==2), tmp2=tmp-salt_clim(:,:,kn); end
    figure(30+ifld), hold off, subplot(1,1,1), 
    subplot(3,1,1), pcolor(long,lat,tmp2'); shading flat, colorbar
    hold on, contour(long,lat,tmp2',[-1000 0 1000],'k-');
    title([titstr,' minus ',str,'_{lev}, at ',num2str(depthm(kn)),'m']);
   end % if 
% **********************PLOT FGN*******************************************
   if (ifld==1 | ifld==2), figure(fgn+2), subplot(3,2,ifld), 
        pcolor(long,lat,tmp'); caxis([cax1 cax2]); shading flat, colorbar
        hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
        set(han,'LineWidth',1.5);
        contour(long,lat,tmp',[cint1:cint:cint2],'k-');
        title([str,' at ',num2str(depthm(kn)),'m']);
        axis([ax1 ax2 ax3 ax4]);
        mintr=min(min(tmp')); maxtr=max(max(tmp'));
        han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
        set(han,'FontSize',6);
% title for page fgn+2
       if (ifld==2), 
         han=text(ax1-8*inton,ax4+3*inton,[titstr]); set(han,'FontSize',12);
         text(ax1-4*inton,ax4+2*inton,['year ', num2str(timeyear)]);
       end
   end
   if (ifld==6)  
       figure(fgn), subplot(3,2,6),pcolor(long,lat,tmp'); 
       caxis([cax1 cax2]); shading flat, colorbar
       hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
       set(han,'LineWidth',1.5);
       contour(long,lat,tmp',[cint1:cint:cint2],'k-');
       title([str]); axis([ax1 ax2 ax3 ax4]);
       mintr=min(min(tmp')); maxtr=max(max(tmp'));
       han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
       set(han,'FontSize',6);
   end % if ifld
% subplot
   if (ifld==1 | ifld==2 | ifld==5); 
     inum=ifld+2; if (ifld==5), inum=2; end
     clear tmp_sp, tmp_sp=field(:,:,knsp); 
     clear fi, fi=find(pmask(:,:,knsp)~=1); tmp_sp(fi)=NaN;
     figure(fgn+3), subplot(2,2,inum), pcolor(long,lat,tmp_sp');
     shading flat, colorbar
     hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
     set(han,'LineWidth',1.5);
% title for page 
     if (ifld==5), 
        han=text(ax1_sp-8*inton_sp,ax4_sp+3*inton_sp,[titstr]); 
        set(han,'FontSize',12);
        text(ax1_sp-4*inton_sp,ax4_sp+2*inton_sp,['year ', num2str(timeyear)]);
     end
     if (ifld==5), title([str,' at ',num2str(depthm(knsp)),'m']);
     else, title([str,' at ',num2str(depth(knsp)),'m']); end
     axis([ax1_sp ax2_sp ax3_sp ax4_sp]);
     mintr=min(min(tmp')); maxtr=max(max(tmp'));
     han=text(ax1_sp,ax3_sp-inton,[num2str(mintr),' to ',num2str(maxtr)]);
     set(han,'FontSize',6);
   end % if ifld
% ***************************************************************************
% plot velocity vectors
   if (ifld==4),
     for ip=1:3,
      if (ip==1), figure(12), hold off, subplot(1,1,1); knn=kn; vmultn=vmult; end
      if (ip==2), figure(fgn), subplot(2,1,1);  knn=kn; vmultn=vmult;end
      if (ip==3), figure(fgn+3), subplot(2,2,1);  knn=knsp; vmultn=vmult_sp; end
      clear tmp, tmp=temp_pmask(:,:,knn);   
      contour(longc,latc,tmp',1,'k-'); hold on     % draw land
      if (nx>105), nv=2; else, nv=1; end
      for i=1:nv:nx, for j=1:nv:ny1,
         if (i<nx1), ix=i+1; else ix=1; end
         u(i,j)=(uvel(i,j,knn)+uvel(ix,j,knn))/2;
         v(i,j)=(vvel(i,j,knn)+vvel(i,j+1,knn))/2;
         spd(i,j)=sqrt( u(i,j)^2 + v(i,j)^2 );
         if (spd(i,j)~=0 ) 
          if (u(i,j)~=0),
           dir=(180*atan(v(i,j)/u(i,j))/pi) ;
           if (u(i,j) < 0), dir = 180 + dir ; end,
          else
           if (v(i,j)>0) dir=90; else, dir=270; end
          end
          if (spd(i,j)>0),
           velvct(longc(i),latc(j),spd(i,j),dir,vmultn,'k'),
          end,
        end % if
     end, end %for i j
     title([titstr,' Velocity at ',num2str(depthm(knn)),'m']); 
     velvct(50,-75,.2,0,vmultn,'k');
     text(50,-70,'20cm/s');
%*********************** PLOT FGN *********************************************
     if (ip==2),
       axis([ax1 ax2 ax3 ax4]);
       title(['Velocity at ',num2str(depthm(kn)),'m']); 
       han=text(ax1+(ax2-ax1)/8,ax3-2*inton,[titstr]); set(han,'FontSize',18);
       han=text(ax1+(ax2-ax1)/8,ax3-3*inton,['year ', num2str(timeyear)]);
           set(han,'FontSize',14);
       text(ax1+(ax2-ax1)/8,ax3-3.5*inton,['dt=',num2str(dtday),' day']); 
       maxtr=max(max(spd(:,:))); han=text(ax1,ax3-inton,['max ',num2str(maxtr)]);
       set(han,'FontSize',6);
       c=axis; inton1=c(2)-c(1); inton2=c(4)-c(3);
% write transport numbers
       text(ax1+2*(ax2-ax1)/4,ax3-1.2*inton,'TRANSPORTS (Sv)')
       text(ax1+2*(ax2-ax1)/4,ax3-1.5*inton,'---------------')
       text(ax1+2*(ax2-ax1)/4,ax3-2*inton, ...
            ['ACC: ',num2str(transp_accd/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-2.5*inton, ...
            ['Arctic Exchange: ',num2str(transp_arcex/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-3*inton, ...
            ['Equatorial Crossing (Atlantic): ',num2str(transp_eqatl/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-3.5*inton, ...
            ['Gulf Stream (24N): ',num2str(transp_gstre/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-4*inton, ...
            ['Kuroshio (20N): ',num2str(transp_kuros/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-4.5*inton, ...
            ['Indonesian Throughflow: ',num2str(transp_indtf/1e6)]);
       text(ax1+2*(ax2-ax1)/4,ax3-5*inton, ...
            ['Mediterranean Exchange: ',num2str(transp_med/1e6)]);
     end % if %ip 2
     if (ip==3),
       axis([ax1_sp ax2_sp ax3_sp ax4_sp]);
       title(['Velocity at ',num2str(depthm(knsp)),'m']);
       maxtr=max(max(spd(:,:))); 
       han=text(ax1_sp,ax3_sp-inton_sp,['max ',num2str(maxtr)]);
       set(han,'FontSize',6);
     end % if
%*******************************************************************************
    end % for ip
   end % if ifld 4
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% meridional transects
 if (ifld<6 | ifld==7),
  num=0; if (ocmipornot==1), nocean=5; else, nocean=6; end
         if (bathynum==11), nocean=5; end
  clear tmp2 tmp3 tmp4
  tmp2=field.*volc; tmp2(filand_temp)=0;
  field(filand_temp)=NaN;
  theta_clim(filand_temp)=NaN;
  salt_clim(filand_temp)=NaN;
  for iocean=1:nocean,
      if (iocean==1), ioc=0; tstr=['Global']; end
      if (iocean==2), ioc=3; tstr=['Atlantic']; end
      if (iocean==3), ioc=2; tstr=['Pacific']; end
      if (iocean==4), ioc=1; tstr=['Indian']; end
      if (iocean==5), ioc=4; tstr=['Southern']; end
      if (iocean==6), ioc=5; tstr=['Arctic']; end
      if (iocean==7), ioc=6; tstr=['Med']; end
      for j=1:ny, for k=1:nz,
        clear fi,
        if (iocean==1), fi=find(tocean(:,j)>ioc & temp_pmask(:,j,k)==1); end;  
        if (iocean>1 & iocean <5), 
              fi=find(tocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
        if (iocean>4) fi=find(cocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
        if (isempty(fi)==0),
         trfield(j,k)=sum(field(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
        else, 
         trfield(j,k)=NaN; 
        end
      end, end % for j k
      figure(20+ifld), 
      if (iocean<5), subplot(3,2,iocean); else, subplot(3,3,iocean+2); end
      pcolor(lat,-depth,trfield'), shading flat, colorbar,
      hold on, c=contour(lat,-depth,trfield','k-'); clabel(c);
      [c,han]=contour(lat,-depth,trfield',[-1000 0 1000],'k-');
      set(han,'LineWidth',1.5); 
      if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
      if (iocean==5), axis([latg(1) 0 -5300 0]); end
      if (iocean==6), axis([60 latg(ny) -5300 0]); end
      if (iocean==7), axis([30 50 -5300 0]); end 
      title([tstr,': ',str]);
% plot Levitus and differences from Levitus
     if (ifld==1 | ifld==2), clear tmp2 lev,
      for j=1:ny, for k=1:nz,
        clear fi,
        if (iocean==1), fi=find(tocean(:,j)>ioc & temp_pmask(:,j,k)==1); end;
        if (iocean>1 & iocean <5),
              fi=find(tocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
        if (iocean>4) fi=find(cocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
        if (isempty(fi)==0),
         if (ifld==1), 
           lev(j,k)=sum(theta_clim(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
         else
           lev(j,k)=sum(salt_clim(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
         end
        else,
         lev(j,k)=NaN;
        end
      end, end % for j k
      figure(2+ifld), subplot(3,4,iocean+4);
      pcolor(lat,-depth,lev'), shading flat, colorbar,
      hold on, [c,han]=contour(lat,-depth,lev',[-1000 0 1000],'k-');
      set(han,'LineWidth',1.5); 
      if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
      if (iocean==5), axis([latg(1) 0 -5300 0]); end
      if (iocean==6), axis([60 latg(ny) -5300 0]); end
      if (iocean==7), axis([30 50 -5300 0]); end
      title([tstr,': ',str,'_{lev}']);
%
      figure(30+ifld), subplot(3,4,iocean+4);
      pcolor(lat,-depth,(trfield-lev)'), shading flat, colorbar,
      hold on, c=contour(lat,-depth,(trfield-lev)','k-'); clabel(c);
      [c,han]=contour(lat,-depth,(trfield-lev)',[-1000 0 1000],'k-');
      set(han,'LineWidth',1.5);
      if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
      if (iocean==5), axis([latg(1) 0 -5300 0]); end
      if (iocean==6), axis([60 latg(ny) -5300 0]); end
      if (iocean==7), axis([30 50 -5300 0]); end
      title([tstr]);
     end % if
% find single profile of differences
     if (ifld==1 | ifld==2), 
      clear prof_field profclim
      for k=1:nz, clear tmp3 tmp4 fi,
       field(filand_temp)=NaN;
       tmp3=volc(:,:,k); tmp4=field(:,:,k).*volc(:,:,k);
       if (ifld==1), tmp5=theta_clim(:,:,k).*volc(:,:,k); end
       if (ifld==2), tmp5=salt_clim(:,:,k).*volc(:,:,k); end
       if (iocean==1), fi=find(tocean>ioc & isnan(tmp4)==0); else,
                       fi=find(cocean==ioc & isnan(tmp4)==0); end
       prof_field(ifld,iocean,k)=sum(tmp4(fi))/sum(tmp3(fi));
       prof_clim(ifld,iocean,k)=sum(tmp5(fi))/sum(tmp3(fi));
       if (k==nz),
        figure(32+ifld), if (iocean==1), hold off, subplot(1,1,1); end 
        subplot(2,7,iocean), clear tmp, tmp=squeeze(prof_clim(ifld,iocean,:));
               plot(tmp,-depthm,'k--');
               hold on,      clear tmp, tmp=squeeze(prof_field(ifld,iocean,:));
               plot(tmp,-depthm,'r-');
        title([tstr]);
        clear tmp, 
        tmp=squeeze(prof_field(ifld,iocean,:)-prof_clim(ifld,iocean,:));
        subplot(2,7,iocean+7), plot(tmp,-depthm,'k-'); hold on,
                        plot([0 0],[-5300 0],':');
%****************************** PLOTFGN ******************************
        if (iocean==1), 
          figure(fgn+2), subplot(3,2,4+ifld);
          plot(tmp,-depthm,'k-'); hold on,
                        plot([0 0],[-5300 0],':');
          title([tstr,': ',str,'-',str,'_{lev}']);
        end % if iocean
%*********************************************************************
       end % if k nz
      end, % for k
     end % if ifld
  
%%******************PLOT FGN******************************
      if (iocean==1 & (ifld==1 | ifld==2) ),
       if (ifld==1), figure(fgn+2), pnum=3; end
       if (ifld==2), figure(fgn+2), pnum=4; end
          subplot(3,2,pnum), pcolor(lat,-depth',trfield'); 
          if (ifld==2), caxis([32 36]); else
                        caxis([cax1 cax2]); end
          shading flat, colorbar, hold on
          [c,han]=contour(lat,-depth',trfield',[-1000 0 1000],'k-');
          set(han,'LineWidth',1.5);
          contour(lat,-depth',trfield',[cint1:cint:cint2],'k');
          axis([latg(1) latg(ny) -5300 0]);
          mintr=min(min(trfield')); maxtr=max(max(trfield'));
          han=text(latg(1),-5500,[num2str(mintr),' to  ',num2str(maxtr)]);
          set(han,'FontSize',6);
          title(['Global, zonal ',str]);
     end % if
%***************************************************************
    end % for iocean
  end % if ifld
 end % for ifld
end % if choice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% merdional overturning
  dxelk=dxcg(1,:); dyelk=dyg(:,1);
  if (choice==1),
         str=['V'];
  else,
         str=['vVeltave'];
         str2=['GM_Kwy-T'];
  end
  field1=vvel;
  fid=fopen([run1str,'/',str2,'.',timestr,'.data'],'r','b');
  clear field3, field3=fread(fid,'float32');
  field3=reshape(field3,nx,ny,nz); 
  if (atlanornot==1), field3=turn3(olongv,field3); end
%
% calculate bolus velocity
  psiby=field3*kgm/2;
  for k=1:nz-1,
       vbc(:,:,k)=-(psiby(:,:,k)-psiby(:,:,k+1))/dz(k);
  end % for k
  vbc(:,:,nz)=(psiby(:,:,nz)-0)/dz(nz);
  for i=1:nx,
     for j=1:ny-1,
       vbg(i,j,:)=(vbc(i,j,:)+vbc(i,j+1,:))/2;
     end, % for j
      vbg(i,ny,:)=(vbc(i,ny,:)+vbc(i,1,:))/2;
   end % for i
   field3=vbg;
   if (atlanornot==0), splocean=9; else, splocean=9; end
 for iv=1:2,
  for k=1:nz, vocean(:,:,k)=cocean.*vmask(:,:,k); end
  for iocean=1:4,
      clear dxtab field
      if (iocean==1), ioc=0; tstr=['Global']; cax1=-30; cax2=30;  end
      if (iocean==2), ioc=3; tstr=['Atlantic']; end
      if (iocean==3), ioc=2; tstr=['Pacific']; end
      if (iocean==4), ioc=1; tstr=['Indian']; end
      clear fi fj,
      if (iocean==1), [fi]=find(vocean>ioc); else,
                                  [fi]=find(vocean==ioc); end
      mask=zeros(size(vocean)); mask(fi)=1;
      if (iv==1), clear field2, field2=field1.*mask; end
      if (iv==2), clear field2, field2=field3.*mask; end 
      clear fi fj fi1 fj1
      clear minfi maxfi minfj maxfj minfi1 maxfi1 mfi mfj mfi1 mfj1
      if (ioc~=splocean)
       [fi,fj]=find(vocean(:,:,1)==ioc);
       minfi=min(fi); maxfi=max(fi); mfi=maxfi-minfi+1;
       minfj=min(fj); maxfj=max(fj); mfj=maxfj-minfj+1;
       field=field2(minfi:maxfi,minfj:maxfj,:);
       lat=latc(minfj:maxfj);
       for k=1:nz,
         dxtab(:,:,k) = dxcg(minfi:maxfi,minfj:maxfj);
       end % for k
       dzglob2 =  ones(size(1:mfj))' * dz;
      else
       hnx=ceil(nx/2);
       [fi,fj]=find(vocean(hnx:nx,:,1)==ioc);
       [fi1,fj1]=find(vocean(1:hnx,:,1)==ioc);
       if (length(fi)>0),
        minfi=min(fi)+hnx-1; maxfi=max(fi)+hnx-1; mfi=maxfi-minfi+1;
        minfj=min(fj); maxfj=max(fj); mfj=maxfj-minfj+1;
       else
        minfi=NaN; maxfi==NaN; mfi=0;
       end % if
       if (length(fi1)>0),
        minfi1=min(fi1); maxfi1=max(fi1); mfi1=maxfi1-minfi1+1;
        minfj=min(min(fj1),minfj); maxfj=max(max(fj1),maxfj);
                                   mfj=maxfj-minfj+1;
       end % if
       field=field2(minfi:maxfi,minfj:maxfj,:);
       for k=1:nz,
          dxtab(:,:,k) = dxcg(minfi:maxfi,minfj:maxfj);
       end % for k
       if (length(fi1)>0),
         field(mfi+1:mfi+mfi1,1:mfj,:)= ...
                        field2(minfi1:maxfi1,minfj:maxfj,:);
         for k=1:nz,
          dxtab(mfi+1:mfi+mfi1,1:mfj,k) = dxcg(minfi1:maxfi1,minfj:maxfj);
         end % for k
       end % if isempty
       lat=latc(minfj:maxfj);
       dzglob2 =  ones(size(1:mfj))' * dz;
      end % ioc 3
% calculate streamfunction
%elk% =============================
%elk% 1. Integrate along longitude
%elk% =============================
%    tabx = dxtab.* squeeze( sum(field, 1) );
     tabx = squeeze( sum(field.*dxtab,1) );
%elk% ===================================
%elk% 2.  Integrate along depth
%elk% ===================================
     psi = NaN*zeros(size(tabx));
     for iz = nz:-1:1
         psi(:, iz) =  ...
                 -sum( dzglob2(:,nz:-1:iz) .* tabx(:,nz:-1:iz),2 );
     end
     psi(find(psi==0))=psi(find(psi==0))*NaN;
     psi=psi/1e6;
%%%%%%%
     figure(40+iv), if (iocean==1), hold off, subplot(1,1,1), end
     subplot(2,2,iocean);
     pcolor(lat,-depth',psi'); caxis([-30 30]), shading flat, 
                    colorbar, hold on
     c=contour(lat,-depth',psi',[-50:5:0],'k--'); clabel(c);
     hold on
     c=contour(lat,-depth',psi',[0:5:50],'k-'); clabel(c);
     [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
     set(han,'LineWidth',1.5);
     axis([latg(1) latg(ny) -5300 0]);
     if (iv==1), title([tstr,' Eulerian Overturning (Sv)']); end
     if (iv==2), title([tstr,' Bolus Overturning (Sv)']); end
     if (iv==1), eval(['psi_eu',num2str(iocean),'=psi;']);  end
     if (iv==2), psi_bo=psi;  end
%*******************PLOT FGN*****************************
     if (iocean==1),
       if (iv==1),ipl=1; end
       if (iv==2),ipl=3; end
       figure(fgn+1), subplot(3,2,ipl), pcolor(lat,-depth',psi'); 
         caxis([cax1 cax2]), shading flat,
                    colorbar, hold on
         contour(lat,-depth',psi',[-200:5:0],'k--'); 
         hold on
         contour(lat,-depth',psi',[0:5:200],'k-'); 
         [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
         set(han,'LineWidth',1.5);
         axis([latg(1) latg(ny) -5300 0]);
         mintr=min(min(psi')); maxtr=max(max(psi'));
         han=text(latg(1),-5500,[num2str(mintr),' to  ',num2str(maxtr)]);
         set(han,'FontSize',6);
         if (iv==1), title([tstr,' Eulerian Overturning (Sv)']); end
         if (iv==2), title([tstr,' Bolus Overturning (Sv)']); end
      end % if icean 1
      if (iv==2), eval(['psi=psi_bo+psi_eu',num2str(iocean),';']);
           if (iocean==1), ipl=5; end
           if (iocean==2), ipl=2; end
           if (iocean==3), ipl=4; end
           if (iocean==4), ipl=6; end
           figure(fgn+1), subplot(3,2,ipl), pcolor(lat,-depth',psi');
           caxis([cax1 cax2]), shading flat,
                    colorbar, hold on
          contour(lat,-depth',psi',[-200:5:0],'k--');
          hold on
          contour(lat,-depth',psi',[0:5:200],'k-');
          [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
          set(han,'LineWidth',1.5);
          axis([latg(1) latg(ny) -5300 0]);
          mintr=min(min(psi')); maxtr=max(max(psi'));
          han=text(latg(1),-5500,[num2str(mintr),' to  ',num2str(maxtr)]);
          set(han,'FontSize',6);
          title([tstr,' Total Overturning (Sv)']); 
      end %if iv 2

%********************************************************
%
   end % for iocean
  end % for iv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% horizontal streamfunctions
  dyelk=dyg(:,2); 
  field1=uvel;
% include Bolus transport
    str2=['GM_Kwx-T'];
    fid=fopen([run1str,'/',str2,'.',timestr,'.data'],'r','b');
    clear field3, field3=fread(fid,'float32');
    field3=reshape(field3,nx,ny,nz); 
    if (atlanornot==1), field3=turn3(olongu,field3); end
    psibx=field3*kgm/2;
    for k=1:nz-1,
        ubc(:,:,k)=-(psibx(:,:,k)-psibx(:,:,k+1))/dz(k);
    end % for k
    ubc(:,:,nz)=(psibx(:,:,nz)-0)/dz(nz);
    for j=1:ny,
       for i=1:nx-1,
        ubg(i,j,:)=(ubc(i,j,:)+ubc(i+1,j,:))/2;
       end, % for i
       ubg(nx,j,:)=(ubc(nx,j,:)+ubc(1,j,:))/2;
    end % for j
    field3=ubg;
    field1=field1+field3;
% end % if
  clear field2, field2=field1.*umask;
% calculate barotropic streamfunction
  dztab = dyelk* ones(size(1:ny));
%elk% =============================
%elk% 1. Integrate along depth
%elk% =============================
     tabz=zeros(nx,ny);
     for k=1:nz, tabz=tabz+field2(:,:,k)*dz(k); end
%elk% ===================================
%elk% 2.  Integrate along meridian
%elk% ===================================
     psi = NaN*zeros(size(tabz));
     for iy = ny:-1:1
         psi(:, iy) =  ...
                 sum( dyg(:,ny:-1:iy) .* tabz(:,ny:-1:iy),2 );
     end
     clear fi, fi=find(umask(:,:,1)==0); psi(fi)=NaN;
     psi= psi/1e6;  %psi(1,:)=psi(nx,:); 
%%%%%%
    figure(43); pcolor(longc,latc,psi');  shading flat,
                     colorbar, hold on
    c=contour(longc,latc,psi',[-200:10:0],'k--'); clabel(c);
     hold on
    c=contour(longc,latc,psi',[0:10:200],'k-'); clabel(c);
     [c,han]=contour(longc,latc,psi',[1000 0 1000],'k-'); 
     set(han,'LineWidth',1.5);
     title(['Barotropic Streamfunction (Sv)']);
% *****************PLOTFGN*****************************
   figure(fgn), subplot(3,2,5), pcolor(longc,latc,psi');  
                caxis([cax1 cax2]); shading flat,
                     colorbar, hold on
     contour(longc,latc,psi',[-200:10:0],'k--'); hold on
     contour(longc,latc,psi',[0:10:200],'k-'); 
     [c,han]=contour(longc,latc,psi',[1000 0 1000],'k-');
     set(han,'LineWidth',1.5);
     title(['Barotropic Streamfunction ']);
     axis([ax1 ax2 ax3 ax4]);
     mintr=min(min(psi')); maxtr=max(max(psi'));
     han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
     set(han,'FontSize',6);
% *****************************************************
% heat flux
%*******************************************************************
    figure(44), hold off, subplot(1,1,1);
    for iocean=1:4,
      if (iocean==1), ioc=0; tstr=['Global']; end
      if (iocean==2), ioc=3; tstr=['Atlantic']; end
      if (iocean==3), ioc=2; tstr=['Pacific']; end
      if (iocean==4), ioc=1; tstr=['Indian']; end
% from adcroft
      for j=1:ny;
       jm1=max([j-1 1]);
       Fadv(:,j,:)=vvel(:,j,:).*(273.15+(temp(:,jm1,:)+temp(:,j,:))/2);
      end
      for j=1:ny;
       jm1=max([j-1 1]);
       Fpar(:,j,:)=vbg(:,j,:).*(temp(:,j,:)-temp(:,jm1,:));
      end
% end from adcroft
      clear fi,
      tmp=zeros(size(temp_pmask)); for k=1:nz, tmp(:,:,k)=cocean; end
      if (iocean==1), fi=find(tmp>ioc); else, fi=find(tmp==ioc); end
      clear tmp2, tmp2=zeros(size(tmp)); tmp2(fi)=1;
      area=area_xzg.*temp_pmask.*tmp2; 
      for j=1:ny-1, for k=1:nz,
        trfield(j,k)=sum((Fadv(:,j,k).*area(:,j,k)),1);
%       trfield(j,k)=sum((Fpar(:,j,k)+Fadv(:,j,k)).*area(:,j,k),1);
      end, end % for j k
%%%%
      clear fi, fi=find(isnan(trfield)==1); trfield(fi)=0;
      rho=1024; cp=4200.;
      tottrfield=sum(trfield,2)*rho*cp/1e15;
      clear fi, fi=find(tottrfield==0); tottrfield(fi)=NaN;
      subplot(2,2,iocean); plot(latc,tottrfield,'k-');
      title([tstr,': meridional heat transport']);
      xlabel('latitude'); ylabel('PW');
     end % for iocean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some statistics
clear tmp, tmp=temp-theta_clim;
max(max(max(tmp))), min(min(min(tmp))),
clear fi, fi=find(tmp<0); tmp(fi)=-tmp(fi);
clear fi, fi=find(tmp>0);
clear tmp1, tmp1=tmp(fi).*volc(fi); stmp1=sum(tmp1);
clear tmp2, tmp2=volc(fi); stmp2=sum(tmp2);
tdiff=stmp1/stmp2
clear tmp, tmp=salt-salt_clim;
max(max(max(tmp))), min(min(min(tmp))),
clear fi, fi=find(tmp<0); tmp(fi)=-tmp(fi);
clear fi, fi=find(tmp>0);
clear tmp1, tmp1=tmp(fi).*volc(fi); stmp1=sum(tmp1);
clear tmp2, tmp2=volc(fi); stmp2=sum(tmp2);
sdiff=stmp1/stmp2
%
clear tmp, tmp=temp.*volc;
clear fi, fi =find(temp==0); tmp(fi)=NaN;
clear fi, fi=find(isnan(tmp)==0);
clear tmp2 tmp3, tmp2=tmp(fi); tmp3=volc(fi);
theta_mean=sum(tmp2)/sum(tmp3)
%
clear tmp, tmp=salt.*volc;
clear fi, fi =find(salt==0); tmp(fi)=NaN;
clear fi, fi=find(isnan(tmp)==0);
clear tmp2 tmp3, tmp2=tmp(fi); tmp3=volc(fi);
salt_mean=sum(tmp2)/sum(tmp3)
%
temp_max=max(max(max(temp)))
temp_min=min(min(min(temp)))
salt_max=max(max(max(salt)))
clear tmp fi, tmp=salt; fi=find(salt==0); tmp(fi)=NaN; 
salt_min=min(min(min(tmp)))