0,0 → 1,745 |
|
PROGRAM def_anomaly |
|
c ************************************************************************ |
c * Define a PV anomaly and set it up for inversion. Additionally, the * |
c * lateral and upper/lower boundary values for potential temperature * |
c * and wind are set * |
c * Michael Sprenger / Spring 2005 * |
c ************************************************************************ |
|
implicit none |
|
c ------------------------------------------------------------------------ |
c Declaration of subroutine parameters |
c ------------------------------------------------------------------------ |
|
c Numerical epsilon |
real eps |
parameter (eps=0.01) |
|
c Variables for input Z file : height level |
character*80 in_zfn,in_varname |
real in_varmin(4),in_varmax(4),in_stag(4) |
integer in_vardim(4) |
real in_mdv |
integer in_ndim |
integer in_nx,in_ny,in_nz |
real in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy |
integer in_ntimes |
real in_aklev(1000),in_bklev(1000) |
real in_aklay(1000),in_bklay(1000) |
|
real in_time |
real in_pollon,in_pollat |
integer in_idate(5) |
integer in_nvars |
character*80 in_vnam(100) |
real,allocatable, dimension (:,:,:) :: in_field |
real,allocatable, dimension (:,:,:) :: in_filt |
real,allocatable, dimension (:,:,:) :: in_anom |
real,allocatable, dimension (:,:) :: in_x,in_y |
real,allocatable, dimension (:,:,:) :: th_anom |
real,allocatable, dimension (:,:,:) :: uu_anom,vv_anom |
|
c Filtering |
real f_xmin,f_xmax |
real f_ymin,f_ymax |
real f_zmin,f_zmax |
integer nfilt |
real boundxy,boundz |
real,allocatable, dimension (:,:) :: tmp |
|
c Auxiliary variables |
integer i,j,k,l,n |
integer cdfid,cstid |
integer ierr,stat |
character*80 cfn |
real xpt,ypt,zpt |
integer isok |
character*80 varname,cdfname |
character*80 parfile |
real distxy,distz,distz1,distz2 |
integer indx,indy,indz |
real sum |
integer count |
character*80 name |
|
c Externals |
real kink |
external kink |
|
c ------------------------------------------------------------------------ |
c Preparations |
c ------------------------------------------------------------------------ |
|
print*,'*********************************************************' |
print*,'* def_anomaly *' |
print*,'*********************************************************' |
|
c Set name of the PV variable |
in_varname='PV' |
|
c Read parameters |
open(10,file='fort.10') |
read(10,*) in_zfn |
read(10,*) name,f_xmin |
if ( trim(name).ne.'BOX_XMIN') stop |
read(10,*) name,f_xmax |
if ( trim(name).ne.'BOX_XMAX') stop |
read(10,*) name,f_ymin |
if ( trim(name).ne.'BOX_YMIN') stop |
read(10,*) name,f_ymax |
if ( trim(name).ne.'BOX_YMAX') stop |
read(10,*) name,f_zmin |
if ( trim(name).ne.'BOX_ZMIN') stop |
read(10,*) name,f_zmax |
if ( trim(name).ne.'BOX_ZMAX') stop |
read(10,*) name,nfilt |
if ( trim(name).ne.'NFILTER' ) stop |
read(10,*) name,boundxy |
if ( trim(name).ne.'BOUND_XY') stop |
read(10,*) name,boundz |
if ( trim(name).ne.'BOUND_Z') stop |
close(10) |
|
print*,'Filter domain, x: ', f_xmin,f_xmax |
print*,' y: ', f_ymin,f_ymax |
print*,' z: ', f_zmin,f_zmax |
print*,'Nfilt: ',nfilt |
print*,'Boundxy,z: ',boundxy,boundz |
|
c Get grid description for Z file : height levels |
call cdfopn(in_zfn,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
call getcfn(cdfid,cfn,ierr) |
if (ierr.ne.0) goto 998 |
call cdfopn(cfn,cstid,ierr) |
if (ierr.ne.0) goto 998 |
call getvars(cdfid,in_nvars,in_vnam,ierr) |
if (ierr.ne.0) goto 998 |
call getdef(cdfid,in_varname,in_ndim,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 998 |
in_nx =in_vardim(1) |
in_ny =in_vardim(2) |
in_nz =in_vardim(3) |
in_xmin=in_varmin(1) |
in_ymin=in_varmin(2) |
call getlevs(cstid,in_nz,in_aklev,in_bklev,in_aklay,in_bklay,ierr) |
if (ierr.ne.0) goto 998 |
call getgrid(cstid,in_dx,in_dy,ierr) |
if (ierr.ne.0) goto 998 |
in_xmax=in_xmin+real(in_nx-1)*in_dx |
in_ymax=in_ymin+real(in_ny-1)*in_dy |
call gettimes(cdfid,in_time,in_ntimes,ierr) |
if (ierr.ne.0) goto 998 |
call getstart(cstid,in_idate,ierr) |
if (ierr.ne.0) goto 998 |
call getpole(cstid,in_pollon,in_pollat,ierr) |
if (ierr.ne.0) goto 998 |
call clscdf(cstid,ierr) |
if (ierr.ne.0) goto 998 |
call clscdf(cdfid,ierr) |
|
c Allocate memory for grid description |
allocate(in_x(in_nx,in_ny),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array in_x ***' |
allocate(in_y(in_nx,in_ny),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array in_y ***' |
|
c Allocate memory for input and output fields |
allocate(in_field(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array in_field ***' |
allocate(in_anom(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array in_anom ***' |
allocate(in_filt(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array in_filt ***' |
allocate(th_anom(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array th_anom ***' |
allocate(uu_anom(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array uu_anom ***' |
allocate(vv_anom(in_nx,in_ny,in_nz),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array vvnom ***' |
|
c Allocate memory for filtering |
allocate(tmp(in_nx,in_ny),stat=stat) |
if (stat.ne.0) print*,'*** error allocating array tmp ***' |
|
c Read variables from file |
call cdfopn(in_zfn,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
varname=trim(in_varname) |
print*,'R ',trim(varname),' ',trim(in_zfn) |
call gettimes(cdfid,in_time,in_ntimes,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,in_time,0,in_field,ierr) |
if (ierr.ne.0) goto 998 |
|
varname='X' |
print*,'R ',trim(varname),' ',trim(in_zfn) |
call gettimes(cdfid,in_time,in_ntimes,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,in_time,0,in_x,ierr) |
if (ierr.ne.0) goto 998 |
|
varname='Y' |
print*,'R ',trim(varname),' ',trim(in_zfn) |
call gettimes(cdfid,in_time,in_ntimes,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,in_time,0,in_y,ierr) |
if (ierr.ne.0) goto 998 |
|
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
c ------------------------------------------------------------------------ |
c Get the average along the x axis: "zonal" average |
c ------------------------------------------------------------------------ |
|
print*,'A ',trim(in_varname) |
do k=1,in_nz |
do j=1,in_ny |
|
c Get the mean along x axis |
sum=0. |
count=0 |
do i=1,in_nx |
if (abs(in_field(i,j,k)-in_mdv).gt.eps) then |
sum=sum+in_field(i,j,k) |
count=count+1 |
endif |
enddo |
if (count.ge.1) then |
sum=sum/real(count) |
else |
sum=in_mdv |
endif |
|
c Get the difference relative to the "zonal" mean |
do i=1,in_nx |
xpt=in_x(i,j) |
ypt=in_y(i,j) |
zpt=in_aklay(k) |
if ( (xpt.ge.f_xmin).and. |
> (xpt.le.f_xmax).and. |
> (ypt.ge.f_ymin).and. |
> (ypt.le.f_ymax).and. |
> (zpt.ge.f_zmin).and. |
> (zpt.le.f_zmax)) then |
in_filt(i,j,k)=sum |
else |
in_filt(i,j,k)=in_field(i,j,k) |
endif |
enddo |
|
enddo |
enddo |
|
c ------------------------------------------------------------------------ |
c Get the anomaly field |
c ------------------------------------------------------------------------ |
c |
c Determine the difference between filtered and unfiltered field |
do i=1,in_nx |
do j=1,in_ny |
do k=1,in_nz |
if ( (abs(in_field(i,j,k)-in_mdv).gt.eps).and. |
> (abs(in_filt (i,j,k)-in_mdv).gt.eps)) then |
in_anom(i,j,k)=in_field(i,j,k)-in_filt(i,j,k) |
else |
in_anom(i,j,k)=in_mdv |
endif |
enddo |
enddo |
enddo |
|
c Take only the positive part of the difference (mdv is set to 0) |
do i=1,in_nx |
do j=1,in_ny |
do k=1,in_nz |
if ( (abs(in_anom(i,j,k)-in_mdv).lt.eps) ) then |
in_anom(i,j,k)=0. |
else if ( in_anom(i,j,k).lt.0.) then |
in_anom(i,j,k)=0. |
endif |
enddo |
enddo |
enddo |
|
c Smooth transition to zero at boundaries |
do i=1,in_nx |
do j=1,in_ny |
do k=1,in_nz |
if ( (abs(in_anom (i,j,k)-in_mdv).gt.eps) ) then |
|
in_anom(i,j,k)=in_anom(i,j,k)* |
> kink(f_zmax -in_aklay(k),boundz )* |
> kink(in_aklay(k)-f_zmin ,boundz )* |
> kink(in_x(i,j) -f_xmin ,boundxy)* |
> kink(f_xmax -in_x(i,j) ,boundxy)* |
> kink(in_y(i,j) -f_ymin ,boundxy)* |
> kink(f_ymax -in_y(i,j) ,boundxy) |
|
endif |
enddo |
enddo |
enddo |
|
c ------------------------------------------------------------------------ |
c Apply a digital filter to the anomaly field |
c ------------------------------------------------------------------------ |
|
print*,'F ',trim(in_varname) |
do k=1,in_nz |
|
do i=1,in_nx |
do j=1,in_ny |
tmp(i,j)=in_anom(i,j,k) |
enddo |
enddo |
|
do n=1,nfilt |
call filt2d (tmp,tmp,in_nx,in_ny,1.,in_mdv,0,0,0,0) |
enddo |
|
do i=1,in_nx |
do j=1,in_ny |
in_anom(i,j,k)=tmp(i,j) |
enddo |
enddo |
|
enddo |
|
c ------------------------------------------------------------------------ |
c Get the modified PV field (original minus anomaly) |
c ------------------------------------------------------------------------ |
|
do i=1,in_nx |
do j=1,in_ny |
do k=1,in_nz |
if ( (abs(in_anom (i,j,k)-in_mdv).gt.eps).and. |
> (abs(in_field(i,j,k)-in_mdv).gt.eps) ) then |
in_filt(i,j,k)=in_field(i,j,k)-in_anom(i,j,k) |
else |
in_filt(i,j,k)=in_mdv |
endif |
enddo |
enddo |
enddo |
|
c ----------------------------------------------------------------- |
c Define the boundary values for wind and potential temperature |
c ----------------------------------------------------------------- |
|
do i=1,in_nx |
do j=1,in_ny |
do k=1,in_nz |
th_anom(i,j,k)=0. |
uu_anom(i,j,k)=0. |
vv_anom(i,j,k)=0. |
enddo |
enddo |
enddo |
|
c ----------------------------------------------------------------- |
c Write output |
c ----------------------------------------------------------------- |
|
c Open output netcdf file |
call cdfwopn(in_zfn,cdfid,ierr) |
if (ierr.ne.0) goto 997 |
|
c Write the 3d filtered field |
isok=0 |
varname=trim(in_varname)//'_FILT' |
print*,'Write ',trim(varname),' ',trim(in_zfn) |
call check_varok(isok,varname,in_vnam,in_nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,4,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 997 |
endif |
call putdat(cdfid,varname,in_time,0,in_filt,ierr) |
if (ierr.ne.0) goto 997 |
|
c Write the 3d anomaly field |
isok=0 |
varname=trim(in_varname)//'_ANOM' |
print*,'Write ',trim(varname),' ',trim(in_zfn) |
call check_varok(isok,varname,in_vnam,in_nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,4,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 997 |
endif |
call putdat(cdfid,varname,in_time,0,in_anom,ierr) |
if (ierr.ne.0) goto 997 |
|
c Write lower and upper boundary |
isok=0 |
varname='TH_ANOM' |
print*,'Write ',trim(varname),' ',trim(in_zfn) |
call check_varok(isok,varname,in_vnam,in_nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,4,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 997 |
endif |
call putdat(cdfid,varname,in_time,0,th_anom,ierr) |
if (ierr.ne.0) goto 997 |
|
c Write lateral boundary |
isok=0 |
varname='UU_ANOM' |
print*,'Write ',trim(varname),' ',trim(in_zfn) |
call check_varok(isok,varname,in_vnam,in_nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,4,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 997 |
endif |
call putdat(cdfid,varname,in_time,0,uu_anom,ierr) |
if (ierr.ne.0) goto 997 |
isok=0 |
varname='VV_ANOM' |
print*,'Write ',trim(varname),' ',trim(in_zfn) |
call check_varok(isok,varname,in_vnam,in_nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,4,in_mdv,in_vardim, |
> in_varmin,in_varmax,in_stag,ierr) |
if (ierr.ne.0) goto 997 |
endif |
call putdat(cdfid,varname,in_time,0,vv_anom,ierr) |
if (ierr.ne.0) goto 997 |
|
c Close output netcdf file |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 997 |
|
|
c ----------------------------------------------------------------- |
c Exception handling and format specs |
c ----------------------------------------------------------------- |
|
stop |
|
998 print*,'Problems with input netcdf file' |
stop |
|
997 print*,'Problems with output netcdf file' |
stop |
|
END |
|
|
c ****************************************************************** |
c * 2D FILTERING * |
c ****************************************************************** |
|
c ----------------------------------------------------------------- |
c Horizontal filter |
c ----------------------------------------------------------------- |
|
subroutine filt2d (a,af,nx,ny,fil,misdat, |
& iperx,ipery,ispol,inpol) |
|
c Apply a conservative diffusion operator onto the 2d field a, |
c with full missing data checking. |
c |
c a real inp array to be filtered, dimensioned (nx,ny) |
c af real out filtered array, dimensioned (nx,ny), can be |
c equivalenced with array a in the calling routine |
c f1 real workarray, dimensioned (nx+1,ny) |
c f2 real workarray, dimensioned (nx,ny+1) |
c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1 |
c corresponds to one application of the linear filter. |
c misdat real inp missing-data value, a(i,j)=misdat indicates that |
c the corresponding value is not available. The |
c misdat-checking can be switched off with with misdat=0. |
c iperx int inp periodic boundaries in the x-direction (1=yes,0=no) |
c ipery int inp periodic boundaries in the y-direction (1=yes,0=no) |
c inpol int inp northpole at j=ny (1=yes,0=no) |
c ispol int inp southpole at j=1 (1=yes,0=no) |
c |
c Christoph Schaer, 1993 |
|
c argument declaration |
integer nx,ny |
real a(nx,ny),af(nx,ny),fil,misdat |
integer iperx,ipery,inpol,ispol |
|
c local variable declaration |
integer i,j,is |
real fh |
real f1(nx+1,ny),f2(nx,ny+1) |
|
c compute constant fh |
fh=0.125*fil |
|
c compute fluxes in x-direction |
if (misdat.eq.0.) then |
do j=1,ny |
do i=2,nx |
f1(i,j)=a(i-1,j)-a(i,j) |
enddo |
enddo |
else |
do j=1,ny |
do i=2,nx |
if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then |
f1(i,j)=0. |
else |
f1(i,j)=a(i-1,j)-a(i,j) |
endif |
enddo |
enddo |
endif |
if (iperx.eq.1) then |
c do periodic boundaries in the x-direction |
do j=1,ny |
f1(1,j)=f1(nx,j) |
f1(nx+1,j)=f1(2,j) |
enddo |
else |
c set boundary-fluxes to zero |
do j=1,ny |
f1(1,j)=0. |
f1(nx+1,j)=0. |
enddo |
endif |
|
c compute fluxes in y-direction |
if (misdat.eq.0.) then |
do j=2,ny |
do i=1,nx |
f2(i,j)=a(i,j-1)-a(i,j) |
enddo |
enddo |
else |
do j=2,ny |
do i=1,nx |
if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then |
f2(i,j)=0. |
else |
f2(i,j)=a(i,j-1)-a(i,j) |
endif |
enddo |
enddo |
endif |
c set boundary-fluxes to zero |
do i=1,nx |
f2(i,1)=0. |
f2(i,ny+1)=0. |
enddo |
if (ipery.eq.1) then |
c do periodic boundaries in the x-direction |
do i=1,nx |
f2(i,1)=f2(i,ny) |
f2(i,ny+1)=f2(i,2) |
enddo |
endif |
if (iperx.eq.1) then |
if (ispol.eq.1) then |
c do south-pole |
is=(nx-1)/2 |
do i=1,nx |
f2(i,1)=-f2(mod(i-1+is,nx)+1,2) |
enddo |
endif |
if (inpol.eq.1) then |
c do north-pole |
is=(nx-1)/2 |
do i=1,nx |
f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny) |
enddo |
endif |
endif |
|
c compute flux-convergence -> filter |
if (misdat.eq.0.) then |
do j=1,ny |
do i=1,nx |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
enddo |
enddo |
else |
do j=1,ny |
do i=1,nx |
if (a(i,j).eq.misdat) then |
af(i,j)=misdat |
else |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
endif |
enddo |
enddo |
endif |
|
end |
|
|
c ****************************************************************** |
c * AUXILIARY SUBROUTINES * |
c ****************************************************************** |
|
C ----------------------------------------------------------------- |
c Rene's KINK function (for smoothing at bounadries) |
C ----------------------------------------------------------------- |
|
real function kink (x,a) |
|
implicit none |
|
c declaration of parameters |
real x,a |
|
c parameters |
real pi |
parameter (pi=3.1415926535) |
|
if (x.lt.0.) then |
kink=0. |
elseif (x.gt.a) then |
kink=1. |
else |
kink=(1.+cos(pi*(x-a)/a))/2. |
endif |
|
return |
end |
|
|
c --------------------------------------------------------- |
c Subroutines to write the netcdf output file |
c --------------------------------------------------------- |
|
subroutine writecdf2D(cdfname, |
> varname,arr,time, |
> dx,dy,xmin,ymin,nx,ny, |
> crefile,crevar) |
|
c Create and write to the netcdf file <cdfname>. The variable |
c with name <varname> and with time <time> is written. The data |
c are in the two-dimensional array <arr>. The list <dx,dy,xmin, |
c ymin,nx,ny> specifies the output grid. The flags <crefile> and |
c <crevar> determine whether the file and/or the variable should |
c be created. 1: create / 0: not created |
|
IMPLICIT NONE |
|
c Declaration of input parameters |
character*80 cdfname,varname |
integer nx,ny |
real arr(nx,ny) |
real dx,dy,xmin,ymin |
real time |
integer crefile,crevar |
|
c Further variables |
real varmin(4),varmax(4),stag(4) |
integer ierr,cdfid,ndim,vardim(4) |
character*80 cstname |
real mdv |
integer datar(14),stdate(5) |
integer i |
|
C Definitions |
varmin(1)=xmin |
varmin(2)=ymin |
varmin(3)=1050. |
varmax(1)=xmin+real(nx-1)*dx |
varmax(2)=ymin+real(ny-1)*dy |
varmax(3)=1050. |
cstname=trim(cdfname)//'_cst' |
ndim=4 |
vardim(1)=nx |
vardim(2)=ny |
vardim(3)=1 |
stag(1)=0. |
stag(2)=0. |
stag(3)=0. |
mdv=-999.98999 |
|
C Create the file |
if (crefile.eq.0) then |
call cdfwopn(cdfname,cdfid,ierr) |
if (ierr.ne.0) goto 906 |
else if (crefile.eq.1) then |
call crecdf(cdfname,cdfid, |
> varmin,varmax,ndim,cstname,ierr) |
if (ierr.ne.0) goto 903 |
|
C Write the constants file |
datar(1)=vardim(1) |
datar(2)=vardim(2) |
datar(3)=int(1000.*varmax(2)) |
datar(4)=int(1000.*varmin(1)) |
datar(5)=int(1000.*varmin(2)) |
datar(6)=int(1000.*varmax(1)) |
datar(7)=int(1000.*dx) |
datar(8)=int(1000.*dy) |
datar(9)=1 |
datar(10)=1 |
datar(11)=0 ! data version |
datar(12)=0 ! cstfile version |
datar(13)=0 ! longitude of pole |
datar(14)=90000 ! latitude of pole |
do i=1,5 |
stdate(i)=0 |
enddo |
call wricst(cstname,datar,0.,0.,0.,0.,stdate) |
endif |
|
c Write the data to the netcdf file, and close the file |
if (crevar.eq.1) then |
call putdef(cdfid,varname,ndim,mdv, |
> vardim,varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 904 |
endif |
call putdat(cdfid,varname,time,0,arr,ierr) |
if (ierr.ne.0) goto 905 |
call clscdf(cdfid,ierr) |
|
return |
|
c Error handling |
903 print*,'*** Problem to create netcdf file ***' |
stop |
904 print*,'*** Problem to write definition ***' |
stop |
905 print*,'*** Problem to write data ***' |
stop |
906 print*,'*** Problem to open netcdf file ***' |
stop |
|
END |
|
|
c ------------------------------------------------------------------ |
c Auxiliary routines |
c ------------------------------------------------------------------ |
|
subroutine check_varok (isok,varname,varlist,nvars) |
|
c Check whether the variable <varname> is in the list <varlist(nvars)>. |
c If this is the case, <isok> is incremented by 1. Otherwise <isok> keeps |
c its value. |
|
implicit none |
|
c Declaraion of subroutine parameters |
integer isok |
integer nvars |
character*80 varname |
character*80 varlist(nvars) |
|
c Auxiliary variables |
integer i |
|
c Main |
do i=1,nvars |
if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
enddo |
|
end |
Property changes: |
Added: svn:executable |