Rev 7 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
c ************************************************************
c * This package provides input routines to read the wind *
c * and other fields from IVE necdf files. The routines are *
c * *
c * 1) input_open : to open a data file *
c * 2) input_grid : to read the grid information, including *
c * the vertical levels *
c * 3) input_wind : to read the wind components *
c * 4) input_close : to close an input file *
c * *
c * The file is characterised by an filename <filename> and *
c * a file identifier <fid>. The horizontal grid is given by *
c * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the *
c * rotated grid is given by <pollon,pollat>. The vertical *
c * grid is characterised by the surface pressure <ps> and *
c * the pressure at staggered <slev> and unstaggered <ulev> *
c * levels. The number of levels is given by <nz>. Finally, *
c * the retrieval of the wind <field> with name <fieldname> *
c * is characterised by a <time> and a missing data value *
c * <mdv>. *
c * *
c * Author: Michael Sprenger, Autumn 2008 *
c ************************************************************
c ------------------------------------------------------------
c Open input file
c ------------------------------------------------------------
subroutine input_open (fid,filename)
c Open the input file with filename <filename> and return the
c file identifier <fid> for further reference.
use netcdf
implicit none
c Declaration of subroutine parameters
integer fid ! File identifier
character*80 filename ! Filename
c Declaration of auxiliary variables
integer ierr
c Open netcdf file
ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
c Exception handling
return
900 print*,'Cannot open input file ',trim(filename)
stop
end
c ------------------------------------------------------------
c Read information about the grid
c ------------------------------------------------------------
subroutine input_grid
> (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
> time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
> timecheck)
c Read grid information at <time> from file with identifier <fid>.
c The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
c with pole position at <pollon,pollat> and grid dimension <nx,ny>.
c The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
c on the staggered or unstaggered grid (with <stagz> as the flag).
c The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
c only the grid dimensions and grid parameters (xmin...pollat,nz) are
c determined and returned (this is needed for dynamical allocation of
c memory).
use netcdf
implicit none
c Declaration of subroutine parameters
integer fid ! File identifier
real xmin,xmax,ymin,ymax ! Domain size
real dx,dy ! Horizontal resolution
integer nx,ny,nz ! Grid dimensions
real pollon,pollat ! Longitude and latitude of pole
real p3(nx,ny,nz) ! Staggered levels
real ps(nx,ny) ! Surface pressure
real time ! Time of the grid information
real ak(nz),bk(nz) ! Ak and Bk for layers or levels
real stagz ! Vertical staggering (0 or -0.5)
character*80 fieldname ! Variable from which to take grid info
character*80 timecheck ! Either 'yes' or 'no'
c Numerical and physical parameters
real eps ! Numerical epsilon
parameter (eps=0.001)
c Netcdf variables
integer vardim(4)
real varmin(4),varmax(4)
real mdv
real stag(4)
integer ndim
character*80 cstfile
integer cstid
real times
integer ntimes
real aklay(nz),bklay(nz),aklev(nz),bklev(nz)
integer nvars
character*80 vars(100)
integer dimids (nf90_max_var_dims)
character*80 dimname(nf90_max_var_dims)
c Auxiliary varaibles
integer ierr
integer i,j,k
integer isok
real tmp(200)
character*80 varname
real rtime
integer varid
integer cdfid
real rmax,rmin
c Set file identifier
if (fid.lt.0) then
cdfid = -fid
else
cdfid = fid
endif
c Get varid
ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get numebr of dimensions -> ndim
ierr = nf90_inquire_variable(cdfid, varid, ndims = ndim)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get dim IDs
ierr = nf90_inquire_variable(cdfid, varid,
> dimids = dimids(1:ndim))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get dimensions -> vardim(1:ndim)
do i=1,ndim
ierr = nf90_inquire_dimension(cdfid, dimids(i),
> name = dimname(i) )
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
enddo
c Get domain min -> varmin(1:3)
ierr = nf90_get_att(cdfid, varid, "xmin", varmin(1))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_att(cdfid, varid, "ymin", varmin(2))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_att(cdfid, varid, "zmin", varmin(3))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get domain max -> varmax(1:3)
ierr = nf90_get_att(cdfid, varid, "xmax", varmax(1))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_att(cdfid, varid, "ymax", varmax(2))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_att(cdfid, varid, "zmax", varmax(3))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get vertical staggering -> stagz
ierr = nf90_get_att(cdfid, varid, "zstag", stagz)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get missing data value -> mdv
ierr = nf90_get_att(cdfid, varid, "missing_value", mdv)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Set the grid dimensions and constants
nx = vardim(1)
ny = vardim(2)
nz = vardim(3)
xmin = varmin(1)
ymin = varmin(2)
xmax = varmax(1)
ymax = varmax(2)
dx = (xmax-xmin)/real(nx-1)
dy = (ymax-ymin)/real(ny-1)
c We want the longitudes within -180 ... + 180
if ( xmin.lt.-180.) then
xmin = xmin + 360.
xmax = xmax + 360.
else if ( xmax.gt.360. ) then
xmin = xmin - 360.
xmax = xmax - 360.
endif
c Get name of constants file -> cstfile
ierr = nf90_get_att(cdfid,nf90_global,
> "constants_file_name ", cstfile )
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get pole position from constants file
ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
varname='pollon'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,pollon)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
varname='pollat'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,pollat)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = NF90_CLOSE(cstid)
IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Skip the rest if fid < 0:
if ( fid.lt.0 ) goto 900
c Get ak and bk from constants file
ierr = NF90_OPEN(TRIM(cstfile),nf90_nowrite, cstid)
IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
varname='aklev'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,aklev)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
varname='bklev'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,bklev)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
varname='aklay'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,aklay)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
varname='bklay'
ierr = NF90_INQ_VARID(cstid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cstid,varid,bklay)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = NF90_CLOSE(cstid)
IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c print*,'AKLEV : '
c print*, (aklev(i),i=1,20)
c print*,'BKLEV : '
c print*,(bklev(i),i=1,20)
c print*,'AKLAY : '
c print*,(aklay(i),i=1,20)
c print*,'BKLAY : '
c print*,(bklay(i),i=1,20)
c Get time information (check if time is correct)
varname = 'time'
ierr = NF90_INQ_VARID(cdfid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cdfid,varid,times)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
isok=0
if (abs(time-times).lt.eps) then
isok = 1
rtime = times
elseif (timecheck.eq.'no') then
isok = 1
rtime = times
endif
if ( isok.eq.0 ) then
print*,' ERROR: time ',rtime,' not found on netCDF file'
stop
endif
c print*,'TIME : ',rtime
c Read surface pressure
varname='PS'
ierr = NF90_INQ_VARID(cdfid,varname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_get_var(cdfid,varid,ps)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Check that PS is in hPa - quick and dirty
rmax = -1.e9
rmin = +1.e9
do i=1,nx
do j=1,ny
if ( (ps(i,j).gt.rmax).and.abs(ps(i,j)-mdv).gt.eps) then
rmax = ps(i,j)
endif
if ( (ps(i,j).lt.rmin).and.abs(ps(i,j)-mdv).gt.eps) then
rmin = ps(i,j)
endif
enddo
enddo
if ( (rmin.lt.200.).or.(rmax.gt.1500.) ) then
print*,' ERROR: PS must be in hPa... Stop'
stop
endif
c Calculate layer and level pressures
do i=1,nx
do j=1,ny
do k=1,nz
if ( abs(stagz).lt.eps ) then
p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
else
p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
endif
enddo
enddo
enddo
c Set the ak and bk for the vertical grid
do k=1,nz
if ( abs(stagz).lt.eps ) then
ak(k)=aklev(k)
bk(k)=bklev(k)
else
ak(k)=aklay(k)
bk(k)=bklay(k)
endif
enddo
c Exit point
900 continue
return
end
c ------------------------------------------------------------
c Read wind information
c ------------------------------------------------------------
subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
> timecheck)
c Read the wind component <fieldname> from the file with identifier
c <fid> and save it in the 3d array <field>. The vertical staggering
c information is provided in <stagz> and gives the reference to either
c the layer or level field from <input_grid>. A consistency check is
c performed to have an agreement with the grid specified by <xmin,xmax,
c ymin,ymax,dx,dy,nx,ny,nz>.
use netcdf
implicit none
c Declaration of variables and parameters
integer fid ! File identifier
character*80 fieldname ! Name of the wind field
integer nx,ny,nz ! Dimension of fields
real field(nx,ny,nz) ! 3d wind field
real stagz ! Staggering in the z direction
real mdv ! Missing data flag
real xmin,xmax,ymin,ymax ! Domain size
real dx,dy ! Horizontal resolution
real time ! Time
character*80 timecheck ! Either 'yes' or 'no'
c Numerical and physical parameters
real eps ! Numerical epsilon
parameter (eps=0.001)
real notimecheck ! 'Flag' for no time check
parameter (notimecheck=7.26537)
c Netcdf variables
integer ierr
character*80 varname
integer vardim(4)
real varmin(4),varmax(4)
real stag(4)
integer ndim
real times(10)
integer ntimes
character*80 cstfile
integer cstid
real aklay(200),bklay(200),aklev(200),bklev(200)
real ps(nx,ny)
integer dimids (nf90_max_var_dims)
character*80 dimname(nf90_max_var_dims)
integer varid
integer cdfid
c Auxiliary variables
integer isok
integer i,j,k
integer nz1
real rtime
c Get varid
ierr = NF90_INQ_VARID(fid,fieldname,varid)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c Get number of vertical levels -> vardim(3)
ierr = nf90_inquire_variable(fid, varid, ndims = ndim)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_inquire_variable(fid, varid,
> dimids = dimids(1:ndim))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
do i=1,ndim
ierr = nf90_inquire_dimension(fid, dimids(i),
> name = dimname(i) )
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
enddo
c Read data
varname=fieldname
ierr = nf90_get_var(fid,varid,field)
IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
c If the field is 2D, expand it to 3D - simple handling of 2D tracing
if ( vardim(3).eq.1 ) then
do i=1,nx
do j=1,ny
do k=1,nz
field(i,j,k) = field(i,j,1)
enddo
enddo
enddo
endif
c Exit point
return
end
c ------------------------------------------------------------
c Close input file
c ------------------------------------------------------------
subroutine input_close(fid)
c Close the input file with file identifier <fid>.
use netcdf
implicit none
c Declaration of subroutine parameters
integer fid
c Auxiliary variables
integer ierr
c Close file
ierr = NF90_CLOSE(fid)
IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
end