! +-======-+ ! Copyright (c) 2003-2007 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ program dyndiff !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS ! !----------------------------------------------------------------------- !BOP ! ! !ROUTINE: dyndiff: find the dynamic vector difference between two hdf files ! ! !USAGE: see the routine usage() below ! ! !USES: ! use m_dyn implicit NONE ! !DESCRIPTION: Uses statistic subroutine dyn_stat to check the dynamic ! vector dfference between two alike hdf files. ! ! !REVISION HISTORY: ! ! 2002.01.07 E. Yeh: Initial code. ! 2005.03.28 Elena N.: Added an option to save resulting dynamic vector in a file. ! 01May2007 Todling Support to new GEOS-5 dyn-vector files ! 06May2007 Todling Add opt to print ave/min/max when one file is given as input. ! 05Mar2009 Todling Add land fractions ! !------------------------------------------------------------------------- !EOP character(len=*), parameter :: myname = 'dyndiff' ! File names ! ---------- integer, parameter :: MFILES = 2 ! max. number of input files character(len=255) :: files(MFILES), binfiles(MFILES), etafile, st, hroot integer :: nfiles ! actual no. of input files character(len=255) :: dyn_dout ! difference output file name ! Dynamics/simulator vectors ! -------------------------- type(dyn_vect) dyn(MFILES) ! dynamics vector in eta ! Locals ! ------ integer, parameter :: READ_ONLY = 1 integer fid, nvars, ngatts integer ios, rc, iopt, ifile integer ntimes, n, freq, nymd, nhms, prec integer freq_d, nymd_d, nhms_d, prec_d !Timetag for newly created diff in *.hdf format integer im, jm, km, system, dyntype logical dominmax ! ******* ! files(1) = "/scratch1/eyeh/TT13b4/OUTPUT/ana/Y1997/M12/TT13b4.bkg.eta.19971215.hdf" ! files(2) = "/scratch1/eyeh/TD13b4/OUTPUT/ana/Y1997/M12/TD13b4.bkg.eta.19971215.hdf" ! files(1) = "/scratch1/eyeh/NG/bfv125_01_y99.bkg.eta.20000730.hdf" ! files(2) = "/scratch1/eyeh/NG/bfv_7_30.bkg.eta.20000730.hdf" ! files(1) = "/scratch1/eyeh/TI13b4/OUTPUT/ana/Y1997/M12/TI13b4.bkg.eta.19971215.hdf" ! files(2) = "/scratch1/eyeh/TD13b6/OUTPUT/ana/Y1997/M12/TD13b6.bkg.eta.19971215.hdf" ! Initialize ! ---------- call Init_ ( dyntype, mfiles, files, dominmax ) call getenv("HOME", hroot) nfiles = 1 ! Loop over input eta files ! ------------------------- do ifile = 1, nfiles etafile = files(ifile) ! Determine how many time levels on file ! -------------------------------------- call GFIO_Open ( etafile, READ_ONLY, fid, rc ) if ( rc .ne. 0 ) then call die(myname,'cannot open GFIO file '//trim(etafile)) end if call GFIO_DimInquire ( fid, im, jm, km, ntimes, nvars, ngatts, rc) if ( rc .ne. 0 ) then call die(myname,'problems getting dimensions' ) end if call GFIO_Close ( fid, rc ) ! For each time on file... ! ------------------------ do n = 1, ntimes ! Get ETA data for this time ! -------------------------- call dyn_get ( etafile, nymd, nhms, dyn(1), rc, timidx=n, freq=freq, vectype=dyntype ) nymd_d = nymd !Newly created diff file with nymd from first file dyn(1) nhms_d = nhms !Newly created diff file with nhms from first file dyn(1) if ( rc .ne. 0 ) then call die(myname,'cannot read dynamics vector file') end if ! print *, "nymd, nhms:", nymd, nhms ! print *, "timidx, freq:", n, freq call dyn_get ( files(2), nymd, nhms, dyn(2), rc, timidx=n, freq=freq, vectype=dyntype ) if ( rc .ne. 0 ) then call die(myname,'cannot read dynamics vector file') end if !bin_diff ! binfiles(1) = trim(files(1)) // "_bin" ! binfiles(2) = trim(files(2)) // "_bin" !_RT binfiles(1) = trim(hroot) // "/.#_dyndiff-A.bin-01" !_RT binfiles(2) = trim(hroot) // "/.#_dyndiff-A.bin-02" !_RT open(unit=21, file=trim(binfiles(1)), form="unformatted", action="write", iostat=rc) !_RT if(rc .ne. 0) call die(myname,'cannot open output file unit 21') !_RT open(unit=22, file=trim(binfiles(2)), form="unformatted", action="write", iostat=rc) !_RT if(rc .ne. 0) call die(myname,'cannot open output file unit 22') !_RT write (21) dyn(1)%ps,dyn(1)%ts,dyn(1)%phis,dyn(1)%lwi,dyn(1)%hs_stdv, & !_RT dyn(1)%frland,dyn(1)%frlandice,dyn(1)%frlake,dyn(1)%frocean,dyn(1)%frseaice,& !_RT dyn(1)%delp,dyn(1)%u,dyn(1)%v,dyn(1)%pt,dyn(1)%q !_RT write (22) dyn(2)%ps,dyn(2)%ts,dyn(2)%phis,dyn(2)%lwi,dyn(2)%hs_stdv, & !_RT dyn(2)%frland,dyn(2)%frlandice,dyn(2)%frlake,dyn(2)%frocean,dyn(2)%frseaice,& !_RT dyn(2)%delp,dyn(2)%u,dyn(2)%v,dyn(2)%pt,dyn(2)%q !_RT st = "diff " // trim(binfiles(1)) // " " // trim(binfiles(2)) !_RT close(21) !_RT close(22) rc = system(st) ! rc == 0 for zero diff !_RT if (rc .eq. 0 .and. trim(files(1)).ne.trim(files(2))) then !_RT print *, "> nymd, nhms: ", nymd, nhms, " (zero diff)" !_RT else print *, "> nymd, nhms: ", nymd, nhms, " (diff)" if ( .not. dominmax ) then dyn(1)%ps = dyn(1)%ps - dyn(2)%ps dyn(1)%ts = dyn(1)%ts - dyn(2)%ts dyn(1)%phis = dyn(1)%phis - dyn(2)%phis dyn(1)%lwi = dyn(1)%lwi - dyn(2)%lwi dyn(1)%frland = dyn(1)%frland - dyn(2)%frland dyn(1)%frlandice = dyn(1)%frlandice - dyn(2)%frlandice dyn(1)%frlake = dyn(1)%frlake - dyn(2)%frlake dyn(1)%frocean = dyn(1)%frocean - dyn(2)%frocean dyn(1)%frseaice = dyn(1)%frseaice - dyn(2)%frseaice dyn(1)%hs_stdv = dyn(1)%hs_stdv - dyn(2)%hs_stdv dyn(1)%delp = dyn(1)%delp - dyn(2)%delp dyn(1)%u = dyn(1)%u - dyn(2)%u dyn(1)%v = dyn(1)%v - dyn(2)%v dyn(1)%pt = dyn(1)%pt - dyn(2)%pt dyn(1)%q = dyn(1)%q - dyn(2)%q endif call dyn_stat(6, dyn(1), rc) if ( rc .ne. 0 ) then call die(myname,'cannot process dyn_stat') end if ! If requested write *.hdf file with a header from dyn(1) ! ------------- if ( trim(dyn_dout) .ne. 'NONE' ) then call dyn_put ( trim(dyn_dout), nymd_d, nhms_d, 0, dyn(1), rc, freq=freq, vectype=dyntype ) endif ! ------------- !_RT end if ! Clean up mess ! ------------- call dyn_clean ( dyn(1) ) call dyn_clean ( dyn(2) ) end do end do ! loop over files st = "rm -f " // trim(binfiles(1)) // " " // trim(binfiles(2)) rc = system(st) ! rc == 0 for success if (rc .ne. 0 ) then print *, "Unable to remove binary files." end if ! All done ! -------- call exit(0) CONTAINS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS ! !----------------------------------------------------------------------- !BOP ! !IROUTINE: Init_ --- Initialize dyn2dyn ! ! !DESCRIPTION: parses command line. ! ! !INTERFACE: ! subroutine Init_ ( dyntype, mfiles, files, dominmax ) implicit NONE integer, intent(out) :: dyntype ! 4=geos4, 5=geos5 integer, intent(in) :: mfiles ! max. number of eta files ! dynamics file names (eta) character*255, intent(out) :: files(mfiles) logical, intent(out) :: dominmax ! ! !REVISION HISTORY: ! 2002.01.07 E. Yeh Initial code. ! !EOP !BOC character*4, parameter :: myname = 'init' integer iret, i, iarg, argc, iargc character(len=255) :: etafile, argv logical dout dout = .false. dyn_dout = 'NONE' dyntype = 4 ! default is GEOS-4 files dominmax = .false. print * print *, ' ---------------------------------------------------------' print *, ' dyndiff - dynamic vector difference between two hdf files' print *, ' ---------------------------------------------------------' print * ! Parse command line ! ------------------ argc = iargc() if ( argc .lt. 1 ) call usage() iarg = 0 nfiles = 0 do i = 1, 32767 iarg = iarg + 1 if ( iarg .gt. argc ) exit call GetArg ( iarg, argv ) select case (argv) case ("-g5") dyntype = 5 case ("-o") dout = .true. if ( iarg+1 .gt. argc ) call usage() iarg = iarg + 1 call GetArg ( iArg, dyn_dout ) case ("-h") if ( iarg+1 .gt. argc ) call usage() case default nfiles = nfiles + 1 if ( nfiles .gt. mfiles ) call usage() files(nfiles) = argv end select end do if ( nfiles .lt. 1 ) call usage() if ( nfiles .eq. 1 ) then nfiles = 2 files(2) = files(1) dominmax = .true. endif if ( nfiles .gt. 2 ) call usage() ! Echo the parameters ! ------------------- print * print *, '------------------------------------------------------------------' print *, ' Eta Dynamics state files: ' do i = 1, nfiles print *, i, ': ', trim(files(i)) end do end subroutine Init_ !................................................................. subroutine usage() print * print *,'Usage: ' print * print *,' dyndiff.x [-h] etafile_1 etafile_2 [-o diff_file]' print * print *, 'where' print * print *, '-h Help (optional)' print *, '-g5 Treats files as GEOS-5 files' print * print * print *, 'Where etafile_1 and etafile_2 are two (required)' print *, 'input dynamics vector files in hybrid (eta) coordinates' print * print *, '-o diff_file Optional to save binary output' print *, ' where diff_file - output dynamics vector' print *, ' file in hybrid (eta) coordinates' print *, ' with timestamp from the etafile_1' print *, ' and 32-bit precision' print *, ' output = file1 - file2' print *, ' (default: difference is shown at the screen in ascii)' print * print *, 'Note: User can save the ascii output by the following command:' print *, ' dyndiff.x etafile_1 etafile_2 > OutPutFileName' call exit(1) end subroutine usage !................................................................. subroutine die ( myname, msg ) character(len=*) :: myname, msg write(*,'(a)') trim(myname) // ': ' // trim(msg) call exit(1) end subroutine die !................................................................. end program dyndiff