module nc_gwforcing ! ! B. Foster: March, 2009 ! Perform i/o on netcdf files containing Fx,Fy forcing from ! Sharon Vadas' ray trace program. ! implicit none #include save ! ! Spatial dimensions: integer,parameter :: | nlon=111, ! number of longitudes (nx) | nlat= 55, ! number of latitudes (ny) | nlev= 51 ! number of levels (nz) ! ! Spatial coordinate arrays: real :: lon(nlon), lat(nlat), lev(nlev) ! ! Forcing fu,fv (nlon,nlat,nlev,ntime) real,allocatable,save :: | fu(:,:,:,:), ! Force in zonal direction (fx) | fv(:,:,:,:) ! Force in meridional direction (fy) real,allocatable,save :: time(:) contains !----------------------------------------------------------------------- subroutine init_data(ntime) implicit none ! ! Allocate and initialize data for ntime times. ! ! Args: integer,intent(in) :: ntime ! number of times ! ! Local: integer :: istat real,parameter :: init_val=0. ! if (ntime < 1) then write(6,"('>>> init_data: ntime=',i5,' but must be >= 1')") | ntime stop 'init_data' endif if (allocated(fu)) deallocate(fu) allocate(fu(nlon,nlat,nlev,ntime),stat=istat) if (istat /= 0) then write(6,"('>>> init_data: error allocating fu: istat=', | i5)") istat stop 'allocate fu' endif if (allocated(fv)) deallocate(fv) allocate(fv(nlon,nlat,nlev,ntime),stat=istat) if (istat /= 0) then write(6,"('>>> init_data: error allocating fv: istat=', | i5)") istat stop 'allocate fv' endif fu = init_val fv = init_val write(6,"('init_data: allocated fu,fv. Init_val=',e12.4)") | init_val end subroutine init_data !----------------------------------------------------------------------- subroutine nc_write(filename,ntime) implicit none ! ! Create, define, and load data to netcdf file filename. ! Assumes that fu,fv have been allocated, and contain data for ntime times. ! ! Args: character(len=*),intent(in) :: filename integer,intent(in) :: ntime ! number of times in fu,fv ! ! Local: integer :: ncid,istat,itime integer :: id_nlon,id_nlat,id_nlev,id_ntime ! dimension ids integer :: id_lon ,id_lat ,id_lev ,id_time ! coord var ids integer :: id_fu,id_fv ! var ids integer :: ids1(1), ids4(4) integer :: start_4d(4), count_4d(4) real :: time_local(ntime) real*4 :: real4(nlon,nlat,nlev) ! ! Create new netcdf dataset: istat = nf_create(filename,NF_CLOBBER,ncid) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error from nf_create',1) write(6,"(/,'nc_write: Created netcdf file ',a)") trim(filename) ! ! Define dimensions: istat = nf_def_dim(ncid,"time",NF_UNLIMITED,id_ntime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining time dimension',1) istat = nf_def_dim(ncid,"nlon",nlon,id_nlon) istat = nf_def_dim(ncid,"nlat",nlat,id_nlat) istat = nf_def_dim(ncid,"nlev",nlev,id_nlev) ! ! Define coord variables: ids1(1) = id_ntime istat = nf_def_var(ncid,"time",NF_DOUBLE,1,ids1,id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var time',1) ids1(1) = id_nlon istat = nf_def_var(ncid,"lon",NF_DOUBLE,1,ids1,id_lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lon',1) ids1(1) = id_nlat istat = nf_def_var(ncid,"lat",NF_DOUBLE,1,ids1,id_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lat',1) ids1(1) = id_nlev istat = nf_def_var(ncid,"lev",NF_DOUBLE,1,ids1,id_lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lev',1) ! ! Define 4d data variables fu,fv: ids4(1) = id_nlon ids4(2) = id_nlat ids4(3) = id_nlev ids4(4) = id_ntime istat = nf_def_var(ncid,"FU",NF_FLOAT,4,ids4,id_fu) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining data var FU',1) istat = nf_def_var(ncid,"FV",NF_FLOAT,4,ids4,id_fv) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining data var FV',1) ! ! Take out of define mode: istat = nf_enddef(ncid) ! ! Write coord variables to the file: istat = nf_put_var_double(ncid,id_lon,lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing coord var lon',1) istat = nf_put_var_double(ncid,id_lat,lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing coord var lat',1) istat = nf_put_var_double(ncid,id_lev,lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing coord var lev',1) ! ! Write data variables to the file (write fu,fv as 4-byte reals). start_4d(1:3) = 1 count_4d(1) = nlon count_4d(2) = nlat count_4d(3) = nlev count_4d(4) = 1 do itime=1,ntime istat = nf_put_var1_double(ncid,id_time,itime,time(itime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var time',1) start_4d(4) = itime write(6,"('Write fu,fv: itime=',i4,' fu min,max=',2e12.4, | ' fv min,max=',2e12.4)") itime, | minval(fu(:,:,:,itime)),maxval(fu(:,:,:,itime)), | minval(fv(:,:,:,itime)),maxval(fv(:,:,:,itime)) real4 = fu(:,:,:,itime) istat = nf_put_vara_real(ncid,id_fu,start_4d,count_4d,real4) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var FU',1) real4 = fv(:,:,:,itime) istat = nf_put_vara_real(ncid,id_fv,start_4d,count_4d,real4) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var FV',1) enddo istat = nf_close(ncid) end subroutine nc_write !----------------------------------------------------------------------- subroutine nc_read(filename) implicit none ! ! Read netcdf force file. ! Read fu,fv into 8-byte reals (for timegcm), even tho they were written ! as 4-byte reals. ! ! Args: character(len=*),intent(in) :: filename ! ! Local: integer :: istat,ncid,id,ntime ! ! Open file for reading: istat = nf_open(filename,NF_NOWRITE,ncid) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error from nf_open',1) write(6,"(/,'read_hist: Opened file ',a,' for reading.')") | trim(filename) istat = nf_inq_dimid(ncid,"nlon",id) istat = nf_inq_dimlen(ncid,id,nlon) istat = nf_inq_varid(ncid,"lon",id) istat = nf_get_var_double(ncid,id,lon) write(6,"('nc_read: nlon=',i4,' lon=',/,(6e12.4))") nlon,lon istat = nf_inq_dimid(ncid,"nlat",id) istat = nf_inq_dimlen(ncid,id,nlat) istat = nf_inq_varid(ncid,"lat",id) istat = nf_get_var_double(ncid,id,lat) write(6,"('nc_read: nlat=',i4,' lat=',/,(6e12.4))") nlat,lat istat = nf_inq_dimid(ncid,"nlev",id) istat = nf_inq_dimlen(ncid,id,nlev) istat = nf_inq_varid(ncid,"lev",id) istat = nf_get_var_double(ncid,id,lev) write(6,"('nc_read: nlev=',i4,' lev=',/,(6e12.4))") nlev,lev istat = nf_inq_dimid(ncid,"time",id) istat = nf_inq_dimlen(ncid,id,ntime) istat = nf_inq_varid(ncid,"time",id) istat = nf_get_var_double(ncid,id,time) write(6,"('nc_read: ntime=',i4,' time=',/,(6e12.4))") ntime,time istat = nf_inq_varid(ncid,"FU",id) istat = nf_get_var_double(ncid,id,fu) write(6,"('nc_read: fu min,max=',2e12.4)") minval(fu),maxval(fu) istat = nf_inq_varid(ncid,"FV",id) istat = nf_get_var_double(ncid,id,fv) write(6,"('nc_read: fv min,max=',2e12.4)") minval(fv),maxval(fv) istat = nf_close(ncid) end subroutine nc_read !----------------------------------------------------------------------- subroutine read_text(filename,ntime) ! ! Read dimensions and coordinate arrays from text file filename. ! Check dimensions read from the file against module dimension parameters. ! Define module coordinate arrays, and return number of times in ntime. ! ! Args: character(len=*),intent(in) :: filename ! input text file integer,intent(out) :: ntime ! ! Local: integer :: iunit=7,nx,ny,nz,istat ! ! Open text file, read dimensions, and confirm they match parameters: open(unit=iunit,file=filename,status='old') read(iunit,*) nx,ny,nz,ntime write(6,"('read_text: nx=',i4,' ny=',i4,' nz=',i4,' ntime=',i4)") | nx,ny,nz,ntime if (nx /= nlon .or. ny /= nlat .or. nz /= nlev) then write(6,"('>>> read_text: dimensions read from file ',a, | ' do not match parameter dimensions:')") write(6,"('Read dims: nx =',i4,' ny =',i4,' nz =',i4)") | nx,ny,nz write(6,"('Param dims: nlon=',i4,' nlat=',i4,' nlev=',i4)") | nlon,nlat,nlev stop 'read_text' endif ! ! Allocate and init time coordinate array (module data): allocate(time(ntime),stat=istat) if (istat /= 0) then write(6,"('>>> read_text: error allocating time: ntime=',i4, | ' istat=',i5)") ntime,istat stop 'allocate time' endif time = 0. write(6,"('read_text: allocated time(ntime=',i5,')')") ntime ! ! Read coord arrays: read(iunit,*) lon write(6,"('read_text: lon=',/,(6e12.4))") lon read(iunit,*) lat write(6,"('read_text: lat=',/,(6e12.4))") lat read(iunit,*) lev write(6,"('read_text: lev=',/,(6e12.4))") lev read(iunit,*) time write(6,"('read_text: time=',/,(6e12.4))") time end subroutine read_text !----------------------------------------------------------------------- subroutine read_binary(filename,ntime) implicit none ! ! Read fortran binary force file containing data for ntime times. ! Read 4-byte reals from binary file, and transfer to 8-byte module data ! ! Args: character(len=*),intent(in) :: filename integer,intent(in) :: ntime ! ! Local: integer,parameter :: iword=4 integer :: iunit=7,it,k,j,irecnum real*4,dimension(nlon,nlat,nlev,ntime) :: fx,fy open(unit=iunit,access='direct',form='unformatted', | file=filename,status='old',recl=nlon*iword) do it=1,ntime do k=1,nlev do j=1,nlat irecnum = j+(k-1)*nlat + (it-1)*nlat*nlev read(iunit,rec=irecnum) fx(:,j,k,it) irecnum = irecnum + nlat*nlev*ntime read(iunit,rec=irecnum) fy(:,j,k,it) enddo enddo enddo close(iunit) ! ! Transfer to 8-byte module data: fu = fx fv = fy write(6,"('read_binary: ntime=',i5,' fu min,max=',2e12.4)") | ntime, minval(fu),maxval(fu) write(6,"('read_binary: ntime=',i5,' fv min,max=',2e12.4)") | ntime, minval(fv),maxval(fv) end subroutine read_binary !----------------------------------------------------------------------- subroutine handle_ncerr(istat,msg,ifatal) ! ! Handle a netcdf lib error: ! integer,intent(in) :: istat,ifatal character(len=*),intent(in) :: msg ! write(6,"(/72('-'))") write(6,"('>>> Error from netcdf library:')") write(6,"(a)") trim(msg) write(6,"('istat=',i5)") istat write(6,"(a)") nf_strerror(istat) write(6,"(72('-')/)") if (ifatal > 0) stop 'handle_ncerr' return end subroutine handle_ncerr !----------------------------------------------------------------------- end module nc_gwforcing