0,0 → 1,1468 |
PROGRAM z2s |
|
c Calculate secondary fields on z levels |
c Michael Sprenger / Summer 2006 |
|
implicit none |
|
c --------------------------------------------------------------- |
c Declaration of parameters |
c --------------------------------------------------------------- |
|
real tzero |
parameter (tzero=273.15) |
real kappa |
parameter (kappa=0.6078) |
real rd |
parameter (rd=287.) |
|
c --------------------------------------------------------------- |
c Declaration of variables |
c --------------------------------------------------------------- |
|
c Variables for output Z file : height level |
character*80 cfn |
real varmin(4),varmax(4),stag(4) |
integer vardim(4) |
real mdv |
integer ndim |
integer nx,ny,nz |
real xmin,xmax,ymin,ymax,dx,dy |
integer ntimes |
real aklev(1000),bklev(1000) |
real aklay(1000),bklay(1000) |
real time |
real pollon,pollat |
integer idate(5) |
integer nfields |
character*80 field_nam(100) |
real,allocatable, dimension (:,:,:,:) :: field_dat |
real,allocatable, dimension (:,:,:) :: z3 |
real,allocatable, dimension (:,:) :: x2,y2,f2,oro |
real,allocatable, dimension (:,:,:) :: out1,out2 |
real,allocatable, dimension (:,:,:) :: inp |
real,allocatable,dimension (:) :: nsqref |
real,allocatable,dimension (:) :: thetaref |
real,allocatable,dimension (:) :: tref |
real,allocatable,dimension (:) :: rhoref |
real,allocatable,dimension (:) :: pressref |
real,allocatable,dimension (:) :: zref |
real deltax,deltay,deltaz |
integer nvars |
character*80 vnam(100),varname |
integer isok |
integer cdfid,cstid |
character*3 mode |
|
c Parameter file |
character*80 fieldname |
integer nfilt |
character*80 ofn,gri |
|
c Auxiliary variables |
integer ierr,stat |
integer i,j,k,n |
real,allocatable, dimension (:,:) :: tmp |
character*80 vnam2(100) |
integer nvars2 |
|
c --------------------------------------------------------------- |
c Preparations |
c --------------------------------------------------------------- |
|
print*,'*********************************************************' |
print*,'* qvec_analysis *' |
print*,'*********************************************************' |
|
c Read parameter file |
open(10,file='fort.10') |
read(10,*) fieldname |
read(10,*) ofn |
read(10,*) gri |
read(10,*) nfilt |
close(10) |
|
c Decide whether to enter ANO or MOD/ORG mode |
mode = ofn(1:3) |
if ( (mode.ne.'ANO').and. |
> (mode.ne.'MOD').and. |
> (mode.ne.'ORG') ) |
>then |
print*,'Unsupported mode ',trim(mode) |
stop |
endif |
|
c Get grid description for Z file : height level |
call cdfopn(ofn,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 getdef(cdfid,'T',ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
nx =vardim(1) |
ny =vardim(2) |
nz =vardim(3) |
xmin=varmin(1) |
ymin=varmin(2) |
call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr) |
if (ierr.ne.0) goto 998 |
call getgrid(cstid,dx,dy,ierr) |
if (ierr.ne.0) goto 998 |
xmax=xmin+real(nx-1)*dx |
ymax=ymin+real(ny-1)*dy |
call gettimes(cdfid,time,ntimes,ierr) |
if (ierr.ne.0) goto 998 |
call getstart(cstid,idate,ierr) |
if (ierr.ne.0) goto 998 |
call getpole(cstid,pollon,pollat,ierr) |
if (ierr.ne.0) goto 998 |
call getvars(cdfid,nvars,vnam,ierr) |
if (ierr.ne.0) goto 998 |
call clscdf(cstid,ierr) |
if (ierr.ne.0) goto 998 |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
c Get a list of all variables on the GRID file |
if (trim(gri).ne.trim(ofn)) then |
call cdfopn(gri,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
call getvars(cdfid,nvars2,vnam2,ierr) |
if (ierr.ne.0) goto 998 |
do i=1,nvars2 |
vnam(nvars+i)=vnam2(i) |
enddo |
nvars=nvars+nvars2 |
endif |
|
c Check whether calculation of <fieldname> is supported |
if ( (fieldname.ne.'GEO' ).and. |
> (fieldname.ne.'AGEO' ).and. |
> (fieldname.ne.'DIV_UV' ).and. |
> (fieldname.ne.'QVEC' ).and. |
> (fieldname.ne.'DIV_Q' ) ) then |
print*,'Variable not supported ',trim(fieldname) |
stop |
endif |
|
c Set dependencies |
if (fieldname.eq.'GEO') then |
nfields=2 |
field_nam(1)='T' |
field_nam(2)='P' |
elseif (fieldname.eq.'AGEO') then |
nfields=4 |
field_nam(1)='U' |
field_nam(2)='UG' |
field_nam(3)='V' |
field_nam(4)='VG' |
else if (fieldname.eq.'DIV_UV') then |
nfields=2 |
field_nam(1)='U' |
field_nam(2)='V' |
else if (fieldname.eq.'QVEC') then |
nfields=3 |
field_nam(1)='UG' |
field_nam(2)='VG' |
field_nam(3)='TH' |
else if (fieldname.eq.'DIV_Q') then |
nfields=2 |
field_nam(1)='QX' |
field_nam(2)='QY' |
endif |
|
c Allocate memory |
allocate(field_dat(nfields,nx,ny,nz),stat=stat) |
if (stat.ne.0) print*,'error allocating field_dat' |
allocate(out1(nx,ny,nz),stat=stat) |
if (stat.ne.0) print*,'error allocating out1' |
allocate(out2(nx,ny,nz),stat=stat) |
if (stat.ne.0) print*,'error allocating out2' |
allocate(inp(nx,ny,nz),stat=stat) |
if (stat.ne.0) print*,'error allocating inp' |
allocate(z3(nx,ny,nz),stat=stat) |
if (stat.ne.0) print*,'error allocating z3' |
allocate(x2(nx,ny),stat=stat) |
if (stat.ne.0) print*,'error allocating x2' |
allocate(y2(nx,ny),stat=stat) |
if (stat.ne.0) print*,'error allocating y2' |
allocate(f2(nx,ny),stat=stat) |
if (stat.ne.0) print*,'error allocating f2' |
allocate(oro(nx,ny),stat=stat) |
if (stat.ne.0) print*,'error allocating oro' |
|
C Allocate memory for reference profile |
allocate(rhoref (0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating rhoref' |
allocate(pressref(0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating pressref' |
allocate(thetaref(0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating thetaref' |
allocate(nsqref (0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating nsqref' |
allocate(zref (0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating zref' |
allocate(tref (0:2*nz),STAT=stat) |
if (stat.ne.0) print*,'error allocating tref' |
|
c Allocate auxiliary fields |
allocate(tmp(nx,ny),stat=stat) |
if (stat.ne.0) print*,'error allocating tmp' |
|
c Read reference profile and grid parameters |
call read_ref (nsqref,rhoref,thetaref,pressref,zref, |
> nx,ny,nz,deltax,deltay,deltaz,f2,oro,gri) |
|
c Read X grid |
varname='X' |
isok=0 |
call check_varok (isok,varname,vnam,nvars) |
if (isok.eq.0) then |
print*,'Variable ',trim(varname),' missing... Stop' |
stop |
endif |
call cdfopn(gri,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,time,0,x2,ierr) |
if (ierr.ne.0) goto 998 |
print*,'R ',trim(varname),' ',trim(gri) |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
c Read Y grid |
varname='Y' |
isok=0 |
call check_varok (isok,varname,vnam,nvars) |
if (isok.eq.0) then |
print*,'Variable ',trim(varname),' missing... Stop' |
stop |
endif |
call cdfopn(gri,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,time,0,y2,ierr) |
if (ierr.ne.0) goto 998 |
print*,'R ',trim(varname),' ',trim(gri) |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
c Calculate reference profile of temperature |
print*,'C T_ref = TH_ref * ( P_ref/1000)^kappa' |
do i=0,2*nz |
tref(i) = thetaref(i) * ( pressref(i)/1000. ) ** kappa |
enddo |
|
c Init height levels |
do i=1,nx |
do j=1,ny |
do k=1,nz |
z3(i,j,k)=zref(2*k-1) |
enddo |
enddo |
enddo |
|
c Load needed variables |
do n=1,nfields |
|
c Check whether variable is available on file |
varname=field_nam(n) |
isok=0 |
call check_varok (isok,varname,vnam,nvars) |
|
c Variable is available on file |
if (isok.eq.1) then |
|
call cdfopn(ofn,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
call getdat(cdfid,varname,time,0,inp,ierr) |
if (ierr.ne.0) goto 998 |
print*,'R ',trim(varname),' ',trim(ofn) |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
do i=1,nx |
do j=1,ny |
do k=1,nz |
field_dat(n,i,j,k)=inp(i,j,k) |
enddo |
enddo |
enddo |
|
else |
print*,'Variable missing : ',trim(varname) |
stop |
endif |
|
enddo |
|
c Add reference profile if ANO file is provided |
if ( mode.eq.'ANO' ) then |
|
print*,'C T -> T_ano + T_ref' |
n=0 |
do i=1,nfields |
if (field_nam(i).eq.'T') n=i |
enddo |
if (n.ne.0) then |
do i=1,nx |
do j=1,ny |
do k=1,nz |
field_dat(n,i,j,k) = field_dat(n,i,j,k) + tref(2*k-1) |
enddo |
enddo |
enddo |
endif |
|
print*,'C TH -> TH_ano + TH_ref' |
n=0 |
do i=1,nfields |
if (field_nam(i).eq.'TH') n=i |
enddo |
if (n.ne.0) then |
do i=1,nx |
do j=1,ny |
do k=1,nz |
field_dat(n,i,j,k) = field_dat(n,i,j,k)+thetaref(2*k-1) |
enddo |
enddo |
enddo |
endif |
|
print*,'C P -> P_ano + P_ref' |
n=0 |
do i=1,nfields |
if (field_nam(i).eq.'P') n=i |
enddo |
if (n.ne.0) then |
do i=1,nx |
do j=1,ny |
do k=1,nz |
field_dat(n,i,j,k) = field_dat(n,i,j,k)+pressref(2*k-1) |
enddo |
enddo |
enddo |
endif |
|
endif |
|
c Change units: P (hPa -> Pa), T(K -> degC) |
n=0 |
do i=1,nfields |
if (field_nam(i).eq.'P') n=i |
enddo |
if (n.ne.0) then |
do i=1,nx |
do j=1,ny |
do k=1,nz |
field_dat(n,i,j,k)=100.*field_dat(n,i,j,k) |
enddo |
enddo |
enddo |
endif |
|
n=0 |
do i=1,nfields |
if (field_nam(i).eq.'T') n=i |
enddo |
if (n.ne.0) then |
do i=1,nx |
do j=1,ny |
do k=1,nz |
if ( field_dat(n,i,j,1).gt.100. ) then |
field_dat(n,i,j,k)=field_dat(n,i,j,k) - tzero |
endif |
enddo |
enddo |
enddo |
endif |
|
c ---------------------------------------------------------------- |
c Calculations |
c ---------------------------------------------------------------- |
|
c Geostrophic wind (GEO) |
if (fieldname.eq.'GEO') then |
|
print*,'C ',trim(fieldname) |
field_nam(1)='RHO' |
do i=1,nx |
do j=1,ny |
do k=1,nz |
if ( mode.eq.'ANO' ) then |
field_dat(1,i,j,k)=rhoref(2*k-1) |
else |
field_dat(1,i,j,k)= |
> 1./rd*field_dat(2,i,j,k)/(field_dat(1,i,j,k)+tzero) |
endif |
enddo |
enddo |
enddo |
call calc_geo (out1,out2, ! (Ug,Vg) |
> field_dat(1,:,:,:), ! RHO |
> field_dat(2,:,:,:), ! P |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL |
> nx,ny,nz,mdv) |
|
c Ageostrophic wind (AGEO) |
elseif (fieldname.eq.'AGEO') then |
|
print*,'C ',trim(fieldname) |
call calc_ageo (out1,out2, ! (Ua,Va) |
> field_dat(1,:,:,:), ! U |
> field_dat(2,:,:,:), ! UG |
> field_dat(3,:,:,:), ! V |
> field_dat(4,:,:,:), ! VG |
> nx,ny,nz,mdv) |
|
c Divergence of wind (DIV_UV) |
else if (fieldname.eq.'DIV_UV') then |
|
print*,'C ',trim(fieldname) |
call calc_div_uv (out1, ! Div(U,V)) |
> field_dat(1,:,:,:), ! U |
> field_dat(2,:,:,:), ! V |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL |
> nx,ny,nz,mdv) |
|
c Q vector (QVEC) |
else if (fieldname.eq.'QVEC') then |
|
print*,'C ',trim(fieldname) |
call calc_qvec (out1,out2, ! (Qx,Qy) |
> field_dat(1,:,:,:), ! UG |
> field_dat(2,:,:,:), ! VG |
> field_dat(3,:,:,:), ! TH |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL |
> nx,ny,nz,mdv) |
|
c Divergence of Q vector (DIV_Q) |
else if (fieldname.eq.'DIV_Q') then |
|
print*,'C ',trim(fieldname) |
call calc_div_q (out1, ! Div(U,V)) |
> field_dat(1,:,:,:), ! QX |
> field_dat(2,:,:,:), ! QY |
> z3,x2,y2,f2, ! Z,X,Y,CORIOL |
> nx,ny,nz,mdv) |
|
endif |
|
c Horizontal filtering of the resulting fields |
print*,'F ',trim(fieldname) |
do k=1,nz |
|
do i=1,nx |
do j=1,ny |
tmp(i,j)=out1(i,j,k) |
enddo |
enddo |
do n=1,nfilt |
call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0) |
enddo |
do i=1,nx |
do j=1,ny |
out1(i,j,k)=tmp(i,j) |
enddo |
enddo |
|
do i=1,nx |
do j=1,ny |
tmp(i,j)=out2(i,j,k) |
enddo |
enddo |
do n=1,nfilt |
call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0) |
enddo |
do i=1,nx |
do j=1,ny |
out2(i,j,k)=tmp(i,j) |
enddo |
enddo |
|
enddo |
|
c ---------------------------------------------------------------- |
c Save result onto netcdf file |
c ---------------------------------------------------------------- |
|
c Open output file |
call cdfwopn(ofn,cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
c Save geostrophic wind |
if (fieldname.eq.'GEO') then |
isok=0 |
varname='UG' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out1,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
isok=0 |
varname='VG' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out2,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
|
c Save ageostrophic wind |
elseif (fieldname.eq.'AGEO') then |
isok=0 |
varname='UA' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out1,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
isok=0 |
varname='VA' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out2,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
|
c Save divergence of wind field |
else if (fieldname.eq.'DIV_UV') then |
isok=0 |
varname='DIV_UV' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out1,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
|
c Save components of Q vector |
else if (fieldname.eq.'QVEC') then |
isok=0 |
varname='QX' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out1,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
isok=0 |
varname='QY' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out2,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
|
c Save divergence of wind field |
else if (fieldname.eq.'DIV_Q') then |
isok=0 |
varname='DIV_Q' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) then |
call putdef(cdfid,varname,ndim,mdv,vardim, |
> varmin,varmax,stag,ierr) |
if (ierr.ne.0) goto 998 |
endif |
call putdat(cdfid,varname,time,0,out1,ierr) |
if (ierr.ne.0) goto 998 |
print*,'W ',trim(varname),' ',trim(ofn) |
|
endif |
|
c Close output file |
call clscdf(cdfid,ierr) |
if (ierr.ne.0) goto 998 |
|
|
c ---------------------------------------------------------------- |
c Exception handling |
c ---------------------------------------------------------------- |
|
stop |
|
998 print*,'Problem with netcdf file' |
stop |
|
end |
|
|
|
c **************************************************************** |
c * SUBROUTINE SECTION: AUXILIARY ROUTINES * |
c **************************************************************** |
|
c ---------------------------------------------------------------- |
c Check whether variable is found on netcdf file |
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> |
c keeps 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 |
|
|
c **************************************************************** |
c * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS * |
c **************************************************************** |
|
c ---------------------------------------------------------------- |
c Calculate geostrophic wind components |
c ---------------------------------------------------------------- |
|
subroutine calc_geo (ug,vg,rho,p, |
> z3,x2,y2,f2,nx,ny,nz,mdv) |
|
c Calculate the geostrophic wind components (ug,vg) if the temperature |
c (t) and the pressure (p) are given. The grid and the missing data |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real ug (nx,ny,nz) |
real vg (nx,ny,nz) |
real rho(nx,ny,nz) |
real p (nx,ny,nz) |
real z3 (nx,ny,nz) |
real x2 (nx,ny) |
real y2 (nx,ny) |
real f2 (nx,ny) |
real mdv |
|
c Physical parameters and numerical constants |
real eps |
parameter (eps=0.01) |
real g |
parameter (g=9.80616) |
|
|
c Auxiliray variables |
integer i,j,k |
real dpdx(nx,ny,nz) |
real dpdy(nx,ny,nz) |
|
c Calculate horizontal derivatives of pressure |
call deriv(dpdx,p,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv(dpdy,p,'y',z3,x2,y2,nx,ny,nz,mdv) |
|
c Calculation |
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
if ((abs(rho (i,j,k)-mdv).gt.eps).and. |
> (abs(dpdx(i,j,k)-mdv).gt.eps).and. |
> (abs(dpdy(i,j,k)-mdv).gt.eps)) then |
|
ug(i,j,k)=-1./(rho(i,j,k)*f2(i,j))*dpdy(i,j,k) |
vg(i,j,k)= 1./(rho(i,j,k)*f2(i,j))*dpdx(i,j,k) |
|
endif |
|
enddo |
enddo |
enddo |
|
end |
|
c ---------------------------------------------------------------- |
c Calculate ageostrophic wind components |
c ---------------------------------------------------------------- |
|
subroutine calc_ageo (ua,va,u,ug,v,vg,nx,ny,nz,mdv) |
|
c Calculate the geostrophic wind components (ug,vg) if the temperature |
c (t) and the pressure (p) are given. The grid and the missing data |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real ua (nx,ny,nz) |
real va (nx,ny,nz) |
real u (nx,ny,nz) |
real ug (nx,ny,nz) |
real v (nx,ny,nz) |
real vg (nx,ny,nz) |
real mdv |
|
c Parameters |
real eps |
parameter (eps=0.01) |
|
c Auxiliray variables |
integer i,j,k |
|
c Calculation |
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
if ((abs(u (i,j,k)-mdv).gt.eps).and. |
> (abs(ug(i,j,k)-mdv).gt.eps).and. |
> (abs(v (i,j,k)-mdv).gt.eps).and. |
> (abs(vg(i,j,k)-mdv).gt.eps)) then |
|
ua(i,j,k)=u(i,j,k) - ug(i,j,k) |
va(i,j,k)=v(i,j,k) - vg(i,j,k) |
|
endif |
|
enddo |
enddo |
enddo |
|
end |
|
c ---------------------------------------------------------------- |
c Calculate divergence of wind field |
c ---------------------------------------------------------------- |
|
subroutine calc_div_uv (div_uv,u,v, |
> z3,x2,y2,f2,nx,ny,nz,mdv) |
|
c Calculate the divergence (div_uv) of the horizontal wind field if+ |
c the wind components (u,v) are given. The grid and the missing data |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real div_uv(nx,ny,nz) |
real u (nx,ny,nz) |
real v (nx,ny,nz) |
real z3 (nx,ny,nz) |
real x2 (nx,ny) |
real y2 (nx,ny) |
real f2 (nx,ny) |
real mdv |
|
c Physical parameters and numerical constants |
real eps |
parameter (eps=0.01) |
|
c Auxiliray variables |
integer i,j,k |
real rho |
real dudx(nx,ny,nz) |
real dvdy(nx,ny,nz) |
|
c Calculate horizontal derivatives of pressure |
call deriv(dudx,u,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv(dvdy,v,'y',z3,x2,y2,nx,ny,nz,mdv) |
|
c Calculation |
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
if ((abs(dudx(i,j,k)-mdv).gt.eps).and. |
> (abs(dvdy(i,j,k)-mdv).gt.eps)) then |
|
div_uv(i,j,k)=dudx(i,j,k)+dvdy(i,j,k) |
|
endif |
|
enddo |
enddo |
enddo |
|
end |
|
c ---------------------------------------------------------------- |
c Calculate Q vector |
c ---------------------------------------------------------------- |
|
subroutine calc_qvec (qx3,qy3, |
> th3,u3,v3, |
> z3,x2,y2,f2,nx,ny,nz,mdv) |
|
c Calculate teh Q vector components <qx3> and <qy3>, as well as the divergence |
c of the Q vector <divq3> on the model grid. The grid is specified in the horizontal |
c by <xmin,ymin,dx,dy,nx,ny>. The number of vertical levels is <nz>. The input field |
c are: potential temperature <th3>, horizontal wind <u3> and <v3>, pressure <p3>. |
c The calculation follows the one described in "Weather Analysis, Dusan Djuric" |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real qx3(nx,ny,nz) |
real qy3(nx,ny,nz) |
real th3(nx,ny,nz) |
real u3 (nx,ny,nz) |
real v3 (nx,ny,nz) |
real z3 (nx,ny,nz) |
real x2 (nx,ny) |
real y2 (nx,ny) |
real f2 (nx,ny) |
real mdv |
|
c Physical and numerical parameters |
real scale1,scale2 |
parameter (scale1=1.E10,scale2=1.E14) |
real eps |
parameter (eps=0.01) |
real g |
parameter (g=9.80616) |
real tzero |
parameter (tzero=273.16) |
|
c Auxiliary variables |
real dudx(nx,ny,nz) |
real dudy(nx,ny,nz) |
real dvdx(nx,ny,nz) |
real dvdy(nx,ny,nz) |
real dtdx(nx,ny,nz) |
real dtdy(nx,ny,nz) |
integer i,j,k |
|
c Needed derivatives |
call deriv (dudx, u3,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv) |
call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv (dvdy, v3,'y',z3,x2,y2,nx,ny,nz,mdv) |
call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv) |
|
c Calculate vector components of the Q vector |
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
c Evaluate Q vector formula with missing data check |
if ((abs(dudx(i,j,k)-mdv).gt.eps).and. |
> (abs(dudy(i,j,k)-mdv).gt.eps).and. |
> (abs(dvdx(i,j,k)-mdv).gt.eps).and. |
> (abs(dvdy(i,j,k)-mdv).gt.eps).and. |
> (abs(dtdx(i,j,k)-mdv).gt.eps).and. |
> (abs(dtdy(i,j,k)-mdv).gt.eps)) then |
|
qx3(i,j,k) = -g/tzero * (dudx(i,j,k)*dtdx(i,j,k) |
> +dvdx(i,j,k)*dtdy(i,j,k)) |
qy3(i,j,k) = -g/tzero * (dudy(i,j,k)*dtdx(i,j,k) |
> +dvdy(i,j,k)*dtdy(i,j,k)) |
|
else |
qx3(i,j,k)=mdv |
qy3(i,j,k)=mdv |
endif |
|
enddo |
enddo |
enddo |
|
c Scale the output |
do i=1,nx |
do j=1,ny |
do k=1,nz |
if (abs(qx3(i,j,k)-mdv).gt.eps) then |
qx3(i,j,k)=scale1*qx3(i,j,k) |
endif |
if (abs(qy3(i,j,k)-mdv).gt.eps) then |
qy3(i,j,k)=scale1*qy3(i,j,k) |
endif |
enddo |
enddo |
enddo |
|
end |
|
c ---------------------------------------------------------------- |
c Calculate divergence of wind field |
c ---------------------------------------------------------------- |
|
subroutine calc_div_q (div_q,qx,qy, |
> z3,x2,y2,f2,nx,ny,nz,mdv) |
|
c Calculate the divergence (div_q) of the Q vector field if |
c the components (qx,qy) are given. The grid and the missing data |
c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real div_q(nx,ny,nz) |
real qx (nx,ny,nz) |
real qy (nx,ny,nz) |
real z3 (nx,ny,nz) |
real x2 (nx,ny) |
real y2 (nx,ny) |
real f2 (nx,ny) |
real mdv |
|
c Physical parameters and numerical constants |
real eps |
parameter (eps=0.01) |
|
c Auxiliray variables |
integer i,j,k |
real rho |
real dqxdx(nx,ny,nz) |
real dqydy(nx,ny,nz) |
|
c Calculate horizontal derivatives of pressure |
call deriv(dqxdx,qx,'x',z3,x2,y2,nx,ny,nz,mdv) |
call deriv(dqydy,qy,'y',z3,x2,y2,nx,ny,nz,mdv) |
|
c Calculation |
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
if ((abs(dqxdx(i,j,k)-mdv).gt.eps).and. |
> (abs(dqydy(i,j,k)-mdv).gt.eps)) then |
|
div_q(i,j,k)=dqxdx(i,j,k)+dqydy(i,j,k) |
|
endif |
|
enddo |
enddo |
enddo |
|
end |
|
|
c **************************************************************** |
c * SUBROUTINE SECTION: GRID HANDLING * |
c **************************************************************** |
|
c ----------------------------------------------------------------- |
c Horizontal and vertical derivatives for 3d fields |
c ----------------------------------------------------------------- |
|
subroutine deriv (df,f,direction, |
> z3,x2,y2,nx,ny,nz,mdv) |
|
c Calculate horizontal and vertical derivatives of the 3d field <f>. |
c The direction of the derivative is specified in <direction> |
c 'x','y' : Horizontal derivative in x and y direction |
c 'p','z','t','m' : Vertical derivative (pressure, height, theta, model) |
c The 3d field <z3> specifies the isosurfaces along which the horizontal |
c derivatives are calculated or the levels for the vertical derivatives. |
|
implicit none |
|
c Input and output parameters |
integer nx,ny,nz |
real df (nx,ny,nz) |
real f (nx,ny,nz) |
real z3 (nx,ny,nz) |
real x2 (nx,ny) |
real y2 (nx,ny) |
character direction |
real mdv |
|
c Numerical and physical parameters |
real pi180 |
parameter (pi180=3.141592654/180.) |
real deltay |
parameter (deltay=111.1775E3) |
real zerodiv |
parameter (zerodiv=0.00000001) |
real eps |
parameter (eps=0.01) |
|
c Auxiliary variables |
integer i,j,k |
real vmin,vmax |
real scale,lat |
real vu,vl,vuvl,vlvu |
integer o,u,w,e,n,s |
|
c Vertical derivative |
if ((direction.eq.'z').or. |
> (direction.eq.'th').or. |
> (direction.eq.'p').or. |
> (direction.eq.'m').and. |
> (nz.gt.1)) then |
|
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
o=k+1 |
if (o.gt.nz) o=nz |
u=k-1 |
if (u.lt.1) u=1 |
|
if ((abs(f(i,j,o)-mdv).gt.eps).and. |
> (abs(f(i,j,u)-mdv).gt.eps).and. |
> (abs(f(i,j,k)-mdv).gt.eps)) then |
|
vu = z3(i,j,k)-z3(i,j,o) |
vl = z3(i,j,u)-z3(i,j,k) |
vuvl = vu/(vl+zerodiv) |
vlvu = 1./(vuvl+zerodiv) |
|
df(i,j,k) = 1./(vu+vl) |
> * (vuvl*(f(i,j,u)-f(i,j,k)) |
> + vlvu*(f(i,j,k)-f(i,j,o))) |
|
else |
df(i,j,k) = mdv |
endif |
|
enddo |
enddo |
enddo |
|
c Horizontal derivative in the y direction: 3d |
elseif (direction.eq.'y') then |
|
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
s=j-1 |
if (s.lt.1) s=1 |
n=j+1 |
if (n.gt.ny) n=ny |
|
if ((abs(f(i,n,k)-mdv).gt.eps).and. |
> (abs(f(i,j,k)-mdv).gt.eps).and. |
> (abs(f(i,s,k)-mdv).gt.eps)) then |
|
vu = 1000.*(y2(i,j)-y2(i,n)) |
vl = 1000.*(y2(i,s)-y2(i,j)) |
vuvl = vu/(vl+zerodiv) |
vlvu = 1./(vuvl+zerodiv) |
|
df(i,j,k) = 1./(vu+vl) |
> * (vuvl*(f(i,s,k)-f(i,j,k)) |
> + vlvu*(f(i,j,k)-f(i,n,k))) |
|
else |
df(i,j,k) = mdv |
endif |
|
enddo |
enddo |
enddo |
|
c Horizontal derivative in the x direction: 3d |
elseif (direction.eq.'x') then |
|
do i=1,nx |
do j=1,ny |
do k=1,nz |
|
w=i-1 |
if (w.lt.1) w=1 |
e=i+1 |
if (e.gt.nx) e=nx |
|
if ((abs(f(w,j,k)-mdv).gt.eps).and. |
> (abs(f(i,j,k)-mdv).gt.eps).and. |
> (abs(f(e,j,k)-mdv).gt.eps)) then |
|
vu = 1000.*(x2(i,j)-x2(e,j)) |
vl = 1000.*(x2(w,j)-x2(i,j)) |
vuvl = vu/(vl+zerodiv) |
vlvu = 1./(vuvl+zerodiv) |
|
df(i,j,k) = 1./(vu+vl) |
> * (vuvl*(f(w,j,k)-f(i,j,k)) |
> + vlvu*(f(i,j,k)-f(e,j,k))) |
|
else |
df(i,j,k) = mdv |
endif |
|
enddo |
enddo |
enddo |
|
c Undefined direction for derivative |
else |
|
print*,'Invalid direction of derivative... Stop' |
stop |
|
endif |
|
end |
|
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 Read refernece profile from netcdf |
c -------------------------------------------------------------------------------- |
|
SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref, |
> nx,ny,nz,deltax,deltay,deltaz,coriol,oro, |
> pvsrcfile) |
|
c Read the reference profile from file |
c |
c thetaref : Reference potential temperature (K) |
c pressref : Reference pressure (Pa) |
c rhoref : Reference density (kg/m^3) |
c nsqref : Stratification (s^-1) |
c zref : Reference height (m) |
c nx,nny,nz : Grid dimension in x,y,z direction |
c deltax,deltay,deltaz : Grid spacings used for calculations (m) |
c coriol : Coriolis parameter (s^-1) |
c oro : Height of orography (m) |
c pvsrcfile : Input file |
|
implicit none |
|
c Declaration of subroutine parameters |
integer nx,ny,nz |
real nsqref (0:2*nz) |
real thetaref(0:2*nz) |
real rhoref (0:2*nz) |
real pressref(0:2*nz) |
real zref (0:2*nz) |
real deltax,deltay,deltaz |
real coriol (0:nx,0:ny) |
real oro (0:nx,0:ny) |
character*80 pvsrcfile |
|
c Numerical and physical parameters |
real eps |
parameter (eps=0.01) |
|
c Auxiliary variables |
integer cdfid,stat |
integer vardim(4) |
real misdat |
integer ndimin |
real varmin(4),varmax(4),stag(4) |
integer i,j,k,nf1 |
integer ntimes |
real time(200) |
character*80 vnam(100),varname |
integer nvars |
integer isok,ierr |
real x(0:nx,0:ny),y(0:nx,0:ny) |
real mean,count |
|
c Get grid description from topography |
call cdfopn(pvsrcfile,cdfid,stat) |
if (stat.ne.0) goto 997 |
call getvars(cdfid,nvars,vnam,stat) |
if (stat.ne.0) goto 997 |
isok=0 |
varname='ORO' |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdef(cdfid,varname,ndimin,misdat,vardim, |
> varmin,varmax,stag,stat) |
if (stat.ne.0) goto 997 |
time(1)=0. |
call gettimes(cdfid,time,ntimes,stat) |
if (stat.ne.0) goto 997 |
call clscdf(cdfid,stat) |
if (stat.ne.0) goto 997 |
|
c Open output netcdf file |
call cdfopn(pvsrcfile,cdfid,stat) |
if (stat.ne.0) goto 997 |
|
c Create the variable if necessary |
call getvars(cdfid,nvars,vnam,stat) |
if (stat.ne.0) goto 997 |
|
c Read data from netcdf file |
isok=0 |
varname='NSQREF' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,nsqref,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='RHOREF' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,rhoref,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='THETAREF' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,thetaref,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='PREREF' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,pressref,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='ZREF' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,zref,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='CORIOL' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,coriol,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='ORO' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,oro,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='X' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,x,stat) |
if (stat.ne.0) goto 997 |
|
isok=0 |
varname='Y' |
print*,'R ',trim(varname),' ',trim(pvsrcfile) |
call check_varok(isok,varname,vnam,nvars) |
if (isok.eq.0) goto 997 |
call getdat(cdfid,varname,time(1),0,y,stat) |
if (stat.ne.0) goto 997 |
|
c Close netcdf file |
call clscdf(cdfid,stat) |
if (stat.ne.0) goto 997 |
|
c Determine the grid spacings <deltax, deltay, deltaz> |
mean=0. |
count=0. |
do i=1,nx |
do j=0,ny |
mean=mean+abs(x(i)-x(i-1)) |
count=count+1. |
enddo |
enddo |
deltax=mean/count |
|
mean=0. |
count=0. |
do j=1,ny |
do i=0,nx |
mean=mean+abs(y(j)-y(j-1)) |
count=count+1. |
enddo |
enddo |
deltay=mean/count |
|
mean=0. |
count=0. |
do k=1,nz-1 |
mean=mean+abs(zref(k+1)-zref(k-1)) |
count=count+1. |
enddo |
deltaz=mean/count |
|
return |
|
c Exception handling |
997 print*,'Read_Ref: Problem with input netcdf file... Stop' |
stop |
|
end |
Property changes: |
Added: svn:executable |