[MITgcm-support] PV analysis in mitgcm

[Baylor Fox-Kemper]

Hi Michael,
   At the moment, there is not.  I was working on it a while back,  
and it is not the easiest thing to do consistently.  Essentially,  
since the model's vorticity equation and tracer equations are not  
discretized the same way, it is not always easy to figure out the  
fluxes.
   However, if you don't care about fluxes, you can calculate the pv  
offline:
   I've attached a matlab program that does an OK job.  The  
'tileopen' routine is also supplied, which allows you to open a  
snapshot of a variable across all multi-processor tiles, at least for  
simple cases and mnc output.  If you don't use mnc, you can just read  
in UVW, and density.
  The variable hydrofac switches between hydrostatic (hydrofac=0) and  
nonhydrostatic (hydrofac=1) versions of the PV.
   I haven't checked this into CVS because there are still a few  
funny things, and a few things that are specific to my model setup,  
but the rough idea is here.
Cheers,
    -Baylor

% This calculates the pv, relative vorticity, and small Ro pv for a  
given day
function  
[pv,pvloro,omegax,omegay,omegaz,Bx,By,Bz,xfoundb,yfoundb,zfoundb,xfoundu 
, yfoundv,zfoundw,U,V,W,B]=pvcalc(getday);

g=9.81;
tAlpha=2e-4;
rho0=1035;
f0=7.29e-5;

hydrofac=0;  %This determines whether the nonhydrostatic terms are  
included=1, else 0

[U,dayfoundu,xfoundu,yfoundu,zfoundu]=tileopen(getday,'U');
[V,dayfoundv,xfoundv,yfoundv,zfoundv]=tileopen(getday,'V');
[W,dayfoundw,xfoundw,yfoundw,zfoundw]=tileopen(getday,'W');

U=U(:,:,1:end-1);
xfoundu=xfoundu(1:end-1);
V=V(:,1:end-1,:);
yfoundv=yfoundv(1:end-1);

[nzu,nyu,nxu]=size(U);
[nzv,nyv,nxv]=size(V);
[nzw,nyw,nxw]=size(W);

zfoundw=[zfoundw-zfoundw(1);-300];

nz=length(zfoundw);
W(nz,:,:)=0;
[nzw,nyw,nxw]=size(W);


[B,dayfoundb,xfoundb,yfoundb,zfoundb]=tileopen(getday,'Temp');

B=g*tAlpha*B;

[nzb,nyb,nxb]=size(B);

% Do Derivatives

omegaz=zeros(nzu,nyv,nxu); %omegaz on zfoundu, yfoundv, xfoundu
omegay=zeros(nzw,nyw,nxu); %omegay on zfoundw, yfoundw, xfoundu
omegax=zeros(nzw,nyv,nxw); %omegax on zfoundw, yfoundv, xfoundw

Bz=zeros(nzw,nyb,nxb); %Bz on zfoundw, yfoundb, xfoundb
By=zeros(nzb,nyv,nxb); %By on zfoundb, yfoundv, xfoundb
Bx=zeros(nzb,nyb,nxu); %Bx on zfoundb, yfoundb, xfoundu

%y-derivatives
for jj=1:(length(yfoundw)-1)
omegax(:,jj+1,:)=hydrofac*(W(:,jj+1,:)-W(:,jj,:))/(yfoundw(jj+1)- 
yfoundw(jj));
omegaz(:,jj+1,:)=-(U(:,jj+1,:)-U(:,jj,:))/(yfoundu(jj+1)-yfoundu(jj));
By(:,jj+1,:)=(B(:,jj+1,:)-B(:,jj,:))/(yfoundb(jj+1)-yfoundb(jj));
end

%Slip,Insulating Wall in Y
omegax(:,1,:)=0;
omegaz(:,1,:)=0;
By(:,1,:)=0;
omegax(:,end,:)=0;
omegaz(:,end,:)=0;
By(:,end,:)=0;

% x-derivatives
for ii=1:(length(yfoundw)-1)
omegay(:,:,ii+1)=-hydrofac*(W(:,:,ii+1)-W(:,:,ii))/(xfoundw(ii+1)- 
xfoundw(ii));
omegaz(:,:,ii+1)=omegaz(:,:,ii+1)+(V(:,:,ii+1)-V(:,:,ii))/(xfoundv(ii 
+1)-xfoundv(ii));
Bx(:,:,ii+1)=(B(:,:,ii+1)-B(:,:,ii))/(xfoundb(ii+1)-xfoundb(ii));
end
%Periodicity in X
omegay(:,:,1)=-hydrofac*(W(:,:,1)-W(:,:,end))/(xfoundw(2)-xfoundw(1));
omegaz(:,:,1)=omegaz(:,:,1)+(V(:,:,1)-V(:,:,end))/(xfoundv(2)-xfoundv 
(1));
Bx(:,:,1)=(B(:,:,1)-B(:,:,end))/(xfoundb(2)-xfoundb(1));

%z-derivatives
for kk=1:(length(zfoundu)-1)
omegax(kk+1,:,:)=omegax(kk+1,:,:)-(V(kk+1,:,:)-V(kk,:,:))/(zfoundv(kk 
+1)-zfoundv(kk));
omegay(kk+1,:,:)=omegay(kk+1,:,:)+(U(kk+1,:,:)-U(kk,:,:))/(zfoundu(kk 
+1)-zfoundu(kk));
Bz(kk+1,:,:)=(B(kk+1,:,:)-B(kk,:,:))/(zfoundb(kk+1)-zfoundb(kk));
end

%Slip Walls in Z
omegax(1,:,:)=0;
omegay(1,:,:)=0;
Bz(1,:,:)=0;
omegax(end,:,:)=0;
omegay(end,:,:)=0;
Bz(end,:,:)=0;

% Regridding and multiplying
%omegaz on zfoundu, yfoundv, xfoundu  => zfoundb, yfoundb, xfoundb

newvar=zeros(size(B));
for ii=1:length(xfoundb)
for jj=1:length(yfoundb)
   newvar(:,jj,ii)=0.25*(omegaz(:,jj,mod(ii,nxb)+1)+omegaz(:,mod 
(jj,nyb)+1,mod(ii,nxb)+1)+omegaz(:,jj,ii)+omegaz(:,mod(jj,nyb)+1,ii));
end
end
omegaz=newvar;

%omegay on zfoundw, yfoundb, xfoundu  => zfoundb, yfoundb, xfoundb
newvar=zeros(size(B));
for ii=1:length(xfoundb)
for kk=1:length(zfoundb)
   newvar(kk,:,ii)=0.25*(omegay(kk+1,:,mod(ii,nxb)+1)+omegay(kk 
+1,:,ii)+omegay(kk,:,mod(ii,nxb)+1)+omegay(kk,:,ii));
end
end
omegay=newvar;

%omegax on zfoundw, yfoundv, xfoundb  => zfoundb, yfoundb, xfoundb
newvar=zeros(size(B));
for jj=1:length(yfoundb)
for kk=1:length(zfoundb)
   newvar(kk,jj,:)=0.25*(omegax(kk+1,mod(jj,nyb)+1,:)+omegax(kk 
+1,jj,:)+omegax(kk,mod(jj,nyb)+1,:)+omegax(kk,jj,:));
end
end
omegax=newvar;

%Bz on zfoundw, yfoundb, xfoundb  => zfoundb, yfoundb, xfoundb
newvar=zeros(size(B));
for kk=1:length(zfoundb)
newvar(kk,:,:)=0.5*(Bz(kk+1,:,:)+Bz(kk,:,:));
end
Bz=newvar;

%By on zfoundb, yfoundv, xfoundb  => zfoundb, yfoundb, xfoundb
newvar=zeros(size(B));
for jj=1:length(yfoundb)
newvar(:,jj,:)=0.5*(By(:,mod(jj,nyb)+1,:)+By(:,jj,:));
end
By=newvar;

%Bx on zfoundb, yfoundb, xfoundu  => zfoundb, yfoundb, xfoundb
newvar=zeros(size(B));
for ii=1:length(xfoundb)
newvar(:,:,ii)=0.5*(Bx(:,:,mod(ii,nxb)+1)+Bx(:,:,ii));
end
Bx=newvar;

pv=Bx.*omegax+By.*omegay+Bz.*(omegaz+f0);
pvloro=Bz*f0;






% Helps load data for a simulation on multiple tiles
% It finds the time closest to getday, and then loads the entire
% variable var.
% function [varout,dayfound,xfound,yfound,zfound]=tileopen(getday,var);
function [varout,dayfound,xfound,yfound,zfound]=tileopen(getday,var);

gdir=dir('grid*.nc');
sdir=dir('state.*.nc');
t1dir=dir('state.*.t001.nc');
if length(t1dir)==0
   t1dir=dir('state.*.000001.nc');
end

day=[];

% Time location
for getlength=1:(length(t1dir))
f=netcdf(t1dir(getlength).name);
if isempty(f), error([' ## Bad Netcdf File: ',t1dir 
(getlength).name]), end
daynow=f{'model_time',1}(:)/86400;
if length(daynow)==0
    daynow=f{'T',1}(:)/86400;
end
day(getlength,1:length(daynow))=daynow;
if (length(day(1,:))>length(daynow))
    day(getlength,(length(daynow)+1):end)=nan;
end
close(f)
end

[val,indx]=min(abs(day'-getday));
[dayfound,filenum]=min(val);
indx=indx(filenum);
dayfound=day(filenum,indx);

%Space Location

X=[];
Y=[];
varout=[];

for tile=1:(length(gdir))
f=netcdf(sdir((filenum-1)*length(gdir)+tile).name);

Varin=f{var}(indx,:,:,:);
thedims=dim(f{var});

Xin=f{name(thedims{4})}(:);
Yin=f{name(thedims{3})}(:);
zfound=f{name(thedims{2})}(:);

if length(find(X==Xin(end)))==0
   if length(find(X==Xin(1)))==0
    X=[X;Xin];
   else
    X=[X;Xin(2:end)];
   end
end
Xbeg=find(X==Xin(1));
Xend=find(X==Xin(end));

if length(find(Y==Yin(end)))==0
   if length(find(Y==Yin(1)))==0
    Y=[Y;Yin];
   else
    Y=[Y;Yin(2:end)];
   end
end
Ybeg=find(Y==Yin(1));
Yend=find(Y==Yin(end));

varout(1:length(zfound),Ybeg:Yend,Xbeg:Xend)=Varin;

close(f);
end

xfound=X;
yfound=Y;










On Mar 15, 2006, at 8:37 AM, Michael Beck wrote:

> good morning all
>
> Is there a way to diagnose or output Ertel's pv using the standard  
> version of the model using the diagnosis package?
>
> thanks
>
> michael beck
>
> ----------------------------------------------------------------
> This message was sent using IMP, the Internet Messaging Program.
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support