text:
	SUBROUTINE BFRTCH ( trange, time, ok, iret )
C************************************************************************
C* BFRTCH								*
C*									*
C* This subroutine checks for a station within time bounds.		*
C*									*
C* BFRTCH  ( TRANGE, TIME, OK, IRET )					*
C*									*
C* Input parameters:							*
C*	TRANGE (2)	CHAR*		Time range			*
C*	TIME		CHAR*		Station time			*
C*									*
C* Output parameters:							*
C*	OK		LOGICAL		Flag (true if station is in)	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal			*
C**									*
C* Log:									*
C* K. Brill/EMC		 2/97						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	trange(2), time
	LOGICAL		ok
C*
C-----------------------------------------------------------------------
	iret = 0
	ok = .false.
C*
	IF ( time .ge. trange (1) .and. time .le. trange (2) ) THEN
	    ok = .true.
	    RETURN
	END IF
	IF ( time .le. trange (1) .and. time .ge. trange (2) ) THEN
	    ok = .true.
	    RETURN
	END IF
	RETURN
	END


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
module psb_d_serial_mod
  use psb_const_mod
  use psb_error_mod

  interface psb_amax
    function psb_damax_s(n, x) result(val)
      import :: psb_ipk_, psb_dpk_
      integer(psb_ipk_), intent(in)  :: n
      real(psb_dpk_), intent(in)  :: x(:)
      real(psb_dpk_)                 :: val
    end function psb_damax_s
  end interface psb_amax

  interface psb_asum
    function psb_dasum_s(n, x) result(val)
      import :: psb_ipk_, psb_dpk_
      integer(psb_ipk_), intent(in)  :: n
      real(psb_dpk_), intent(in)  :: x(:)
      real(psb_dpk_)                 :: val
    end function psb_dasum_s
  end interface psb_asum

  interface psb_spspmm
    subroutine psb_dspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_dspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_dspmat_type), intent(in)  :: a,b
      type(psb_dspmat_type), intent(out) :: c
      integer(psb_ipk_), intent(out)                :: info
    end subroutine psb_dspspmm
    subroutine psb_dcsrspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_d_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_d_csr_sparse_mat), intent(in) :: a,b
      type(psb_d_csr_sparse_mat), intent(out) :: c
      integer(psb_ipk_), intent(out)          :: info
    end subroutine psb_dcsrspspmm
    subroutine psb_dcscspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_d_csc_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_d_csc_sparse_mat), intent(in) :: a,b
      type(psb_d_csc_sparse_mat), intent(out) :: c
      integer(psb_ipk_), intent(out)          :: info
    end subroutine psb_dcscspspmm
  end interface

  interface psb_symbmm
    subroutine psb_dsymbmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_dspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_dspmat_type), intent(in)  :: a,b
      type(psb_dspmat_type), intent(out) :: c
      integer(psb_ipk_), intent(out)                :: info
    end subroutine psb_dsymbmm
    subroutine psb_dbase_symbmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_d_base_sparse_mat, psb_d_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_d_base_sparse_mat), intent(in) :: a,b
      type(psb_d_csr_sparse_mat), intent(out)  :: c
      integer(psb_ipk_), intent(out)                     :: info
    end subroutine psb_dbase_symbmm
  end interface psb_symbmm

  interface psb_numbmm
    subroutine psb_dnumbmm(a,b,c)
      use psb_d_mat_mod, only : psb_dspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_dspmat_type), intent(in) :: a,b
      type(psb_dspmat_type), intent(inout)  :: c
    end subroutine psb_dnumbmm
    subroutine psb_dbase_numbmm(a,b,c)
      use psb_d_mat_mod, only : psb_d_base_sparse_mat, psb_d_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_d_base_sparse_mat), intent(in) :: a,b
      type(psb_d_csr_sparse_mat), intent(inout)  :: c
    end subroutine psb_dbase_numbmm
  end interface psb_numbmm

  interface psb_rwextd
    subroutine psb_drwextd(nr,a,info,b,rowscale)
      use psb_d_mat_mod, only : psb_dspmat_type
      import :: psb_ipk_
      implicit none
      integer(psb_ipk_), intent(in)                 :: nr
      type(psb_dspmat_type), intent(inout)        :: a
      integer(psb_ipk_),intent(out)                 :: info
      type(psb_dspmat_type), intent(in), optional :: b
      logical,intent(in), optional                 :: rowscale
    end subroutine psb_drwextd
    subroutine psb_dbase_rwextd(nr,a,info,b,rowscale)
      use psb_d_mat_mod, only : psb_d_base_sparse_mat
      import :: psb_ipk_
      implicit none
      integer(psb_ipk_), intent(in)                        :: nr
      class(psb_d_base_sparse_mat), intent(inout)        :: a
      integer(psb_ipk_),intent(out)                        :: info
      class(psb_d_base_sparse_mat), intent(in), optional :: b
      logical,intent(in), optional                       :: rowscale
    end subroutine psb_dbase_rwextd
    subroutine psb_ldrwextd(nr,a,info,b,rowscale)
      use psb_d_mat_mod, only : psb_ldspmat_type
      import :: psb_ipk_, psb_lpk_
      implicit none
      integer(psb_lpk_), intent(in)                  :: nr
      type(psb_ldspmat_type), intent(inout)        :: a
      integer(psb_ipk_),intent(out)                  :: info
      type(psb_ldspmat_type), intent(in), optional :: b
      logical,intent(in), optional                 :: rowscale
    end subroutine psb_ldrwextd
    subroutine psb_ldbase_rwextd(nr,a,info,b,rowscale)
      use psb_d_mat_mod, only : psb_ld_base_sparse_mat
      import :: psb_ipk_, psb_lpk_
      implicit none
      integer(psb_lpk_), intent(in)                         :: nr
      class(psb_ld_base_sparse_mat), intent(inout)        :: a
      integer(psb_ipk_),intent(out)                         :: info
      class(psb_ld_base_sparse_mat), intent(in), optional :: b
      logical,intent(in), optional                       :: rowscale
    end subroutine psb_ldbase_rwextd
  end interface psb_rwextd

  interface psb_geprt
    subroutine psb_dgeprtn2(fname,a,head)
      use psb_const_mod, only : psb_dpk_, psb_ipk_
      character(len=*), intent(in)  :: fname
      real(psb_dpk_), intent(in)    :: a(:,:)
      character(len=*), optional    :: head
    end subroutine psb_dgeprtn2
    subroutine psb_dgeprtn1(fname,a,head)
      use psb_const_mod, only : psb_dpk_, psb_ipk_
      character(len=*), intent(in)  :: fname
      real(psb_dpk_), intent(in)    :: a(:)
      character(len=*), optional    :: head
    end subroutine psb_dgeprtn1
    subroutine psb_dgeprt2(iout,a,head)
      use psb_const_mod, only : psb_dpk_, psb_ipk_
      integer(psb_ipk_), intent(in)            :: iout
      real(psb_dpk_), intent(in)     :: a(:,:)
      character(len=*), optional     :: head
    end subroutine psb_dgeprt2
    subroutine psb_dgeprt1(iout,a,head)
      use psb_const_mod, only : psb_dpk_, psb_ipk_
      integer(psb_ipk_), intent(in)            :: iout
      real(psb_dpk_), intent(in)     :: a(:)
      character(len=*), optional     :: head
    end subroutine psb_dgeprt1
  end interface psb_geprt

  interface psb_csprt
    module procedure psb_dcsprt, psb_dcsprtn
  end interface psb_csprt

  interface psb_spdot_srtd
    function psb_d_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot)
      use psb_const_mod, only : psb_ipk_, psb_dpk_
      integer(psb_ipk_), intent(in) :: nv1,nv2
      integer(psb_ipk_), intent(in) :: iv1(*), iv2(*)
      real(psb_dpk_), intent(in) :: v1(*),v2(*)
      real(psb_dpk_)      :: dot
    end function psb_d_spdot_srtd
  end interface psb_spdot_srtd

  interface psb_spge_dot
    function psb_d_spge_dot(nv1,iv1,v1,v2) result(dot)
      use psb_const_mod, only : psb_ipk_, psb_dpk_
      integer(psb_ipk_), intent(in) :: nv1
      integer(psb_ipk_), intent(in) :: iv1(*)
      real(psb_dpk_), intent(in) :: v1(*),v2(*)
      real(psb_dpk_)      :: dot
    end function psb_d_spge_dot
  end interface psb_spge_dot

  interface psb_aspxpby
    subroutine psb_d_aspxpby(alpha, nx, ix, x, beta, y, info)
      use psb_const_mod, only : psb_ipk_, psb_dpk_
      integer(psb_ipk_), intent(in)               :: nx
      integer(psb_ipk_), intent(in)               :: ix(:)
      real(psb_dpk_), intent (in)       :: x(:)
      real(psb_dpk_), intent (inout)    :: y(:)
      real(psb_dpk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_d_aspxpby
  end interface psb_aspxpby

  interface psb_spspmm
    subroutine psb_ldspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_ldspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_ldspmat_type), intent(in)  :: a,b
      type(psb_ldspmat_type), intent(out) :: c
      integer(psb_ipk_), intent(out)                :: info
    end subroutine psb_ldspspmm
    subroutine psb_ldcsrspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_ld_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_ld_csr_sparse_mat), intent(in) :: a,b
      type(psb_ld_csr_sparse_mat), intent(out) :: c
      integer(psb_ipk_), intent(out)          :: info
    end subroutine psb_ldcsrspspmm
    subroutine psb_ldcscspspmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_ld_csc_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_ld_csc_sparse_mat), intent(in) :: a,b
      type(psb_ld_csc_sparse_mat), intent(out) :: c
      integer(psb_ipk_), intent(out)          :: info
    end subroutine psb_ldcscspspmm
  end interface psb_spspmm

  interface psb_symbmm
    subroutine psb_ldsymbmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_ldspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_ldspmat_type), intent(in)  :: a,b
      type(psb_ldspmat_type), intent(out) :: c
      integer(psb_ipk_), intent(out)                :: info
    end subroutine psb_ldsymbmm
    subroutine psb_ldbase_symbmm(a,b,c,info)
      use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_ld_base_sparse_mat), intent(in) :: a,b
      type(psb_ld_csr_sparse_mat), intent(out)  :: c
      integer(psb_ipk_), intent(out)                     :: info
    end subroutine psb_ldbase_symbmm
  end interface psb_symbmm

  interface psb_numbmm
    subroutine psb_ldnumbmm(a,b,c)
      use psb_d_mat_mod, only : psb_ldspmat_type
      import :: psb_ipk_
      implicit none
      type(psb_ldspmat_type), intent(in) :: a,b
      type(psb_ldspmat_type), intent(inout)  :: c
    end subroutine psb_ldnumbmm
    subroutine psb_ldbase_numbmm(a,b,c)
      use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat
      import :: psb_ipk_
      implicit none
      class(psb_ld_base_sparse_mat), intent(in) :: a,b
      type(psb_ld_csr_sparse_mat), intent(inout)  :: c
    end subroutine psb_ldbase_numbmm
  end interface psb_numbmm

contains

  subroutine psb_dcsprt(iout,a,iv,head,ivr,ivc)
    use psb_d_mat_mod, only : psb_dspmat_type
    integer(psb_ipk_), intent(in)       :: iout
    type(psb_dspmat_type), intent(in) :: a
    integer(psb_lpk_), intent(in), optional :: iv(:)
    character(len=*), optional    :: head
    integer(psb_lpk_), intent(in), optional :: ivr(:),ivc(:)

    call a%print(iout,iv,head,ivr,ivc)

  end subroutine psb_dcsprt

  subroutine psb_dcsprtn(fname,a,iv,head,ivr,ivc)
    use psb_d_mat_mod, only : psb_dspmat_type
    character(len=*), intent(in)  :: fname
    type(psb_dspmat_type), intent(in) :: a
    integer(psb_lpk_), intent(in), optional :: iv(:)
    character(len=*), optional    :: head
    integer(psb_lpk_), intent(in), optional :: ivr(:),ivc(:)

    call a%print(fname,iv,head,ivr,ivc)

  end subroutine psb_dcsprtn

end module psb_d_serial_mod


text:
	SUBROUTINE SNDCMP  ( nparms, rlevel, nlev, levtyp, lvert,
     +			     sndata, numlev, parms, wavlen, wavspd,
     +			     storm, crdrot, tropht, trpint, cldhgt,
     +			     mxdpth, squall, delz, filtyp, spline,
     +			     nlun, lun, dattim, nlevel, hdata, iret )
C************************************************************************
C* SNDCMP								*
C*									*
C* This subroutine gets the requested sounding data and computes the	*
C* parameters at all the station levels.				*
C*									*
C* SNDCMP  ( NPARMS, RLEVEL, NLEV, LEVTYP, LVERT, SNDATA, NUMLEV, 	*
C*	     PARMS, WAVLEN, WAVSPD, STORM, CRDROT, TROPHT, TRPINT, 	*
C*	     CLDHGT, MXDPTH, SQUALL, DELZ, FILTYP, SPLINE, NLUN, LUN,	*
C*	     DATTIM, NLEVEL, HDATA, IRET )				*
C*									*
C* Input parameters:							*
C*	NPARMS		INTEGER		Number of parameters		*
C*	RLEVEL (LLMXLV)	REAL		Data levels			*
C*	NLEV		INTEGER		Number of levels		*
C*	LEVTYP		INTEGER		Level type			*
C*	LVERT		INTEGER		Vertical coordinate		*
C*	SNDATA (LLMXLV)	REAL		Original sound data		*
C*	NUMLEV		INTEGER		# of original sounding levels	*
C*	PARMS (NPARMS)	CHAR*		Parameters			*
C*	WAVLEN		CHAR*		Length of wave			*
C*	WAVSPD		CHAR*		Speed of wave			*
C*	STORM		CHAR*		Storm speed and direction	*
C*	CRDROT		CHAR*		Coordinate rotations		*
C*	TROPHT		CHAR*		Tropopause height		*
C*	TRPINT		CHAR*		Interval above and below trop	*
C*	CLDHGT		CHAR*		Height of cloud base		*
C*	MXDPTH		CHAR*		Height of mixed layer depth	*
C*	SQUALL		CHAR*		Length of squall line		*
C*	DELZ		CHAR*		Height interval for new sounding*
C*	FILTYP		CHAR*		Filter type			*
C*	SPLINE		LOGICAL		Whether to use splines or not	*
C*	NLUN		INTEGER		Number of file numbers		*
C*	LUN (NLUN)	INTEGER		Output file number		*
C*	DATTIM		CHAR*		Time				*
C*									*
C* Output parameters:							*
C*	NLEVEL		INTEGER		Number of new sounding levels	*
C*	HDATA (LLMXLV)	REAL		New sounding data		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal			*
C**									*
C* Log:									*
C* S. Jacobs/SSAI	 4/92		Copied from SNOCMP		*
C* J. Whistler/SSAI	 4/93		Cleaned up			*
C* J. Whistler/SSAI	 5/93		Filter all of the data		*
C* J. Whistler/SSAI	 6/93		Added check to see if DWPT is	*
C*					interpolated higher than TEMP	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'sndiag.prm'
C*
	PARAMETER	( NUMPTL = 8, NUMPTH = 15 )
	PARAMETER	( NMARL = NUMPTL+1, NMARH = NUMPTH+1 )
C*
	CHARACTER*(*)	wavlen, wavspd, storm, crdrot, tropht, trpint,
     +			cldhgt, mxdpth, squall, delz,  filtyp, parms(*),
     +			dattim
	REAL		rlevel(*), sndata(*), hdata(*)
	LOGICAL		spline
	INTEGER		lun(*)
C*
	CHARACTER	cdata(LLMXDT)
	REAL		hghts(LLMXDT), rdata(LLMXDT),
     +			thtv(LLMXLV), uwnd(LLMXLV), vwnd(LLMXLV),
     +			uabs(LLMXLV), vabs(LLMXLV),
     +			filarl(NMARL), filarh(NMARH)
	LOGICAL		chrflg(MMPARM), cmpflg(MMPARM)
C*
	DATA		filarl /  .3200156553,  .2643210644,
     +				  .1343767864,  .0114831529,
     +				 -.0458900942, -.0377520979,
     +				 -.0064753717,  .0111929411,
     +				  .0087357913 /
C
	DATA		filarh /  .0974726419,  .0954676702,
     +				  .0896329731,  .0804876892,
     +				  .0688283154,  .0556356004,
     +				  .0419626250,  .0288191067,
     +				  .0170670479,  .0073407153,
     +				 0.0000000000, -.0048789007,
     +				 -.0074769889, -.0081766107,
     +				 -.0074864031, -.0059591606 /
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret   = 0
	levout = 0
	knt    = 1
C
C*	Check for list of levels or range.
C
	IF  ( levtyp .eq. 1 )  THEN
C
C*	    Always compute surface data first.
C
	    CALL PC_CMVR  ( 0., lvert, sndata, rdata (knt),
     +			    cdata (knt), ier )
	    IF  ( rlevel (1) .eq. 0. )  THEN
		istart = 2
	      ELSE
		istart = 1
	    END IF
C
C*	    Loop through the list of levels.
C
	    levout = 1
	    knt    = nparms + 1
	    DO  lev = istart, nlev
		CALL PC_CMVR  ( rlevel (lev), lvert, sndata,
     +				rdata (knt), cdata (knt), ier )
		IF  ( .not. ERMISS (rdata (knt) ) )  THEN
		    levout = levout + 1
		    knt    = knt + nparms
		END IF
	    END DO
C
C*	    Check for surface as only data.
C
	    IF  ( ( levout .eq. 1 ) .and. ( rlevel (1) .ne. 0. ) )  THEN
		levout = 0
	    END IF
C
C*	    Check for range with increment.
C
	ELSE IF  ( levtyp .eq. 2 )  THEN
C
C*	    Find bottom level.
C
	    CALL PC_FLVL  ( rlevel (1), lvert, sndata, begin, lvnum1,
     +			    lvnum2, levloc, ier )
	    IF  ( levloc .eq. 1 )  THEN
		ilev = lvnum1
	      ELSE IF  ( levloc .eq. 2 )  THEN
		ilev = lvnum2
	      ELSE IF  ( levloc .eq. 3 )  THEN
		ilev = 1
	      ELSE
		ilev = numlev + 1
	    END IF
C
C*	    Find top level.
C
	    CALL PC_FLVL  ( rlevel (2), lvert, sndata, end, lvnum1,
     +			    lvnum2, levloc, ier )
	    IF  ( levloc .eq. 1 )  THEN
		jlev = lvnum1
	      ELSE IF  ( levloc .eq. 2 )  THEN
		jlev = lvnum1
	      ELSE IF  ( levloc .eq. 4 )  THEN
		jlev = numlev
	      ELSE
		jlev = 0
	    END IF
C
C*	    Always compute surface data first.
C
	    CALL PC_CMVR  ( 0., lvert, sndata, rdata (knt),
     +			    cdata (knt), ier )
	    IF  ( ilev .eq. 1 )  THEN
		iilev = 2
	      ELSE
		iilev = ilev
	    END IF
C
C*	    Compute data at each level.
C
	    levout = 1
	    knt    = nparms + 1
	    DO  lev = iilev, jlev
		CALL PC_CMLV  ( lev, sndata, rdata (knt), cdata (knt),
     +				ier )
		IF  ( .not. ERMISS  ( rdata (knt) ) )  THEN
		    levout = levout + 1
		    knt    = knt + nparms
		END IF
	    END DO
C
C*	    Check for surface as only data.
C
	    IF  ( ( levout .eq. 1 ) .and. ( ilev .ne. 1 ) )  THEN
		levout = 0
	    END IF
	END IF
C
C*	Change missing dewpoints to 40 degrees less than the TEMP.
C
	ll = 1
	DO  n = 1, levout
	    IF  ( ERMISS(rdata((n-1)*nparms+ITEMP)) .and.
     +		  ERMISS(rdata((n-1)*nparms+IDWPT)) )  THEN
	    ELSE
		rdata((ll-1)*nparms+IHGHT) = rdata((n-1)*nparms+IHGHT)
		rdata((ll-1)*nparms+IPRES) = rdata((n-1)*nparms+IPRES)
		rdata((ll-1)*nparms+ITEMP) = rdata((n-1)*nparms+ITEMP)
		rdata((ll-1)*nparms+IDWPT) = rdata((n-1)*nparms+IDWPT)
		rdata((ll-1)*nparms+IUWND) = rdata((n-1)*nparms+IUWND)
		rdata((ll-1)*nparms+IVWND) = rdata((n-1)*nparms+IVWND)
		rdata((ll-1)*nparms+ITHTV) = rdata((n-1)*nparms+ITHTV)
c	write(6,'(7f9.2)' ) rdata((ll-1)*nparms+IHGHT),
c     +			rdata((ll-1)*nparms+IPRES),
c     +			rdata((ll-1)*nparms+ITEMP),
c     +			rdata((ll-1)*nparms+IDWPT),
c     +			rdata((ll-1)*nparms+IUWND),
c     +			rdata((ll-1)*nparms+IVWND),
c     +			rdata((ll-1)*nparms+ITHTV)
	    	IF  ( ERMISS(rdata((ll-1)*nparms+IDWPT)) )
     +		    rdata((ll-1)*nparms+IDWPT) =
     +			rdata((ll-1)*nparms+ITEMP) - 40.
		ll = ll + 1
	    END IF
	END DO
	levout = ll - 1
C
C*	Find average DELZ and put data at these levels.
C*	Also give the user a chance to exit.
C
	CALL SNDARR ( delz, nparms, levout, rdata,
     +		      nlevel, hghts, iret )
	IF  ( iret .ne. 0 )  RETURN
C
C*	Initialize the output array to missing.
C
	DO  i = 1, nlevel*nparms
	    hdata(i) = RMISSD
	END DO
C
	IF  ( spline )  THEN
C
C*	    Fit cubic splines to the data...
C
	    DO  n = 1, ITHTV
		CALL SNDSPL ( nparms, n, levout, rdata,
     +			      nlevel, hghts, hdata, iret )
	    END DO
	    IF  ( nparms .gt. NVAR-2 )  THEN
		CALL SNDSPL ( nparms, ILATI, levout, rdata,
     +			      nlevel, hghts, hdata, iret )
		CALL SNDSPL ( nparms, ILONI, levout, rdata,
     +			      nlevel, hghts, hdata, iret )
	    END IF
	ELSE
C
C*	    ...or interpolate to the new height surfaces.
C
	    CALL SNDINT ( nparms, levout, rdata,
     +			  nlevel, hghts, hdata, iret )
	END IF
C
C*	Check to see if Dewpoint is less than temperature.
C
	DO i = 1, nlevel
	    IF ((( .not. ERMISS ( hdata((i-1)*nparms+IDWPT) )) .and.
     +		 ( .not. ERMISS ( hdata((i-1)*nparms+ITEMP) ))) .and.
     +		( hdata((i-1)*nparms+IDWPT) .gt.
     +		  hdata((i-1)*nparms+ITEMP) ) )
     +		  hdata((i-1)*nparms+IDWPT) = hdata((i-1)*nparms+ITEMP)
	END DO
C
C*	Filter all of the data, if requested. (See below.)
C
	CALL ST_LCUC ( filtyp, filtyp, ier )
	IF  ( filtyp .eq. 'LIGHT' )  THEN
C
C*	    Apply a LIGHT or MESO-PASS filter.
C
C*	    Filter UWND.
C
	    CALL SNDFLT ( NUMPTL, filarl, IUWND, nparms, nlevel,
     +			  hdata, uabs, iret )
C
C*	    Filter VWND.
C
	    CALL SNDFLT ( NUMPTL, filarl, IVWND, nparms, nlevel,
     +			  hdata, vabs, iret )
C
	ELSE IF  ( filtyp .eq. 'HEAVY' )  THEN
C
C*	    Apply a HEAVY or MACRO-PASS filter.
C
C*	    Filter UWND.
C
	    CALL SNDFLT ( NUMPTH, filarh, IUWND, nparms, nlevel,
     +			  hdata, uabs, iret )
C
C*	    Filter VWND.
C
	    CALL SNDFLT ( NUMPTH, filarh, IVWND, nparms, nlevel,
     +			  hdata, vabs, iret )
C
	ELSE
C
C*	    Use NO filter.
C
	    DO  i = 1, nlevel
		uabs(i) = hdata((i-1)*nparms+IUWND)
		vabs(i) = hdata((i-1)*nparms+IVWND)
	    END DO
C
	END IF
	DO  i = 1, nlevel
	    hdata((i-1)*nparms+IUABS) = uabs(i)
	    hdata((i-1)*nparms+IVABS) = vabs(i)
	END DO
C
C*	Calculate the components of the shear for the interpolated
C*	winds, before computing the relative winds.
C
	CALL SNDSHR ( nparms, nlevel, delz, hdata, ier )
C
C*	Calculate relative winds.
C
	CALL SNDWND ( storm, crdrot, nparms, nlevel, hdata, iret )
C
C*	Filter all of the data.
C*	The reference for this procedure is:
C*	    Madden, et al., NCAR Tech note [NCAR-TN/STR-55] (1971).
C
	CALL ST_LCUC ( filtyp, filtyp, ier )
	DO iii = 1, 7
	    IF  ( filtyp .eq. 'LIGHT' )  THEN
C
C*	        Apply a LIGHT or MESO-PASS filter.
C
C*	        Filter THTV.
C
	        CALL SNDFLT ( NUMPTL, filarl, iii, nparms, nlevel,
     +			      hdata, thtv, iret )
	    ELSE IF  ( filtyp .eq. 'HEAVY' )  THEN
C
C*	        Apply a HEAVY or MACRO-PASS filter.
C
C*	        Filter THTV.
C
	        CALL SNDFLT ( NUMPTH, filarh, iii, nparms, nlevel,
     +			      hdata, thtv, iret )
C
	    ELSE
C
C*	        Use NO filter.
C
	        DO  i = 1, nlevel
		    thtv(i) = hdata((i-1)*nparms+iii)
	        END DO
C
	    END IF
C
C*	    Set the data in the large array to the filtered values.
C
	    DO  i = 1, nlevel
	        hdata((i-1)*nparms+iii) = thtv(i)
	    END DO
	END DO
C
C*	Check to see if Dewpoint is less than temperature.
C
	DO i = 1, nlevel
	    IF ((( .not. ERMISS ( hdata((i-1)*nparms+IDWPT) )) .and.
     +		 ( .not. ERMISS ( hdata((i-1)*nparms+ITEMP) ))) .and.
     +		( hdata((i-1)*nparms+IDWPT) .gt.
     +		  hdata((i-1)*nparms+ITEMP) ) )
     +		  hdata((i-1)*nparms+IDWPT) = hdata((i-1)*nparms+ITEMP)
	END DO
C
	jnumlv = nlevel
C
C*	Calculate the:	Intrinsic Frequency		INFQ
C*			Brunt-Vaisala/Intrinsic		BVIN
C*			Scorer Parameter		SCOR
C*			Term 1 of Scorer Parameter	SCR1
C*			Term 2 of Scorer Parameter	SCR2
C
	CALL SNDSCR ( nparms, nlevel, wavlen, wavspd, delz,
     +		      hdata, ier )
C
C*	Compute the mixed layer quantities.
C
	CALL SNDMIX ( nlun, lun, mxdpth, nparms, nlevel, hdata,
     +		      tlcl, plcl, thtwlcl, hlcl, iret )
C
C*	Compute the air and moisture flux.
C
	CALL SNDFLX ( squall, nparms, nlevel, hdata, iret )
C
C*	Compute the wind and lapse rate above and below the trop.
C
	CALL SNDTRP ( tropht, trpint, delz, nlun, lun, nparms, nlevel,
     +		      hdata, iret )
C
C*	Compute the layer properties of the wind.
C
	CALL SNDLYR ( tropht, delz, cldhgt, hlcl, nlun, lun,
     +		      nparms, nlevel, hdata, uabs, vabs, iret )
C
C*	Compute an estimation for the tropopause temp and hght.
C
	CALL SNDEST ( delz, nlun, lun, nparms, nlevel, hdata, iret )
C
C*	Compute the helicity at each level.
C
	CALL SNDHLC ( nparms, nlevel, hdata, iret )
C
C*	Compute the lid strength index, etc.
C
C	CALL SNDLID ( nlun, lun, dattim, nparms, levout, rdata, iret )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!OPTIONS: -flogical-abbreviations -fxor-operator

! Like m4 in resolve63 but compiled with different options.
! Alternate operators are enabled so treat these as intrinsic.
module m4
contains
  subroutine s1(x, y, z)
    logical :: x
    real :: y, z
    !ERROR: Operands of .AND. must be LOGICAL; have REAL(4) and REAL(4)
    x = y .a. z
    !ERROR: Operands of .OR. must be LOGICAL; have REAL(4) and REAL(4)
    x = y .o. z
    !ERROR: Operand of .NOT. must be LOGICAL; have REAL(4)
    x = .n. y
    !ERROR: Operands of .NEQV. must be LOGICAL; have REAL(4) and REAL(4)
    x = y .xor. z
    !ERROR: Operands of .NEQV. must be LOGICAL; have REAL(4) and REAL(4)
    x = y .x. y
  end
end

! Like m4 in resolve63 but compiled with different options.
! Alternate operators are enabled so treat .A. as .AND.
module m5
  interface operator(.A.)
    logical function f1(x, y)
      integer, intent(in) :: x, y
    end
  end interface
  interface operator(.and.)
    logical function f2(x, y)
      real, intent(in) :: x, y
    end
  end interface
contains
  subroutine s1(x, y, z)
    logical :: x
    complex :: y, z
    !ERROR: No user-defined or intrinsic .A. operator matches operand types COMPLEX(4) and COMPLEX(4)
    x = y .and. z
    !ERROR: No user-defined or intrinsic .A. operator matches operand types COMPLEX(4) and COMPLEX(4)
    x = y .a. z
  end
end


text:
program struct_tx

use, intrinsic :: iso_c_binding, only : c_int, c_bool, c_char, c_null_char

implicit none (type, external)

type, bind(C) :: my_struct
!! order and length must match in Fortran and C
integer(c_int) :: my_int
logical(c_bool) :: my_bool
character(kind=c_char) :: my_char(1000)
!! character(kind=c_char) in bind(c) type cannot be allocatable. Just have to make it "long enough"
!! or use iso_c_binding.h stuff

end type my_struct

interface
subroutine struct_check(p) bind(C)
import my_struct
type(my_struct), intent(in) :: p
end subroutine struct_check
end interface

type(my_struct) :: s
integer :: i
character(:), allocatable :: my_char

my_char = "Hello"

s%my_int = 123
s%my_bool = .true.

do i = 1, len(my_char)
  s%my_char(i) = my_char(i:i)
end do
s%my_char(i+1) = c_null_char

call struct_check(s)

print *, "OK: Fortran => C struct"
end program

max_stars_repo_name: scivision/fortran-c-cpp-interface

max_issues_repo_name: scivision/fortran-c-cpp-interface

max_forks_repo_name: scivision/fortran-c-cpp-interface


!--------+--------+--------+--------+--------+--------+--------+------!
! main program of a set of tools based on Wannier90 TB
! constructed by Q.S.Wu on 4/9/2010
! change      by Q.S.Wu on 4/22/2010
! changed     by Q.S.wu on July/15/2010
! version     HmnR.data  contains soc
! mpi-version is not test yet, take you own risk.
! mpi-version is  tested , please report bugs to QSWU
! Jan 25 2015 by Q.S.Wu at ETH Zurich
! wuquansheng@gmail.com
!--------+--------+--------+--------+--------+--------+--------+------!

  program main

    use para
    use mpi
    implicit none

    ! file existence
    logical  :: exists

     character*4 :: c_temp
     integer :: i_temp
	  real(dp) :: r_temp

     integer :: namelen
     integer :: ierr

     ! initial the environment of mpi
     call mpi_init(ierr)
     call mpi_comm_rank(mpi_cmw,cpuid,ierr)
     call mpi_comm_size(mpi_cmw,num_cpu,ierr)

     ! if mpi initial wrong, alarm
     if (cpuid==0.and.ierr.ne.0)then
        write(stdout,*)'mpi initialize wrong'
        stop
     endif

     call readinput

     ! open file Hmn_R.data to get Num_wann and Nrpts
     if (cpuid.eq.0)then
        write(stdout,*)''
        inquire (file =infilename, EXIST = exists)
        if (exists)then
           if (.not.index(infilename, 'HWR')) then
              write(stdout,'(2x,a,a,a)')'File ',trim(infilename), &
                 ' exist, We are using HmnR from wannier90'
              open(unit=1001,file=infilename,status='old')
              read(1001,*)
              read(1001,'(i)') Num_wann
              read(1001,'(i)') Nrpts
              write(stdout,*)'>> Num_wann', Num_wann
              write(stdout,*)'>> NRPTS', NRPTS
              close(1001)
           else
              write(stdout,'(2x, 3a)')'File ',trim(infilename), &
                 ' exist, We are using HmnR from HWR'
              open(unit=1001,file=infilename,status='old')
              read(1001,*)
              read(1001,'(a26, i3)') c_temp, Num_wann
              read(1001,'(a32,i10)')c_temp,Nrpts
              write(stdout,*)'>> Num_wann', Num_wann
              write(stdout,*)'>> NRPTS', NRPTS
              close(1001)
           endif ! hwr or not
        else
           write(stdout,'(2x,a25)')'>>> Error : no HmnR input'
           stop
        endif ! exists or not
        if ((soc==0 .and. sum(nprojs)/=Num_wann) .or. &
           (soc>0 .and. sum(nprojs)/=Num_wann/2))then
           print *, sum(nprojs), num_wann, num_wann/2
           stop 'projectors are wrong'
        endif
     endif ! cpuid

	  ! broadcast and Nrpts to every cpu
     call MPI_bcast(Num_wann,1,mpi_in,0,mpi_cmw,ierr)
     call MPI_bcast(Nrpts,1,mpi_in,0,mpi_cmw,ierr)

     !> dimension for surface green's function
     Ndim= Num_wann* Np

	  allocate(irvec(3,nrpts))
	  allocate(ndegen(nrpts))
     allocate(HmnR(num_wann,num_wann,nrpts))

     if(cpuid==0)then
       write(stdout,*) 'begin reading Hmn_R.data'
     endif
     call readHmnR()
     if(cpuid==0)then
       write(stdout,*) 'read Hmn_R.data successfully'
     endif

	  ! broadcast data to every cpu
     call MPI_bcast(irvec,size(irvec),mpi_in,0,mpi_cmw,ierr)
     call MPI_bcast(HmnR,size(HmnR),mpi_dc,0,mpi_cmw,ierr)
     call MPI_bcast(ndegen,size(ndegen),mpi_in,0,mpi_cmw,ierr)

     !> import symmetry
     call symmetry
     !> bulk band
	  if(cpuid.eq.0)write(stdout, *)'begin to calculate bulk band'
     if (BulkBand_calc) then
        call ek_bulk
       !call ek_bulk_mirror_x
       !call ek_bulk_mirror_z
       !call psik_bulk
       !call ek_bulk_polar
       !call ek_bulk_fortomas
       !call ek_bulk_spin
       !call dos_calc
       !call ek_bulk2D
       !call ek_bulk2D_spin
       !call fermisurface3D
       !call gapshape
       !call gapshape3D
       !call landau_level_k
       !call landau_level_B
     endif
     if(cpuid.eq.0)write(stdout, *)'end calculate bulk band'

     !> slab band
     if (SlabBand_calc)then
        call ek_slab
        call psik
     endif

     !> wannier center calculate
    !if (wanniercenter_calc)call wannier_center2D
    !if (wanniercenter_calc)call wannier_center2D_alt
    !if (wanniercenter_calc)call wannier_center3D
     if (wanniercenter_calc)then
        call wannier_center3D_plane
       !call wannier_center3D_plane_mirror_plus
       !call wannier_center3D_plane_mirror_minus
     endif
    !if (berry_calc)call berry_curvarture
     if (berry_calc)call berryphase

     !> surface state
	  if(cpuid.eq.0)write(stdout, *)'begin to calculate surface state'
     if (SlabSS_calc)call surfstat
     if(cpuid.eq.0)write(stdout, *)'end calculate surface state'

	  if(cpuid.eq.0)write(stdout, *)'begin to calculate surface state'
     if (WireBand_calc) then
   	  if(cpuid.eq.0)write(stdout, *)'begin to calculate ribbon band'
        call ek_ribbon
        if(cpuid.eq.0)write(stdout, *)'end calculate ribbon band'
     endif

     !> fermi arc
	  if(cpuid.eq.0)write(stdout, *)'begin to calculate fermi arc'
     if (SlabArc_calc)call fermiarc
     if(cpuid.eq.0)write(stdout, *)'end calculate fermi arc'

     !> calculate spin-texture
     if(cpuid.eq.0)write(stdout, *)'begin to calculate spin texture'
     if (SlabSpintexture_calc)call spintext
     if(cpuid.eq.0)write(stdout, *)'end calculate spin texture'

     if (cpuid.eq.0)write(stdout,*)'Congratulations! you finished the calculation.'

     call mpi_finalize(ierr)

  end


text:
        SUBROUTINE FLUSH(IDUM)
        INTEGER IDUM
        RETURN
        END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


module PackageWriterModule

  use ConstantsModule, only: LENFTYPE, LENPACKAGENAME, LINELENGTH, &
                             MAXCHARLEN
  use ConstantsPHMFModule, only: FCUNKNOWN, LENCTYPE
  use FileListModule, only: FileListType
  use FileTypeModule, only: FileType, ConstructFileType
  use FileWriterModule, only: FileWriterType
  use GLOBAL, only: NPER, NCOL, NROW, NLAY, IUNIT, IOUT, DELR, DELC, IBOUND
  use GlobalVariablesModule, only: echo, masteridomain, verbose
  use LineListModule, only: LineListType, same_lines
  use ListModule, only: ListType
  use ModelPackageModule, only: ModelPackageType, ConstructModelPackageType
  use MoverModule, only: MoverType, AddMoverToList
  use ObsWriterModule, only: ObsWriterType
  use SimModule, only: store_error, store_note, ustop
  use SimVariablesModule, only: SimMovers
  use UtilitiesModule, only: ConstantReal2D

  implicit none

  private
  public :: PackageWriterType, CastAsPackageWriterClass

  type, extends(FileWriterType) :: PackageWriterType
    integer, pointer :: ICbc => null()       ! CBC unit number
    integer, pointer :: ipr => null()        ! print flag
    integer, pointer :: NBndPeriod => null() ! number of boundaries used in stress period
    integer, pointer :: IPB => null()
    integer, pointer :: NBNDVL => null()
    integer, pointer :: Nlaynew => null()
    integer, dimension(:), pointer :: Layptr => null()
    integer, dimension(:,:,:), pointer :: Idomain => null()
    integer :: Igrid = 0
    integer :: IuOrig = 0        ! unit number for package input file of MF2005 model
    integer :: MaxActiveBnd = -1
    integer :: NAux = 0          ! number of auxiliary variables
    integer :: nstart = 4        ! default: first element in rlist following k,i,j
    integer :: nstop = 0         ! last element in rlist to be read
    integer :: nvalp1            ! Rlist index for element following minimum required data
    integer :: FieldsPerLine     ! number of fields available per line of array output
    integer :: DataWidth = 14    ! characters available for each array data value
    integer :: NStressDim  ! Number of values after k,i,j required to specify boundary stress
    integer :: Iouts = 0  ! For solute output? Never assigned in mf2005
    integer :: Istcb2 = 0
    integer :: IunitBcf = 0
    integer :: IunitGwt = 0
    integer :: IunitHuf = 0
    integer :: IunitLak = 0
    integer :: IunitLpf = 0
    integer :: IunitSfr = 0
    integer :: IunitUpw = 0
    integer :: IunitUzf = 0
    integer :: iulist = 0
    integer :: Nsol = 0
    character(len=MAXCHARLEN)     :: ModelBasename
    character(len=LENPACKAGENAME) :: PackageName = ''
    character(len=LENPACKAGENAME) :: DefaultBudgetText = ''
    character(len=6) :: source = ''
    character(len=100) :: fmat      ! format for printing boundary data
    character(len=MAXCHARLEN) :: TsFile = ''
    character(len=4) :: OrigPkg = ''
    character(len=16), pointer, dimension(:)  :: Aux => null()
    double precision, pointer, dimension(:,:) :: rlist => null()
    double precision, pointer, dimension(:,:) :: CellArea => null()
    double precision, pointer, dimension(:,:) :: work => null()
    logical :: PrintInput = .true.
    logical :: PrintFlows = .true.
    logical :: SaveFlows = .false.
    logical :: NeedDimensionsBlock = .true.
    logical :: NeedWaterMover = .false.
    logical :: Newton = .false.
    type(FileListType), pointer     :: Mf6Files => null()
    type(ListType), pointer         :: ModelMovers => null()
    type(ModelPackageType), pointer :: ModelPack => null()
    type(LineListType), pointer     :: BlockA => null()
    type(LineListType), pointer     :: BlockB => null()
    type(LineListType), pointer     :: CurrentBlock => null()
    type(LineListType), pointer     :: PreviousBlock => null()
    ! for observations
    logical :: ObsActive = .false.
    class(ObsWriterType), pointer   :: PkgObsWriter => null()
  contains
    procedure :: ProcessAllocate
    procedure :: ProcessStressLoop
    procedure :: MyType
    procedure :: AddMover
    procedure :: AllocatePointers
    procedure :: BuildModelMovers
    procedure :: BuildSimMovers
    procedure :: GetBudgetText
    procedure :: SameType
    procedure :: WriteOptions
    procedure :: WriteDimensions
    procedure :: WriteStressPeriodListData
    procedure :: WriteStressPeriodArrayData
    procedure :: WriteBlockIfNeeded
    procedure :: WriteArrayDataIfNeeded
  end type PackageWriterType

contains

  !=============================================================================
  ! Procedures that need to be overridden by types
  ! that inherit from PackageWriterType.

  subroutine ProcessAllocate(this, igrid)
    class(PackageWriterType) :: this
    integer, intent(in) :: igrid
    stop 'Error: PackageWriterType%ProcessAllocate not overridden in child type.'
  end subroutine ProcessAllocate

  subroutine ProcessStressLoop(this, igrid)
    class(PackageWriterType) :: this
    integer, intent(in) :: igrid
    stop 'Error: PackageWriterType%ProcessStressLoop not overridden in child type.'
  end subroutine ProcessStressLoop

  function MyType(this) result(ctype)
    class(PackageWriterType) :: this
    character(len=30) :: ctype
    ctype = ' '
    stop 'Error: PackageWriterType%MyType not overridden in child type.'
  end function MyType

  !=============================================================================

  subroutine AllocatePointers(this)
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    ! local
    type(FileType), pointer :: newFile
    !
    allocate(this%ICbc)
    this%ICbc = 0
    allocate(this%ipr)
    allocate(this%NBndPeriod)
    allocate(this%IPB)
    allocate(this%NBNDVL)
    allocate(this%Nlaynew)
    !
    allocate(this%BlockA)
    allocate(this%BlockB)
    allocate(this%CurrentBlock)
    allocate(this%PreviousBlock)
    !
    call this%BlockA%InitializeLineList()
    call this%BlockB%InitializeLineList()
    call this%CurrentBlock%InitializeLineList()
    call this%PreviousBlock%InitializeLineList()
    !
    if (.not. associated(this%fileobj)) then
      call ConstructFileType(newFile)
      this%fileobj => newFile
    endif
    !
    if (.not. associated(this%PkgObsWriter)) then
      allocate(this%PkgObsWriter)
    endif
    !
    return
  end subroutine AllocatePointers

  subroutine WriteOptions(this)
    ! Override if needed
    implicit none
    class(PackageWriterType) :: this
    ! local
    integer :: i, iu
    character(len=200) :: ermsg
    character(len=MAXCHARLEN) :: auxline
    ! formats
    5  format()
    20 format(2x,a,2x,a,2x,a)
    50 format(2x,a)
60     format(a)
70 format(a,1x,i0,1x,a)
    !
    iu = this%fileobj%IUnit
    ! Write BEGIN Options
    write(iu,5)
    write(iu,60)'BEGIN Options'
    !
!    ! Aux variable names
!    if (this%NAux > 0) then
!      do i=1,this%NAux
!        write(iu,20) 'AUXILIARY', trim(this%Aux(i))
!      enddo
!    endif
    ! Aux variable names
    if (this%NAux > 0) then
      auxline = 'AUXILIARY'
      do i=1,this%NAux
        auxline = trim(auxline) // ' ' // trim(this%Aux(i))
      enddo
      write(iu,20)trim(auxline)
    endif
    !
    if (echo) then
      if (this%ipr==1) then
        this%PrintInput = .true.
      endif
    endif
    !
    if (associated(this%ICbc)) then
      if (this%ICbc > 0) then
        this%SaveFlows = .true.
      elseif (this%ICbc < 0) then
        this%PrintFlows = .true.
      endif
    endif
    !
    ! Print_Input
    if (echo) then
      if (this%PrintInput) write(iu,50)'PRINT_INPUT'
    endif
    !
    ! Print_Flows
    if (this%PrintFlows) write(iu,50)'PRINT_FLOWS'
    !
    ! Save_Flows
    if (this%SaveFlows) write(iu,50)'SAVE_FLOWS'
    !!
    !! Newton
    !if (.not. this%Newton) write(iu,50)'NO_NEWTON'
    !
    if (this%ObsActive) then
      if (this%PkgObsWriter%fileobj%FName == '') then
        write(ermsg,70)'Programmer error: OBS6 file for package opened on unit', &
                       this%fileobj%IUnit, 'is unnamed.'
        call store_error(ermsg)
        call ustop()
      endif
      write(iu,20)'OBS6', 'FILEIN', trim(this%PkgObsWriter%fileobj%FName)
    endif
    !
    if (this%TsFile /= '') then
      write(iu,20)'TS6 FILEIN',trim(this%TsFile)
    endif
    !
    ! Write END Options
    write(iu,60)'END Options'
    !
    return
  end subroutine WriteOptions

  subroutine WriteDimensions(this)
    implicit none
    class(PackageWriterType) :: this
    ! local
    integer :: iu
    ! formats
    5  format()
    30 format(2x,a,2x,i0)
    60 format(a)
    !
    if (.not. this%NeedDimensionsBlock) return
    !
    iu = this%fileobj%IUnit
    ! Write BEGIN Dimensions
    write(iu,5)
    write(iu,60)'BEGIN Dimensions'
    !
    write(iu,30)'MAXBOUND',this%MaxActiveBnd
    !
    ! Write END Dimensions
    write(iu,60)'END Dimensions'
    !
    return
  end subroutine WriteDimensions

  subroutine WriteBlockIfNeeded(this, kper, option, forceWrite)
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    integer, intent(in) ::  kper
    character(len=*), optional :: option
    logical, optional :: forceWrite
    ! local
    integer :: kcurrent
    logical :: needToWrite, forceWriteLocal
    integer :: i, iu
    character(len=200) :: line
    ! formats
    5  format()
    10 format(a)
    20 format('BEGIN PERIOD ',i0)
    25 format('BEGIN PERIOD ',i0,2x,a)
    30 format('END PERIOD')
    !
    ! Determine if a block is needed for current stress period,
    ! and if it is, write it.
    forceWriteLocal = .false.
    if (present(forceWrite)) then
      forceWriteLocal = forceWrite
    endif
    !
    if (forceWriteLocal) then
      needToWrite = .true.
      kcurrent = this%CurrentBlock%CountLines()
    else
      needToWrite = .false.
      kcurrent = this%CurrentBlock%CountLines()
      if (kper==1) then
        if (kcurrent>0) needToWrite = .true.
      else
        ! kper > 1
        needToWrite = .not. same_lines(this%BlockA, this%BlockB)
      endif
    endif
    !
    if (needToWrite) then
      ! write block of stress-period input for MF6
      iu = this%fileobj%IUnit
      write(iu,5)
      if (present(option)) then
        write(iu,25)kper,trim(option)
      else
        write(iu,20)kper
      endif
      do i=1,kcurrent
        call this%CurrentBlock%GetLine(i, line)
        write(iu,10)trim(line)
      enddo
      write(iu,30)
    endif
    !
    return
  end subroutine WriteBlockIfNeeded

  subroutine WriteStressPeriodListData(this, lineList)
    ! Override if needed
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    type(LineListType), pointer :: lineList
    ! local
    integer :: i, ii, j, k, kold, n
    character(len=MAXCHARLEN) :: line, msg
    ! format
    10 format(a,' boundary removed at (',i0,',',i0,',',i0,')')
    !
    ! Write stress-period data for MF6
    do ii=1,this%NBndPeriod
      kold = nint(this%rlist(1,ii))
      k = this%Layptr(kold)
      i = nint(this%rlist(2,ii))
      j = nint(this%rlist(3,ii))
      if (associated(masteridomain)) then
        if (masteridomain(j,i,k)==0) then
          if (verbose) then
            msg = ''
            write(msg,10)trim(this%PkgType),k,i,j
!            write(*,*)trim(msg)
            call store_note(msg)
          endif
          cycle
        endif
      else
        if (IBOUND(j,i,kold)==0) then
          if (verbose) then
            msg = ''
            write(msg,10)trim(this%PkgType),k,i,j
!            write(*,*)trim(msg)
            call store_note(msg)
          endif
          cycle
        endif
      endif
      line = ' '
      write(line,this%fmat)k, i, j, (this%rlist(n,ii),n=this%nstart,this%nstop)
      call lineList%AddLine(line)
    enddo
    !
    return
  end subroutine WriteStressPeriodListData

  subroutine WriteStressPeriodArrayData(this, array, lineList, convertToFlux)
    implicit none
    ! dummy
    class(PackageWriterType)    :: this
    double precision, dimension(NCOL, NROW) :: array
    type(LineListType), pointer :: lineList
    logical                     :: convertToFlux
    ! local
    integer :: fpl, i, j, jstart, jend, n, nrem, nrep
    character(len=MAXCHARLEN) :: line
    !
    ! Write array data to lineList
    do i=1,NROW
      do j=1,NCOL
        ! If specified, convert flow rates to flux rates by dividing by cell area
        if (convertToFlux) then
          this%work(j,i) = array(j,i)/this%CellArea(j,i)
        else
          this%work(j,i) = array(j,i)
        endif
      enddo
    enddo
    !
    fpl = this%FieldsPerLine
    nrep = NCOL / fpl
    nrem = NCOL - (nrep * this%FieldsPerLine)
    do i=1,NROW
      jstart = 1
      do n=1,nrep
        line = ''
        jend = jstart + fpl - 1
        write(line,this%fmat)(this%work(j,i),j=jstart,jend)
        call lineList%AddLine(line)
        jstart = jend + 1
      enddo
      if (nrem > 0) then
        jend = jstart + nrem - 1
        line = ''
        write(line,this%fmat)(this%work(j,i),j=jstart,jend)
        call lineList%AddLine(line)
      endif
    enddo
    !
    return
  end subroutine WriteStressPeriodArrayData

  subroutine WriteArrayDataIfNeeded(this, iu, kper, array, currentBlock, previousBlock, &
                                    writing, convertToFlux, label)
    implicit none
    ! dummy
    class(PackageWriterType), intent(inout)    :: this
    integer, intent(in)                        :: iu, kper
    double precision, dimension(NCOL, NROW), intent(inout) :: array
    type(LineListType), pointer, intent(inout) :: currentBlock, previousBlock
    logical, intent(in)                        :: convertToFlux
    logical, intent(inout)                     :: writing
    character(len=*), intent(in)               :: label
    ! local
    integer :: i, iprnr, n
    logical :: constant, needToWrite
    double precision    :: rval
    character(len=MAXCHARLEN) :: line
    ! formats
    10 format(a)
    20 format(2x,a)
    30 format(/,'BEGIN PERIOD ',i0)
    60 format(4x,a,2x,g15.8)
    80 format(4x,a,2x,i0)
    !
    if (echo) then
      iprnr = 12
    else
      iprnr = -12
    endif
    ! Determine if max ET rate array differs from previous stress period,
    ! and write array if needed.
    call this%WriteStressPeriodArrayData(array, currentBlock, convertToFlux)
    needToWrite = .not. same_lines(currentBlock, previousBlock)
    if (needToWrite) then
      call ConstantReal2D(NCOL, NROW, this%work, constant, rval)
      ! If kper=1 and rval is zero, no need to write
      if (kper == 1 .and. constant .and. rval == 0.0) needToWrite = .false.
      if (needToWrite) then
        if (.not. writing) then
          write(iu,30)kper
          writing = .true.
        endif
        write(iu,20) trim(label)
        if (constant) then
          write(iu,60)'CONSTANT', rval
        else
          ! write array control line and array to input file
          write(iu,80)'INTERNAL  FACTOR  1.0  IPRN ',iprnr
          n = currentBlock%CountLines()
          do i=1,n
            call currentBlock%GetLine(i, line)
            write(iu,10)trim(line)
          enddo
        endif
      endif
    endif
    !
    return
  end subroutine WriteArrayDataIfNeeded

  function CastAsPackageWriterClass(obj) result(res)
    ! Cast an unlimited polymorphic object as class(PackageWriterType)
    implicit none
    class(*), pointer, intent(inout) :: obj
    class(PackageWriterType), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    class is (PackageWriterType)
      res => obj
    end select
    return
  end function CastAsPackageWriterClass

  subroutine GetBudgetText(this, text)
    ! Return text that will be associated with this
    ! package in MF6 volume budget tables.
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    character(len=*), intent(out) :: text
    !
    if (this%PackageName /= '') then
      text = this%PackageName
    elseif (this%DefaultBudgetText /= '') then
      text = this%DefaultBudgetText
    else
      text = '(undefined)'
    endif
    !
    return
  end subroutine GetBudgetText

  subroutine AddMover(this, mover)
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    type(MoverType), pointer :: mover
    !
    if (mover%ProvModelName == mover%RecModelName) then
      call AddMoverToList(this%ModelMovers, mover)
    else
      call AddMoverToList(SimMovers, mover)
    endif
    !
    return
  end subroutine AddMover

  subroutine BuildModelMovers(this)
    ! Most packages do not support water movers, so this procedure does nothing.
    ! Packages that support water movers must override this procedure.
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    !
    return
  end subroutine BuildModelMovers

  subroutine BuildSimMovers(this)
    ! Most packages do not support water movers, so this procedure does nothing.
    ! Packages that support water movers must override this procedure.
    implicit none
    ! dummy
    class(PackageWriterType) :: this
    !
    return
  end subroutine BuildSimMovers

  logical function SameType(this, pkgwriter)
    ! dummy
    class(PackageWriterType) :: this
    class(PackageWriterType) :: pkgwriter
    !
    character(len=LENCTYPE) :: thistype, othertype
    !
    SameType = .false.
    thistype = this%MyType()
    othertype = pkgwriter%MyType()
    if (thistype == othertype) then
      SameType = .true.
    endif
    !
    return
  end function SameType

end module PackageWriterModule


text:
      subroutine ggridp(ipch_asv,ipcv_asv,itgenf,jpdblo,jpfc_asv,
     $ jptffl,narx_asv,ndat0s,ndb_asv,noutpt,ntid_asv,ntpr_asv,
     $ nttyo,q500fl,uakey)
c
c     This suboutine gets the necessary dimensioning parameters for the
c     temperature grid on which thermodynamic data are represented.
c     These parameters are presently embedded in the data file title,
c     in order to maintain maximum forward and backward compatibility
c     of data files between Versions 7 and 8 of EQ3/6. There is
c     compatibility if the temperature grid matches the old "standard"
c     form and if there are no data grids for volume and enthalpy data.
c
c     This suboutine is called by:
c
c       EQPT/eqpt.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       ndat0s = unit number of the stripped DATA0 file
c       uakey  = string identifying the type of aqueous species
c                activity coefficient model ("SEDH" or "Pitzer")
c
c
c     Principal output:
c
c       ipch_asv = the maximum order of pressure corrections to the
c                    enthalpy functions; ipch_asv = min(1,ipch)
c       ipcv_asv = the maximum order of pressure corrections to the
c       itgenf   = temperature grid enforcement parameter:
c                    -1 = do nothing
c                     0 = warn if any sparsely filled ranges in
c                         "log K" temperature grid data
c                     1 = declare error condition if any sparsely
c                           filled ranges in "log K" temperature
c                           grid data
c       jpdblo   = Pitzer data block organization flag; -1 =
c                    "classical", 0 = new
c                    volume functions; ipcv_asv = min(1,ipcv)
c       jpfc_asv = the number of terms in the temperature function
c                    used to represent Pitzer interaction parameters;
c                    must be at least 1
c       jptffl   = integer flag denoting the temperature function
c                    used to represent Pitzer interaction parameters;
c                    -1 = classical truncated Taylor's series,
c                    0 = LLNL5TERM, 1 = Greenberg and Moller (1989)
c                    eight-term equation
c       narx_asv = the maximum number of points in any range of the
c                    temperature grid
c       ntid_asv = the number of lines in the data file title
c       ntpr_asv = the number of temperature ranges on the temperature
c                    grid (ntpr_asv = ntprt)
c       ndb_asv  = the maximum number of distinct points on the
c                    temperature grid;
c                    ndb_asv = ntpr_asv*(narx_asv - 1) + 1
c       q500fl   = logical flag; if .true., treat instances of "500."
c                    in "log K" temperature grids as a "no data"
c                    condition"
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ndat0s,noutpt,nttyo
c
      integer ipch_asv,ipcv_asv,jpfc_asv,narx_asv,ndb_asv,ntid_asv,
     $ ntpr_asv
c
      integer itgenf,jpdblo,jptffl
c
      character(len=8) uakey
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,j,j2,k,n,narxti,narx_x,ntpr,ntpri
c
      integer ilnobl
c
      logical q500fl,qnarxr,qrderr
c
      character(len=80), dimension(:), allocatable :: ux80ar
c
      character(len=80) ulbufa,ulbufb,uline
      character(len=16) ustr16,uterm
      character(len=8) ustr
c
c-----------------------------------------------------------------------
c
      data uterm  /'+---------------'/
c
c-----------------------------------------------------------------------
c
c     Determine the number of lines in the data file title.
c     Count that first terminator line as the last line of the title.
c
      n = 0
  100 read (ndat0s,1000,end=990,err=995) uline
 1000 format(a)
      n = n + 1
      if (uline(1:16) .ne. uterm(1:16)) go to 100
c
      ntid_asv = n
      rewind(ndat0s)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Allocate an array to hold the lines in the title.
c
      ALLOCATE(ux80ar(ntid_asv))
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Read the title into the holding array.
c
      do n = 1,ntid_asv
        read (ndat0s,1000,end=990,err=995) ux80ar(n)
      enddo
c
      rewind(ndat0s)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Set parameters for the temperature grid used to represent log K
c     values and standard state enthalpy and volume functions.
c
c     Values corresponding to the "classical" EQ3/6 temperature grid
c     will be used as defaults. These are intended to support the use
c     of version 7 format data files by the version 8 and higher
c     software.
c
c       ntprt       = the number of temperature ranges
c       narxt(ntpr) = the number of points or coefficients in the
c                       ntpr-th temperature range
c
c     The old standard grid corresponds to:
c
c       ntprt = 2
c       narxt(1) = 4  (0, 25, 60 100 C)
c       narxt(2) = 5  (100, 150, 200, 250, 300 C)
c       ipch = -1   (no enthalpy function grids are present)
c       ipcv = -1   (no volume function grids are present)
c
c     Note that the last point of the ntpr-th temperature range
c     is also the first point of the following range, if any.
c     This insures continuity.
c
c     Hence the matching default dimensioning is as follows. Note
c     that if ipch < 1, ipch_asv must be set to one, as an array
c     can't have a null or negative dimension. The same applies to
c     ipcv and ipcv_asv.
c
      ntpr_asv = 2
      narx_asv = 5
      ipch_asv = 1
      ipcv_asv = 1
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the number of temperature ranges (ntpr_asv).
c
      do n = 1,ntid_asv
        ulbufa = ' '
        ulbufb = ' '
        j = 0
c
c       Check for a keystring.
c
        i = index(ux80ar(n),'NO. OF TEMPERATURE RANGES')
        if (i .gt. 0) then
          j = i + 25
        else
          i = index(ux80ar(n),'NUMBER OF TEMPERATURE RANGES')
          if (i .gt. 0) then
            j = i + 28
          else
            i = index(ux80ar(n),'NTPRT')
            if (i .gt. 0) j = i + 5
          endif
        endif
c
        if (j .gt. 0) then
c
c         Check for an equal sign following the keystring.
c
          ulbufa = ux80ar(n)(j:80)
          call lejust(ulbufa)
          i = index(ulbufa,'=')
          j = 2
          if (i .ne. 1) j = 0
        endif
c
        if (j .gt. 0) then
c
c         Extract a number matching the keystring.
c
          ulbufb = ulbufa(2:80)
          call lejust(ulbufb)
          ustr = ulbufb(1:8)
          call chrint(ntpr_asv,nttyo,qrderr,ustr)
          if (qrderr) go to 997
          go to 200
        endif
      enddo
c
  200 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the maximum number of points in any temperature range
c     (narx_asv). Note that this is done by finding the number of
c     points in each range and determining the maximum. However, an
c     array to hold the actual values has not yet been allocated.
c     The actual values will be redetermined and stored by rdwttl.f.
c
      qnarxr = .false.
      narx_x = 0
      narxti = 0
      do ntpr = 1,ntpr_asv
        do n = 1,ntid_asv
          ulbufa = ' '
          ulbufb = ' '
          j = 0
c
c         Check for a keystring.
c
          i = index(ux80ar(n),'NO. OF POINTS IN RANGE')
          if (i .gt. 0) then
            j = i + 22
          else
             i = index(ux80ar(n),'NUMBER OF POINTS IN RANGE')
            if (i .gt. 0) then
              j = i + 25
            else
              i = index(ux80ar(n),'NARXT')
              if (i .gt. 0) j = i + 5
            endif
          endif
c
          if (j .gt. 0) then
c
c           Extract a number acting as a subscript to the keystring.
c
            ulbufb = ux80ar(n)(j:80)
            call lejust(ulbufb)
            j = index(ulbufb,' ')
            j = min(j,9)
            k = j - 1
            if (k .gt. 0) then
              ustr = ulbufb(1:k)
              call chrint(ntpri,nttyo,qrderr,ustr)
              if (qrderr) go to 997
              if (ntpri .ne. ntpr) j = 0
            else
              j = 0
            endif
          endif
c
          if (j .gt. 0) then
c
c           Check for an equal sign following the keystring.
c
            ulbufa = ulbufb(j:80)
            call lejust(ulbufa)
            i = index(ulbufa,'=')
            j = 2
            if (i .ne. 1) j = 0
          endif
c
          if (j .gt. 0) then
c
c           Extract a number matching the keystring.
c
            ulbufb = ulbufa(j:80)
            call lejust(ulbufb)
            ustr = ulbufb(1:8)
            call chrint(narxti,nttyo,qrderr,ustr)
            if (qrderr) go to 997
            qnarxr = .true.
            narx_x = max(narx_x,narxti)
            go to 220
          endif
        enddo
  220   continue
c
      enddo
c
c     If a value was read for any temperature range, use it. But
c     do not use a number less than 1, as an array can not have a
c     null or negative dimension.
c
      if (qnarxr) narx_asv = narx_x
      narx_asv = max(1,narx_asv)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Calculate the maximum possible number of points on the temperature
c     grid for the current grid parameters.
c
      ndb_asv = ntpr_asv*(narx_asv - 1) + 1
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Set flag for dealing with sparsely filled ranges in the
c     "log K" temperature grids.
c       itgenf:
c         -1 = do nothing
c          0 = generate warnings for sparsely filled ranges
c          1 = generate errors for sparsely filled ranges
c
      itgenf = -1
      do n = 1,ntid_asv
        ulbufa = ' '
        j = 0
c
c       Check for a keystring.
c
        i = index(ux80ar(n),'SPARSE GRID RANGE CONDITION')
        if (i .gt. 0) j = i + 27
c
        if (j .gt. 0) then
c
c         Check for an equal sign following the keystring.
c
          ulbufa = ux80ar(n)(j:80)
          call lejust(ulbufa)
          i = index(ulbufa,'=')
          j = 2
          if (i .ne. 1) j = 0
        endif
c
        if (j .gt. 0) then
c
c         Extract a string input matching the keystring.
c
          ustr16 = ulbufa(2:17)
          call lejust(ustr16)
          k = index(ustr16,'IGNORE')
          if (k .gt. 0) then
            itgenf = -1
          else
            k = index(ustr16,'WARN')
            if (k .gt. 0) then
              itgenf = 0
            else
              k = index(ustr16,'ERROR')
              if (k .gt. 0) then
                itgenf = 1
              endif
            endif
          endif
          go to 300
        endif
      enddo
  300 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Set flag for dealing with values of "500." on the "log K"
c     temperature grids.
c       q500fl:
c         .false. = do nothing
c          .true. = treat values of 500. as "no data"
c
      q500fl = .true.
      do n = 1,ntid_asv
        ulbufa = ' '
        j = 0
c
c       Check for a keystring.
c
        i = index(ux80ar(n),'INTERPRET 500 AS NO DATA')
        if (i .gt. 0) j = i + 24
c
        if (j .gt. 0) then
c
c         Check for an equal sign following the keystring.
c
          ulbufa = ux80ar(n)(j:80)
          call lejust(ulbufa)
          i = index(ulbufa,'=')
          j = 2
          if (i .ne. 1) j = 0
        endif
c
        if (j .gt. 0) then
c
c         Extract a string input matching the keystring.
c
          ustr16 = ulbufa(2:17)
          call lejust(ustr16)
          k = index(ustr16,'NO')
          if (k .le. 0) k = index(ustr16,'FALSE')
          if (k .gt. 0) then
            q500fl = .false.
          else
            k = index(ustr16,'YES')
            if (k .le. 0) k = index(ustr16,'TRUE')
            if (k .gt. 0) then
              q500fl = .true.
            endif
          endif
          go to 310
        endif
      enddo
  310 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the dimension corresponding to the order of the pressure
c     dependence of the enthalpy functions. An order of -1 means that
c     No enthalpy functions data are present on the data file. An order
c     of 0 means that enthalpy functions data are present, but no data
c     for the pressure dependence of the same are present. However,
c     the dimensioning variable ipch_asv must be set to a value of
c     at least one, as an array can not have a null or negative
c     dimension. Thus, it is necessary here to check only for data
c     corresponding to order 2 or higher.
c
      do n = 1,ntid_asv
        ulbufa = ' '
        j = 0
c
c       Check for a keystring.
c
        i = index(ux80ar(n),'d2H/dP2')
        if (i .gt. 0) j = i + 7
c
        if (j .gt. 0) then
c
c         Check for an equal sign following the keystring.
c
          ulbufa = ux80ar(n)(j:80)
          call lejust(ulbufa)
          i = index(ulbufa,'=')
          j = 2
          if (i .ne. 1) j = 0
        endif
c
        if (j .gt. 0) then
c
c         Extract a string input matching the keystring.
c
          ustr = ulbufa(2:9)
          call lejust(ustr)
          k = index(ustr,'ON')
          if (k .le. 0) k = index(ustr,'PRESENT')
          if (k .le. 0) k = index(ustr,'ACTIVE')
          if (k .gt. 0) ipch_asv = 2
          go to 320
        endif
      enddo
  320 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the dimension corresponding to the order of the pressure
c     dependence of the volume functions. This process and the dimension
c     itself are analogous to those described above for the pressure
c     dependence of the enthalpy functions.
c
      do n = 1,ntid_asv
        ulbufa = ' '
        j = 0
c
c       Check for a keystring.
c
        i = index(ux80ar(n),'d2V/dP2')
        if (i .gt. 0) j = i + 7
c
        if (j .gt. 0) then
c
c         Check for an equal sign following the keystring.
c
          ulbufa = ux80ar(n)(j:80)
          call lejust(ulbufa)
          i = index(ulbufa,'=')
          j = 2
          if (i .ne. 1) j = 0
        endif
c
        if (j .gt. 0) then
c
c         Extract a string input matching the keystring.
c
          ustr = ulbufa(2:9)
          call lejust(ustr)
          k = index(ustr,'ON')
          if (k .le. 0) k = index(ustr,'PRESENT')
          if (k .le. 0) k = index(ustr,'ACTIVE')
          if (k .gt. 0) ipcv_asv = 2
          go to 350
        endif
      enddo
  350 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Set flag for Pitzer data block organization.
c       jpdblo:
c         = -1 Classical: two superblocks, one for pure aqueous
c              electrolytes, one for mixtures of two such electrolytes
c              having a common ion. A given theta parameter can appear
c              multiple times. Data for parameters involving electrically
c              neutral species are folded into this structure. Checks are
c              in place to sort things out.
c         =  0 New: multiple superblocks, each corresponding to a
c              charge-based type of species pair or triplet. This
c              structure is more rational. Theta parameters for example
c              have their own superblock.
c
      jpdblo = -1
      if (uakey(1:8) .eq. 'Pitzer  ') then
        do n = 1,ntid_asv
          ulbufa = ' '
          j = 0
c
c         Check for a keystring.
c
          i = index(ux80ar(n),'PITZER DATA BLOCK ORG.')
          if (i .gt. 0) j = i + 22
c
          if (j .gt. 0) then
c
c           Check for an equal sign following the keystring.
c
            ulbufa = ux80ar(n)(j:80)
            call lejust(ulbufa)
            i = index(ulbufa,'=')
            j = 2
            if (i .ne. 1) j = 0
          endif
c
          if (j .gt. 0) then
c
c           Extract a string input matching the keystring.
c
            ustr16 = ulbufa(2:17)
            call lejust(ustr16)
            k = index(ustr16,'NEW')
            if (k .le. 0) k = index(ustr16,'BYSPECIES')
            if (k .gt. 0) jpdblo = 0
            go to 400
          endif
        enddo
  400   continue
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the identification of the temperature function used to
c     represent Pitzer interaction parameters. Note that the number
c     of terms or coefficients (jpfcmx) is a dimensioning parameter
c     that may also be set by an option embedded in the data file
c     title.
c
c     Temperature function flag.
c       jptffl:
c         = -1 Classic: a Taylor's series centered at 25C and
c              truncated at second order.
c
c                a1 = value at 25C
c                a2 = first derivative at 25C
c                a3 = second derivative at 25C
c
c              This is currently required if the Pitzer data block
c              organization flag is set to "Classical" (jpdblo = -1).
c              The number of terms used (jpfcmx) is fixed at 3.
c
c         =  0 The function:
c
c                x(T) = a1 + a2/(T - Tr) + a3*ln(T/Tr)
c                            a4(T - Tr) + a5(T**2 - Tr**2)
c
c              This is also 25C centric (Tr = 25C).
c              The number of terms used (jpfcmx) may vary from
c              1 to 5. This always refers to the first jpfcmx terms.
c
c         =  1 The Greenberg-Moller (1988) combination temperature
c              function. This is not 25C centric. The number of terms
c              used (mx) is fixed at 8.
c
      jptffl = -1
      if (uakey(1:8) .eq. 'Pitzer  ') then
        do n = 1,ntid_asv
          ulbufa = ' '
          j = 0
c
c         Check for a keystring.
c
          i = index(ux80ar(n),'PITZER TEMP FUNCTION')
          if (i .gt. 0) j = i + 20
c
          if (j .gt. 0) then
c
c           Check for an equal sign following the keystring.
c
            ulbufa = ux80ar(n)(j:80)
            call lejust(ulbufa)
            i = index(ulbufa,'=')
            j = 2
            if (i .ne. 1) j = 0
          endif
c
          if (j .gt. 0) then
c
c           Extract a string input matching the keystring.
c
            ustr16 = ulbufa(2:17)
            call lejust(ustr16)
            k = index(ustr16,'NEW')
            if (k .le. 0) k = index(ustr16,'Livermore')
            if (k .le. 0) k = index(ustr16,'LIVERMORE')
            if (k .le. 0) k = index(ustr16,'LLNL')
            if (k .le. 0) k = index(ustr16,'L5TERM')
            if (k .le. 0) k = index(ustr16,'5TERM')
            if (k .gt. 0) then
              jptffl = 0
              go to 410
            endif
            k = index(ustr16,'TEQUIL')
            if (k .le. 0) k = index(ustr16,'T8TERM')
            if (k .le. 0) k = index(ustr16,'GM8TERM')
            if (k .le. 0) k = index(ustr16,'8TERM')
            if (k .gt. 0) then
              jptffl = 1
              go to 410
            endif
          endif
        enddo
  410   continue
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Get the number of terms or coefficients (jpfc_asv) in the
c     temperature function used to represent Pitzer parameters.
c
      jpfc_asv = 1
      if (uakey(1:8) .eq. 'Pitzer  ') then
c
c       Three is the default value for a Pitzer model.
c
        jpfc_asv = 3
        do n = 1,ntid_asv
          ulbufa = ' '
          j = 0
c
c         Check for a keystring.
c
          i = index(ux80ar(n),'NO. OF PITZER TEMP FUNC TERMS')
          if (i .gt. 0) j = i + 29
c
          if (j .gt. 0) then
c
c           Check for an equal sign following the keystring.
c
            ulbufa = ux80ar(n)(j:80)
            call lejust(ulbufa)
            i = index(ulbufa,'=')
            j = 2
            if (i .ne. 1) j = 0
          endif
c
          if (j .gt. 0) then
c
c           Extract a string input matching the keystring.
c
            ustr = ulbufa(2:9)
            call lejust(ustr)
            read (ustr,'(i8)') jpfc_asv
            go to 420
          endif
        enddo
  420   continue
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Deallocate the array holding the data file title.
c
      DEALLOCATE(ux80ar)
      go to 999
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Write messages for errors associated with reading the DATA0
c     file.
c
  990 write (noutpt,2000)
      write (nttyo,2000)
 2000 format(/' * Error - (EQPT/ggridp) Unexpectedly encountered',
     $ /7x,'end-of-file while reading the title of the DATA0 file.')
      if (n .le. 1) then
        write (noutpt,2010)
        write (nttyo,2010)
 2010   format(/7x,'This occurred while trying to read the first line',
     $  /7x,'of the title.')
      else
        ulbufa = ux80ar(n - 1)
        call lejust(ulbufa)
        j2 = ilnobl(ulbufa)
        j2 = min(j2,70)
        write (noutpt,2020) ulbufa(1:j2)
        write (nttyo,2020) ulbufa(1:j2)
 2020   format(/7x,'This occurred while trying to read the line',
     $  ' following',/7x,'the line:',/7x,'"',a,'"')
      endif
      stop
c
  995 write (noutpt,2030)
      write (nttyo,2030)
 2030 format(/' * Error - (EQPT/ggridp) Encountered a read format',
     $ /7x,'error while reading the DATA0 file.')
      if (n .le. 1) then
        write (noutpt,2010)
        write (nttyo,2010)
      else
        ulbufa = ux80ar(n - 1)
        call lejust(ulbufa)
        j2 = ilnobl(ulbufa)
        j2 = min(j2,70)
        write (noutpt,2020) ulbufa(1:j2)
      endif
      stop
c
  997 ulbufa = ux80ar(n)
      call lejust(ulbufa)
      j2 = ilnobl(ulbufa)
      j2 = min(j2,70)
      write (noutpt,2050) ulbufa(1:j2)
      write (nttyo,2050) ulbufa(1:j2)
 2050 format(/' * Error - (EQPT/ggridp) Encountered a read format',
     $ /7x,'error while reading data embedded in the title of the',
     $ /7x,'DATA0 file. This occurred while attempting to process',
     $ ' the line:',/7x,'"',a,'"')
      stop
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


text:
   MODULE test2008_41_A
     REAL :: X=5.0
   END MODULE test2008_41_A
   MODULE test2008_41_B
     USE test2008_41_A
     PRIVATE :: X               !  X cannot be accessed through module B
     REAL :: C=80, D=50
   END MODULE test2008_41_B
   PROGRAM TEST
     INTEGER :: TX=7
     CALL SUB
     CONTAINS

     SUBROUTINE SUB
     USE test2008_41_B, ONLY : C
     USE test2008_41_B, T1 => C
     USE test2008_41_B, TX => C             !  C is given another local name
     USE test2008_41_A
     PRINT *, TX                !  Value written is 80 because use-associated
                                !  entity overrides host entity
     END SUBROUTINE
   END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab