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(:),ut(:),dday(:) integer,allocatable,save :: doy(:) 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 if (allocated(doy)) deallocate(doy) allocate(doy(ntime),stat=istat) if (istat /= 0) then write(6,"('>>> init_data: error allocating doy: istat=', | i5)") istat stop 'allocate doy' endif if (allocated(ut)) deallocate(ut) allocate(ut(ntime),stat=istat) if (istat /= 0) then write(6,"('>>> init_data: error allocating ut: istat=', | i5)") istat stop 'allocate ut' endif if (allocated(dday)) deallocate(dday) allocate(dday(ntime),stat=istat) if (istat /= 0) then write(6,"('>>> init_data: error allocating dday: istat=', | i5)") istat stop 'allocate dday' endif 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_doy ,id_ut ,id_dday integer :: id_fu,id_fv ! var ids integer :: ids1(1), ids4(4) integer :: start_4d(4), count_4d(4) real :: time_local(ntime) character(len=120) :: char120 ! ! Write 4-byte reals to the nc file: real*4 :: glb_4byte(nlon,nlat,nlev), | lon_4byte(nlon), lat_4byte(nlat), lev_4byte(nlev), | time_4byte(ntime),ut_4byte(ntime),dday_4byte(ntime) ! ! 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,"lon",nlon,id_nlon) istat = nf_def_dim(ncid,"lat",nlat,id_nlat) istat = nf_def_dim(ncid,"lev",nlev,id_nlev) ! ! Define coord variables: ids1(1) = id_ntime istat = nf_def_var(ncid,"day",NF_INT,1,ids1,id_doy) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var time',1) write(char120,"('Day of year (1-365)')") istat = nf_put_att_text(ncid,id_doy,"long_name", | len_trim(char120),trim(char120)) ids1(1) = id_ntime istat = nf_def_var(ncid,"ut",NF_FLOAT,1,ids1,id_ut) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var ut',1) write(char120,"('Universal Time')") istat = nf_put_att_text(ncid,id_ut,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_ut,"units",5,"hours") ids1(1) = id_ntime istat = nf_def_var(ncid,"dday",NF_FLOAT,1,ids1,id_dday) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var dday',1) write(char120,"('Fractional day of year')") istat = nf_put_att_text(ncid,id_dday,"long_name", | len_trim(char120),trim(char120)) ids1(1) = id_ntime istat = nf_def_var(ncid,"time",NF_FLOAT,1,ids1,id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var time',1) write(char120,"('UT since beginning of October 1')") istat = nf_put_att_text(ncid,id_time,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_time,"units",7,"minutes") ids1(1) = id_nlon istat = nf_def_var(ncid,"lon",NF_FLOAT,1,ids1,id_lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lon',1) write(char120,"('geographic longitude (-west, +east)')") istat = nf_put_att_text(ncid,id_lon,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_lon,"units",12,'degrees_east') ids1(1) = id_nlat istat = nf_def_var(ncid,"lat",NF_FLOAT,1,ids1,id_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lat',1) write(char120,"('geographic latitude (-south, +north)')") istat = nf_put_att_text(ncid,id_lat,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_lat,"units",13,'degrees_north') ids1(1) = id_nlev istat = nf_def_var(ncid,"lev",NF_FLOAT,1,ids1,id_lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining coord var lev',1) istat = nf_put_att_text(ncid,id_lev,"long_name",8,"altitude") istat = nf_put_att_text(ncid,id_lev,"units",2,'km') ! ! 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) write(char120,"('Forcing in zonal direction')") istat = nf_put_att_text(ncid,id_fu,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_fu,"units",5,'m/s^2') 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) write(char120,"('Forcing in meridional direction')") istat = nf_put_att_text(ncid,id_fv,"long_name", | len_trim(char120),trim(char120)) istat = nf_put_att_text(ncid,id_fv,"units",5,'m/s^2') ! ! Take out of define mode: istat = nf_enddef(ncid) ! ! Write coord variables to the file: lon_4byte = lon istat = nf_put_var_real(ncid,id_lon,lon_4byte) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing coord var lon',1) lat_4byte = lat istat = nf_put_var_real(ncid,id_lat,lat_4byte) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing coord var lat',1) lev_4byte = lev istat = nf_put_var_real(ncid,id_lev,lev_4byte) 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 time_4byte = time doy(:) = 274 ! Oct 1 do itime=1,ntime istat = nf_put_var1_real(ncid,id_time,itime,time_4byte(itime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var time',1) if (time(itime) <= 1440.) then ut_4byte(itime) = time(itime)/60. dday_4byte(itime) = float(doy(itime))+ut_4byte(itime)/24. else ut_4byte(itime) = (time(itime)-1440.)/60. doy(itime) = doy(1)+1 dday_4byte(itime) = float(doy(itime))+ut_4byte(itime)/24. endif istat = nf_put_var1_int(ncid,id_doy,itime,doy(itime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var doy',1) istat = nf_put_var1_real(ncid,id_ut,itime,ut_4byte(itime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var ut',1) istat = nf_put_var1_real(ncid,id_dday,itime,dday_4byte(itime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var dday',1) start_4d(4) = itime write(6,"('Write fu,fv: itime=',i4,' time=',f8.2,' fu min,max=', | 2e12.4,' fv min,max=',2e12.4)") itime,time(itime), | minval(fu(:,:,:,itime)),maxval(fu(:,:,:,itime)), | minval(fv(:,:,:,itime)),maxval(fv(:,:,:,itime)) glb_4byte = fu(:,:,:,itime) istat = nf_put_vara_real(ncid,id_fu,start_4d,count_4d,glb_4byte) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error writing data var FU',1) glb_4byte = fv(:,:,:,itime) istat = nf_put_vara_real(ncid,id_fv,start_4d,count_4d,glb_4byte) 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,it ! ! 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,"(/,'nc_read: Opened file ',a,' for reading.')") | trim(filename) istat = nf_inq_dimid(ncid,"lon",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,"lat",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,"lev",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) if (allocated(time)) deallocate(time) allocate(time(ntime),stat=istat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'nc_read: Error allocating time',1) if (allocated(doy)) deallocate(doy) ! int day of year (ntime) allocate(doy(ntime),stat=istat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'nc_read: Error allocating doy',1) if (allocated(ut)) deallocate(ut) ! ut (ntime) allocate(ut(ntime),stat=istat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'nc_read: Error allocating ut',1) if (allocated(dday)) deallocate(dday) ! decimal day allocate(dday(ntime),stat=istat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'nc_read: Error allocating dday',1) istat = nf_inq_varid(ncid,"time",id) istat = nf_get_var_double(ncid,id,time) istat = nf_inq_varid(ncid,"day",id) istat = nf_get_var_int(ncid,id,doy) istat = nf_inq_varid(ncid,"ut",id) istat = nf_get_var_double(ncid,id,ut) istat = nf_inq_varid(ncid,"dday",id) istat = nf_get_var_double(ncid,id,dday) write(6,"('nc_read: ntime=',i4,' time=',/,(6e12.4))") ntime,time write(6,"('nc_read: ntime=',i4,' doy =',/,(6i6))") ntime,doy write(6,"('nc_read: ntime=',i4,' ut =',/,(6f9.4))") ntime,ut write(6,"('nc_read: ntime=',i4,' dday=',/,(6f9.4))") ntime,dday 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