C +-======-+ C Copyright (c) 2003-2007 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ module m_GFIO_GetFld ! !USES: implicit none PRIVATE PUBLIC GFIO_GetFld interface GFIO_GetFld; module procedure & GFIO_Get_Fall_, ! GEOS-3 & GFIO_Get_2d_, ! GEOS-3 & GFIO_Get_3d_ ! GEOS-3 end interface interface Fix_Undef_In; module procedure & Fix_Undef_In_ end interface ! Internal Interfaces ! ------------------- interface Hflip_ module procedure HFlip2_ ! Horizontal flip of 2d fields module procedure HFlip3_ ! Horizontal flip of 3d fields end interface character(len=*), parameter :: myname = 'GFIO_GetFld' logical, parameter :: verb = .true. integer, parameter :: stdout = 6 integer, parameter :: stderr = 6 real , parameter :: UNDEF = 1.e15 real , parameter :: rtol = 1.e-3 ! 25Sep2007 Todling: Adjusted tol CONTAINS !------------------------------------------------------------------------- ! NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS ! !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: GFIO_Get_Fall_ --- ! ! !INTERFACE: ! subroutine GFIO_Get_Fall_ ( gfiofile, nymd, nhms, im, jm, mlev, plevs, & nvars2d, vname2d, fld2d, & nvars3d, vname3d, fld3d, & stat ) ! !USES: Implicit NONE ! ! !INPUT PARAMETERS: ! character*(*) gfiofile ! Input GFIO file name ! Synoptic time: integer nymd ! Year-month-day, e.g., 19971012 integer nhms ! hour-minute-sec, e.g., 120000 integer im integer jm integer mlev integer nvars2d integer nvars3d character*(*) vname2d(nvars2d) character*(*) vname3d(nvars3d) real plevs(mlev) ! !OUTPUT PARAMETERS: real fld2d(im,jm,nvars2d) real fld3d(im,jm,mlev,nvars3d) integer, intent(out), optional :: stat ! ! !DESCRIPTION: Read grided analysis fields using GFIO pre- and post- ! analysis files. ! ! !TO DO: ! replace all exits by output return code ! ! !REMARKS: ! ! !REVISION HISTORY: ! ! 21Nov2000 Todling Initial code. ! 24Feb2006 Todling Adapted to handle flipped fields from GEOS-5. ! 11Jul2007 Stassi/RT Fix for lon precision issue between G5 and G4 files. ! 16Jul2008 Todling Lon check now using 64-bin array ! !EOP !------------------------------------------------------------------------- !_BOC character(len=*), parameter :: myname_ = myname//'*'//'GFIO_Get_Fall_' ! Grid Information ! direction integer ier, year, month, day, hour integer nobs_oi, i, j, k, kk, ia, ios integer lm_max_max ! Maximum number of time stamps in file integer ngatts ! number of attributes on file integer nvars, nvars_in integer im_max, jm_max, mlev_max, lm_max integer n2dim, n3dim integer iv, icount integer ncid integer dimSize, dimId integer ncvid character*255 dimName integer kloc(1), levcnt logical usrlevs logical fliplon real*8 mxlon, mnlon real*4, allocatable :: lon32(:) real*8, allocatable :: lon64(:) real, allocatable :: lev(:) real, allocatable :: aux(:,:,:,:) ! Define constants ! ---------------- nvars = nvars2d + nvars3d usrlevs = .true. if ( present(stat) ) stat = 0 ! Open GFIO file ! -------------- im_max = im jm_max = jm mlev_max = mlev call GFIO_Open ( trim(gfiofile), 0, ncid, ier ) if ( ier .eq. 0 ) then call GFIO_DimInquire ( ncid, im_max, jm_max, mlev_max, lm_max, & nvars_in, ngatts, ier) if ( nvars .gt. nvars_in ) then write(stdout,'(a,2(a,i5))') myname_, & ': Warning, no. of variable in file: ', nvars_in, & ' smaller than requested no. of vars: ', nvars end if ! Check consistency of Gfio file ! ------------------------------ if (im_max .ne. im) then write(stderr,'(2a,i5)') myname_, & ': trying to read unallowed longitudinal dim' call exit(7) endif if (jm_max .ne. jm) then write(stderr,'(2a,i5)') myname_, & ': trying to read unallowed latitudinal dim' call exit(7) endif if (mlev_max .ne. mlev) then write(stdout,'(2a,i5)') myname_, & '-WARNING: trying to read unallowed number of levels' usrlevs = .false. endif else write(stdout,'(2a)') myname_, . ': No GFIO file for current date ' call flush (stdout) if ( present(stat) ) then stat = 99 return else call exit(7) end if endif ! OK, let's read data from file ! ============================= icount = 0 n2dim = im * jm n3dim = n2dim * mlev ! Get longitude array to check on GEOS-5 flipped fields ! ----------------------------------------------------- allocate ( lon64(im), stat=ier ) if (ier/=0) then write(stderr,'(2a,i10,2a)') myname_, 'Alloc(lon) error ' call flush (stderr); call exit(7) endif call ncdinq (ncid, 1, dimName, dimSize, ier) if (ier/=0)then write(stderr,'(2a)') myname_, 'Error(ncdinq)' call flush (stderr); call exit(7) endif dimId = ncvid (ncid, dimName, ier) if (ier/=0)then write(stderr,'(2a)') myname_, 'Error(ncvid)' call flush (stderr); call exit(7) endif ! Caution: will try to read as 64-bit first ... ! --------------------------------------------- call ncvgt ( ncid, dimId, 1, im, lon64, ier ) if ( ier .ne. 0 ) then write(stderr,'(2a,i10,2a)') myname_, & ': Gfio_GetVar, ier = ', ier, & ' for variable ', 'lon' call flush(stderr); call exit(7) endif ! ... if absurd lons, read as 32-bit ! ---------------------------------- mnlon = minval(lon64) mxlon = maxval(lon64) if ( mnlon < -200.0 .or. mxlon > 400.0 ) then allocate ( lon32(im), stat=ier ) call ncvgt ( ncid, dimId, 1, im, lon32, ier ) lon64=lon32 deallocate ( lon32 ) endif ! Is orientation of fields shifted by 180 degs from internal def? ! --------------------------------------------------------------- fliplon = .false. if ( abs(lon64(1)+180.0)0) then fld3d(:,:,k:k,iv:iv) = aux(:,:,j:j,1:1) levcnt = levcnt + 1 end if end do if ( levcnt.ne.mlev ) then write(stdout,'(a,2(i10))') myname_, & 'WARNING: Requested levs not found, levcnt = ', & levcnt, ' levs requested ', mlev end if icount = icount + 1 end do deallocate ( aux, lev, stat=ier ) end if ! Close GFIO file now ! ------------------- call Gfio_Close ( ncid, ier ) if ( ier .ne. 0 ) then write(stdout,'(2a,i7)') myname_, & ': WARNING, problems closing GFIO file, ier = ', ier end if write(stdout,'(2a,i12,2a)') myname_, & ': read ', icount, ' variables from ', trim(gfiofile) ! All done ! -------- return end subroutine GFIO_Get_Fall_ !EOC !------------------------------------------------------------------------- ! NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS ! !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: GFIO_Get_2d_ --- ! ! !INTERFACE: ! subroutine GFIO_Get_2d_ ( gfiofile, nymd, nhms, im, jm, & nvars, vname, fld ) ! !USES: Implicit NONE ! ! !INPUT PARAMETERS: ! character*(*) gfiofile ! Input GFIO file name ! Synoptic time: integer nymd ! Year-month-day, e.g., 19971012 integer nhms ! hour-minute-sec, e.g., 120000 integer im integer jm integer nvars character*(*) vname(nvars) ! !OUTPUT PARAMETERS: real fld(im,jm,nvars) ! ! !DESCRIPTION: Read grided analysis fields using GFIO pre- and post- ! analysis files. ! ! !TO DO: ! ! !REMARKS: ! ! !REVISION HISTORY: ! ! 21Nov2000 Todling Initial code. ! 24Feb2006 Todling Adapted to handle flipped fields from GEOS-5. ! 11Jul2007 Stassi/RT Fix for lon precision issue between G5 and G4 files. ! 16Jul2008 Todling Lon check now using 64-bin array ! !EOP !------------------------------------------------------------------------- !_BOC character(len=*), parameter :: myname_ = myname//'*'//'GFIO_Get_2d_' ! Grid Information ! direction integer ier, year, month, day, hour integer nobs_oi, i, j, k, kk, ia, ios integer lm_max_max ! Maximum number of time stamps in file integer ngatts ! number of attributes on file integer nvars_in integer im_max, jm_max, mlev_max, lev_max, lm_max integer n2dim integer iv, icount integer ncid integer ncvid integer dimId, dimSize character*255 dimName real*8 mxlon, mnlon real*4, allocatable :: lon32(:) real*8, allocatable :: lon64(:) logical fliplon ! Open GFIO file ! -------------- im_max = im jm_max = jm call GFIO_Open ( trim(gfiofile), 0, ncid, ier ) if ( ier .eq. 0 ) then call GFIO_DimInquire ( ncid, im_max, jm_max, mlev_max, lm_max, & nvars_in, ngatts, ier) if ( nvars .gt. nvars_in ) then write(stdout,'(a,2(a,i5))') myname_, & ': Warning, no. of variable in file: ', nvars_in, & ' smaller than requested no. of vars: ', nvars end if ! Check consistency of Gfio file ! ------------------------------ if (im_max .ne. im) then write(stderr,'(2a,i5)') myname_, & ': trying to read unallowed longitudinal dim' call exit(7) endif if (jm_max .ne. jm) then write(stderr,'(2a,i5)') myname_, & ': trying to read unallowed latitudinal dim' call exit(7) endif else write(stdout,'(2a)') myname_, . ': No GFIO file for current date ' call flush (stdout) call exit(7) endif ! OK, let's read data from file ! ============================= icount = 0 n2dim = im * jm ! Get longitude array to check on GEOS-5 flipped fields ! ----------------------------------------------------- allocate ( lon64(im), stat=ier ) if (ier/=0) then write(stderr,'(2a,i10,2a)') myname_, 'Alloc(lon) error ' call flush (stdout); call exit(7) endif call ncdinq (ncid, 1, dimName, dimSize, ier) if (ier/=0)then write(stderr,'(2a)') myname_, 'Error(ncdinq)' call flush (stderr); call exit(7) endif dimId = ncvid (ncid, dimName, ier) if (ier/=0)then write(stderr,'(2a)') myname_, 'Error(ncvid)' call flush (stderr); call exit(7) endif ! Caution: will try to read as 64-bit first ... ! --------------------------------------------- call ncvgt ( ncid, dimId, 1, im, lon64, ier ) if ( ier .ne. 0 ) then write(stderr,'(2a,i10,2a)') myname_, & ': Gfio_GetVar, ier = ', ier, & ' for variable ', 'lon' call flush(stderr); call exit(7) endif ! ... if absurd lons, read as 32-bit ! ---------------------------------- mnlon = minval(lon64) mxlon = maxval(lon64) if ( mnlon < -200.0 .or. mxlon > 400.0 ) then allocate ( lon32(im), stat=ier ) call ncvgt ( ncid, dimId, 1, im, lon32, ier ) lon64=lon32 deallocate ( lon32 ) endif ! Is orientation of fields shifted by 180 degs from internal def? ! --------------------------------------------------------------- fliplon = .false. if ( abs(lon64(1)+180.0) 400.0 ) then allocate ( lon32(im), stat=ier ) call ncvgt ( ncid, dimId, 1, im, lon32, ier ) lon64=lon32 deallocate ( lon32 ) endif ! Is orientation of fields shifted by 180 degs from internal def ! ------------------------------------------------------------- fliplon = .false. if ( abs(lon64(1)+180.0)0) then fld(:,:,k:k,iv:iv) = aux(:,:,j:j,1:1) levcnt = levcnt + 1 end if end do if ( levcnt.ne.mlev ) then if(verbose_) write(stderr,'(a,2(i10))') myname_, & 'WARNING: Requested levs not found, levcnt = ', & levcnt, ' levs requested ', mlev end if else fld(:,:,1:mlev,iv:iv) = aux(:,:,1:mlev,1:1) endif ! < present(plevs) > icount = icount + 1 end do deallocate ( aux, lev, stat=ier ) end if ! Close GFIO file now ! ------------------- call Gfio_Close ( ncid, ier ) if ( ier .ne. 0 ) then if(verbose_) write(stdout,'(2a,i7)') myname_, & ': WARNING, problems closing GFIO file, ier = ', ier end if if(verbose_) write(stdout,'(2a,i12,2a)') myname_, & ': read ', icount, ' variables from ', trim(gfiofile) if(verbose_) write(stdout,'(2a,i8.8,a,i6.6)') myname_, & ': nymd, nhms = ', nymd,' ', nhms ! All done ! -------- return end subroutine GFIO_Get_3d_ !EOC !------------------------------------------------------------------------- ! NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS ! !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: Fix_Undef_In_ - Fixes undef values when input is 32 bit ! ! !INTERFACE: ! subroutine Fix_Undef_In_ ( ndim, var, vname, verbose ) ! ! !USES: use m_fpe, only: isnan implicit none ! !INPUT PARAMETERS: integer ndim character*(*) vname logical, intent(in), optional :: verbose ! !INPUT/OUTPUT PARAMETERS: real var(ndim) ! !DESCRIPTION: ! ! !REVISION HISTORY: ! ! 13Apr00 Todling - Initial code. ! 23Jun2010 Kokron - eliminate previous calculation of val_max_in ! because it would fail if the input field ! contained NaNs ! !EOP !------------------------------------------------------------------------- character(len=*), parameter :: myname_ = myname//'*'//'Fix_Undef_In_' integer i integer nonmiss, numnan real tol real val_max_in, val_max_out logical verbose_ tol = 0.001 * undef verbose_ = verb if ( present(verbose) ) verbose_ = verbose ! Make sure there is no precision problem with undef's ! ---------------------------------------------------- nonmiss = 0; numnan = 0.; val_max_in=0. do i= 1,ndim if (isnan(var(i))) then numnan=numnan+1 var(i) = undef elseif ( abs(var(i)-undef) .lt. tol ) then nonmiss = nonmiss + 1 var(i) = undef end if end do val_max_out = maxval(var) if (nonmiss.ne.0) then if(verbose_) & write(stdout,'(2a,i10,3a)') myname_, ': Fixed ', nonmiss, . ' values from input with strange undef', . ' for variable ', trim(vname) call flush(stdout) if ( val_max_out .ne. undef ) then write(stdout,*) ' Largest value on input: ', val_max_in write(stdout,*) ' Largest value on output: ', val_max_out write(stdout,*) ' Undef value spec.: ', undef write(stdout,*) ' Correction not done. Aborting ... ' call exit (7) end if end if if (numnan .ne. 0) then write(stdout,'(2a,i10,3a)') myname_, ': There were ',numnan, . ' NaN values from input', . ' for variable ', trim(vname) call flush(stdout) endif return end subroutine Fix_Undef_In_ subroutine hflip3_ ( q,im,jm,km ) implicit none integer im,jm,km,i,j,k real, intent(inout) :: q(im,jm,km) real, allocatable :: dum(:) allocate ( dum(im) ) do k=1,km do j=1,jm do i=1,im/2 dum(i) = q(i+im/2,j,k) dum(i+im/2) = q(i,j,k) enddo q(:,j,k) = dum(:) enddo enddo deallocate ( dum ) end subroutine hflip3_ subroutine hflip2_ ( q,im,jm ) implicit none integer im,jm,i,j real, intent(inout) :: q(im,jm) real, allocatable :: dum(:) allocate ( dum(im) ) do j=1,jm do i=1,im/2 dum(i) = q(i+im/2,j) dum(i+im/2) = q(i,j) enddo q(:,j) = dum(:) enddo deallocate ( dum ) end subroutine hflip2_ end module m_GFIO_GetFld