[Mitgcm-support] Matlab Reference Geopotential Calculation

[mitgcm-support at dev.mitgcm.org]

Below is Jean-Michel's script for calculating
reference geopotential for AIM_5L setup.
Script integrates up hydrostatic relation for
a reference set of in situ atmospheric temps.

Chris

%-----------------------------------------------------------------------------
xc= -179:2:180 ; yc= -89:2:90 ; zc= 1.:1.:5. ;
xu= -180. : 2. : 178. ; yv = -90. : 2. : 88. ; zi= 0.5 :1: 5.5 ;
yvp= -90 : 2. : 90 ;
Rt=6.370;Cp=.001;degrad=pi/180.;g=9.81;kappa=2./7.;
delY=pi*Rt/64;
delXv=cos(yvp*degrad)*2*pi*Rt;
delXc=cos(yc*degrad)*2*pi*Rt;
%delXv=2*pi*Rt*ones(65,1); delXc=2*pi*Rt*ones(64,1); 
%---------------
delR=[1.E1, 2.5E1, 3.E1, 2.E1, 1.5E1];
Rw(6)=0.;
for kl=5:-1:1, Rw(kl)=Rw(kl+1)+delR(kl); Rc(kl)=0.5*(Rw(kl+1)+Rw(kl));
end
Rtot=Rw(1); zcp=10*Rc;zip=10*Rw;
%-----------------------------------------------------------------------------
for kl=1:5, factP(kl)=(Rc(kl)/Rtot)^kappa ; end
tRef=[281.089  289.321  307.191  331.377  432.413];
%atm_cp=1004.;
  ddRp=1004.*(1.-factP(1)) ;
  phiRef(1)=ddRp*tRef(1) ;
for kl=2:5,
  ddRp=1004.*(factP(kl-1)-factP(kl)) ;
  phiRef(kl)=phiRef(kl-1)+ddRp*(tRef(kl-1)+tRef(kl))*0.5 ;
end