text:
	SUBROUTINE MA_CGMM ( string, mxormn, ier )
C************************************************************************
C* MA_CGMM                                                              *
C*                                                                      *
C* This subroutine gets the maximum or minimum temperature and stores	*
C* it into the interface array in common /rintfv/.			*
C*                                                                      *
C* MA_CGMM  ( STRING, MXORMN, IER )					*
C*                                                                      *
C* Input parameters:                                                    *
C*	RIVALS(IRMNTH)  REAL		Month
C*      STRING          CHAR*		Character string with min/max.	*
C*	MXORMN		INTEGER		Flag to get max or min temp.	*
C*                                                                      *
C* Output parameters:                                                   *
C*	RIVALS(IRMXTM)	REAL		Maximum temperature (degrees C)	*
C*	RIVALS(IRMITM)  REAL		Minimum temperature (degrees C)	*
C*      IER             INTEGER         Return code                     *
C*                                        0 = Normal return             *
C**                                                                     *
C* Log:                                                                 *
C* F. J. Yen/NCEP	 4/01	Created from CG_PRTM.			*
C* D. Kidwell/NCEP	 3/02	Converted temps to degrees C            *
C************************************************************************
	INCLUDE		'macmn.cmn'
C*
	CHARACTER*(*)	string
C------------------------------------------------------------------------
	IF ( mxormn .eq. 1 ) THEN
C
C*	    Get maximum temperature.
C
	    CALL ST_INTG ( string, ist1, ier )
	    rivals(irmxtm) = FLOAT( ist1 )
	    IF ( rivals(irmnth) .gt. 3 .and.
     +		    rivals(irmnth) .lt. 11 ) THEN
C
C*		Check for 100+ temperatures.
C
		IF ( ist1 .lt. 20 ) THEN
		    rivals(irmxtm) = FLOAT( ist1 ) + 100.
		END IF
	      ELSE
C
C*		Check for sub-zero temperatures.
C
		IF ( ist1 .ge. 90 ) THEN
		    rivals(irmxtm) = FLOAT( ist1 ) - 100.
		END IF
	    END IF
	    rivals(irmxtm) = PR_TMFC ( rivals(irmxtm) )
	  ELSE
C
C*	    Get minimum temperature.
C
	    CALL ST_INTG ( string, ist1, ier )
	    rivals(irmitm) = FLOAT( ist1 )
	    IF ( rivals(irmnth) .lt. 5 .or.
     +		        rivals(irmnth) .gt. 10 ) THEN
C
C*		Check for sub-zero temperatures.
C
		IF ( ist1 .gt. 70 .and. rivals(irtmpf) .lt. ist1 ) THEN
		    rivals(irmitm) = FLOAT( ist1 ) - 100.
		END IF
	    END IF
	    rivals(irmitm) = PR_TMFC ( rivals(irmitm) )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program i
  real :: a, b, c, d(10), e(1:11, 23)
end program i

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
	SUBROUTINE GDNARW ( ictyp, line, loci, info, iret )
C************************************************************************
C* GDNARW								*
C*									*
C*									*
C**									*
C* Log:									*
C* S. Jacobs/SSAI	12/91						*
C* T. Lee/SAIC		01/06	Copied from GPNARW			*
C* C. Bailey/HPC	10/06	Changed call to IN_LINE			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	PARAMETER	( HANGLE1 = DTR * 30.0 )
	PARAMETER	( HANGLE2 = DTR * 10.0 )
C*
	CHARACTER*(*)	line, loci, info
C*
	CHARACTER	sys*1, arwpts(6)*24
	REAL		tmpt(2), xarw(2), yarw(2), xfarw(8), yfarw(8)
	REAL		smooth, linfltr
	LOGICAL		scflag
C------------------------------------------------------------------------
C*	Get plot bounds in Normal coordinates.
C
	CALL GQBND ( 'N', xleft, ybot, xrght, ytop, ier )
C
C*	Parse the fill input.
C
	CALL IN_LINE ( line, 1, 1, icolor, itype, iwidth, ilabel,
     +                 smooth, linfltr, scflag, ier )
C
	IF  ( icolor .ne. 0 )  THEN
C
C*	    Set coordinate system to be used.
C
	    sys = 'N'
	    IF  ( loci(1:1) .eq. '@' )  THEN
		sys = 'G'
		loci = loci(2:)
	    ELSE IF  ( loci(1:1) .eq. '#' )  THEN
		sys = 'M'
		loci = loci(2:)
	    END IF
C
C*	    Get the point(s) for the plot.
C
	    CALL ST_RLST ( loci, ';', 0., 2, tmpt, num2, ier )
C
C*	    Separate the shape information.
C
	    CALL ST_CLST ( info, '/', ' ', 6, arwpts, num, ier )
C
C*	    Get the line type.
C
	    CALL ST_NUMB ( arwpts(1), ltype, ier )
	    IF  ( ltype .eq. IMISSD )  ltype = 1
C
C*	    Get the arrow type.
C
	    CALL ST_NUMB ( arwpts(2), iatp, ier )
	    IF  ( iatp .eq. IMISSD )  iatp = 1
C
C*	    Get the arrow size.
C
	    CALL ST_CRNM ( arwpts(3), asize, ier )
	    IF  ( asize .eq. RMISSD )  asize = 1.
C
C*	    Get the head size.
C
	    CALL ST_CRNM ( arwpts(4), hsize, ier )
	    IF  ( hsize .eq. RMISSD )  hsize = 1.
C
C*	    Get the head angle multiplier.
C
	    CALL ST_CRNM ( arwpts(5), hdang, ier )
	    IF  ( hdang .eq. RMISSD )  hdang = 1.
C
C*	    Get the rotation angle.
C
	    CALL ST_CRNM ( arwpts(6), rot, ier )
	    IF  ( rot .eq. RMISSD )  rot = 0.
C
C*	    If using Normal coords, correct for PS device.
C
	    IF  ( sys .eq. 'N' )  THEN
		tmpt(1) = (xrght-xleft) * tmpt(1)
		tmpt(2) = (ytop -ybot ) * tmpt(2)
	    ENDIF
C
C*	    Set attributes.
C
	    CALL GSLINE ( ltype, 0, iwidth, 0, ier )
	    CALL GSCOLR ( icolor, ier )
C
C*	    Find arrow length, head size and head angles.
C
	    arrwln = asize * .1
	    hsize  = hsize * .025
	    rdir   = rot * DTR
	    rang1  = rdir - ( HANGLE1 * hdang )
	    rang2  = rdir + ( HANGLE1 * hdang )
	    rang3  = rdir - ( HANGLE2 * hdang )
	    rang4  = rdir + ( HANGLE2 * hdang )
C
C*	    Choose arrow type. Type 1 is similar to the original
C*	    GEMPAK arrow. Type 2 is a solid arrow, which can be
C*	    filled or open.
C
	    IF  ( iatp .eq. 2 )  THEN
C
C*		Find the 7 points that make a filled arrow.
C
		xfarw(1) = tmpt(1)
		yfarw(1) = tmpt(2)
		xfarw(2) = ( hsize * COS(rang1) ) + tmpt(1)
		yfarw(2) = ( hsize * SIN(rang1) ) + tmpt(2)
		xfarw(3) = ( hsize * .75 * COS(rang3) ) + tmpt(1)
		yfarw(3) = ( hsize * .75 * SIN(rang3) ) + tmpt(2)
		xfarw(4) = ( arrwln * COS(rdir) ) + xfarw(3)
		yfarw(4) = ( arrwln * SIN(rdir) ) + yfarw(3)
		xfarw(7) = ( hsize * COS(rang2) ) + tmpt(1)
		yfarw(7) = ( hsize * SIN(rang2) ) + tmpt(2)
		xfarw(6) = ( hsize * .75 * COS(rang4) ) + tmpt(1)
		yfarw(6) = ( hsize * .75 * SIN(rang4) ) + tmpt(2)
		xfarw(5) = ( arrwln * COS(rdir) ) + xfarw(6)
		yfarw(5) = ( arrwln * SIN(rdir) ) + yfarw(6)
                IF ( ictyp .eq. 1 ) THEN
                    xfarw (8) = xfarw (1)
                    yfarw (8) = yfarw (1)
                    CALL GLINE ( sys, 8, xfarw, yfarw, ier )
                ELSE
                    CALL GFILL ( sys, 7, xfarw, yfarw, ier )
                END IF
C
	    ELSE
C
C*		Draw the arrow body.
C
		xarw(1) = tmpt(1)
		yarw(1) = tmpt(2)
		xarw(2) = ( arrwln * COS(rdir) ) + tmpt(1)
		yarw(2) = ( arrwln * SIN(rdir) ) + tmpt(2)
		CALL GLINE ( sys, 2, xarw, yarw, ier )
C
C*		Draw the arrow head.
C
		CALL GSLINE ( 1, 0, iwidth, 0, ier )
		xarw(2) = tmpt(1) + ( hsize * COS(rang1) )
		yarw(2) = tmpt(2) + ( hsize * SIN(rang1) )
		CALL GLINE ( sys, 2, xarw, yarw, ier )
C
		xarw(2) = tmpt(1) + ( hsize * COS(rang2) )
		yarw(2) = tmpt(2) + ( hsize * SIN(rang2) )
		CALL GLINE ( sys, 2, xarw, yarw, ier )
C
	    END IF
C
	END IF
C*
	RETURN
	END

ext: f


text:
! RUN: %S/test_errors.sh %s %t %f18 -fopenacc

! Check OpenACC clause validity for the following construct and directive:
!   2.6.5 Data
!   2.5.1 Parallel
!   2.5.2 Serial
!   2.5.3 Kernels
!   2.9 Loop
!   2.12 Atomic
!   2.13 Declare
!   2.14.3 Set
!   2.14.4 Update
!   2.15.1 Routine
!   2.10 Cache
!   2.11 Parallel Loop
!   2.11 Kernels Loop
!   2.11 Serial Loop
!   2.14.3 Set
!   2.16.13 Wait

program openacc_clause_validity

  implicit none

  type atype
    real(8), dimension(10) :: arr
    real(8) :: s
  end type atype

  integer :: i, j, b, gang_size, vector_size, worker_size
  integer, parameter :: N = 256
  integer, dimension(N) :: c
  logical, dimension(N) :: d, e
  integer :: async1
  integer :: wait1, wait2
  real :: reduction_r
  logical :: reduction_l
  real(8), dimension(N, N) :: aa, bb, cc
  logical :: ifCondition = .TRUE.
  type(atype) :: t
  type(atype), dimension(10) :: ta

  !ERROR: At least one clause is required on the DECLARE directive
  !$acc declare
  real(8), dimension(N) :: a, f, g, h

  !$acc init
  !$acc init if(.TRUE.)
  !$acc init if(ifCondition)
  !$acc init device_num(1)
  !$acc init device_num(i)
  !$acc init device_type(i)
  !$acc init device_type(2, i, j)
  !$acc init device_num(i) device_type(i, j) if(ifCondition)

  !$acc parallel
  !ERROR: Directive INIT may not be called within a compute region
  !$acc init
  !$acc end parallel

  !$acc serial
  !ERROR: Directive INIT may not be called within a compute region
  !$acc init
  !$acc end serial

  !$acc kernels
  !ERROR: Directive INIT may not be called within a compute region
  !$acc init
  !$acc end kernels

  !$acc parallel
  !$acc loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc serial
  !$acc loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do
  !$acc end serial

  !$acc kernels
  !$acc loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do
  !$acc end kernels

  !$acc parallel loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do

  !$acc serial loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do

  !$acc kernels loop
  do i = 1, N
    !ERROR: Directive INIT may not be called within a compute region
    !$acc init
    a(i) = 3.14
  end do

  !$acc parallel
  !ERROR: Directive SHUTDOWN may not be called within a compute region
  !$acc shutdown
  !$acc end parallel

  !$acc serial
  !ERROR: Directive SHUTDOWN may not be called within a compute region
  !$acc shutdown
  !$acc end serial

  !$acc kernels
  !ERROR: Directive SHUTDOWN may not be called within a compute region
  !$acc shutdown
  !$acc end kernels

  !$acc parallel
  !$acc loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc serial
  !$acc loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do
  !$acc end serial

  !$acc kernels
  !$acc loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do
  !$acc end kernels

  !$acc parallel loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do

  !$acc serial loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do

  !$acc kernels loop
  do i = 1, N
    !ERROR: Directive SHUTDOWN may not be called within a compute region
    !$acc shutdown
    a(i) = 3.14
  end do

  !$acc parallel
  !ERROR: Directive SET may not be called within a compute region
  !$acc set default_async(i)
  !$acc end parallel

  !$acc serial
  !ERROR: Directive SET may not be called within a compute region
  !$acc set default_async(i)
  !$acc end serial

  !$acc kernels
  !ERROR: Directive SET may not be called within a compute region
  !$acc set default_async(i)
  !$acc end kernels

  !$acc parallel
  !$acc loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc serial
  !$acc loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do
  !$acc end serial

  !$acc kernels
  !$acc loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do
  !$acc end kernels

  !$acc parallel loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do

  !$acc serial loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do

  !$acc kernels loop
  do i = 1, N
    !ERROR: Directive SET may not be called within a compute region
    !$acc set default_async(i)
    a(i) = 3.14
  end do

  !ERROR: At least one of DEFAULT_ASYNC, DEVICE_NUM, DEVICE_TYPE clause must appear on the SET directive
  !$acc set

  !ERROR: At least one of DEFAULT_ASYNC, DEVICE_NUM, DEVICE_TYPE clause must appear on the SET directive
  !$acc set if(.TRUE.)

  !ERROR: At most one DEFAULT_ASYNC clause can appear on the SET directive
  !$acc set default_async(2) default_async(1)

  !ERROR: At most one DEFAULT_ASYNC clause can appear on the SET directive
  !$acc set default_async(2) default_async(1)

  !ERROR: At most one DEVICE_NUM clause can appear on the SET directive
  !$acc set device_num(1) device_num(i)

  !ERROR: At most one DEVICE_TYPE clause can appear on the SET directive
  !$acc set device_type(i) device_type(2, i, j)

  !$acc set default_async(2)
  !$acc set default_async(i)
  !$acc set device_num(1)
  !$acc set device_num(i)
  !$acc set device_type(i)
  !$acc set device_type(2, i, j)
  !$acc set device_num(1) default_async(2) device_type(2, i, j)

  !ERROR: At least one of ATTACH, COPYIN, CREATE clause must appear on the ENTER DATA directive
  !$acc enter data

  !ERROR: Modifier is not allowed for the COPYIN clause on the ENTER DATA directive
  !$acc enter data copyin(zero: i)

  !ERROR: Modifier is not allowed for the COPYOUT clause on the EXIT DATA directive
  !$acc exit data copyout(zero: i)

  !ERROR: Only the ZERO modifier is allowed for the CREATE clause on the ENTER DATA directive
  !$acc enter data create(readonly: i)

  !ERROR: Only the ZERO modifier is allowed for the COPYOUT clause on the DATA directive
  !$acc data copyout(readonly: i)
  !$acc end data

  !ERROR: COPYOUT clause is not allowed on the ENTER DATA directive
  !$acc enter data copyin(i) copyout(i)

  !ERROR: At most one IF clause can appear on the DATA directive
  !$acc data copy(i) if(.true.) if(.true.)
  !$acc end data

  !ERROR: At least one of COPYOUT, DELETE, DETACH clause must appear on the EXIT DATA directive
  !$acc exit data

  !ERROR: At least one of USE_DEVICE clause must appear on the HOST_DATA directive
  !$acc host_data
  !$acc end host_data

  !$acc host_data use_device(aa)
  !$acc end host_data

  !$acc host_data use_device(aa) if(.true.)
  !$acc end host_data

  !$acc host_data use_device(aa) if(ifCondition)
  !$acc end host_data

  !$acc host_data use_device(aa, bb) if_present
  !$acc end host_data

  !$acc host_data use_device(aa, bb) if(.true.) if_present
  !$acc end host_data

  !ERROR: At least one of DEFAULT_ASYNC, DEVICE_NUM, DEVICE_TYPE clause must appear on the SET directive
  !$acc set

  !ERROR: At least one of ATTACH, COPY, COPYIN, COPYOUT, CREATE, DEFAULT, DEVICEPTR, NO_CREATE, PRESENT clause must appear on the DATA directive
  !$acc data
  !$acc end data

  !$acc data copy(aa) if(.true.)
  !$acc end data

  !$acc data copy(aa) if(ifCondition)
  !$acc end data

  !$acc data copy(aa, bb, cc)
  !$acc end data

  !$acc data copyin(aa) copyin(readonly: bb) copyout(cc)
  !$acc end data

  !$acc data copyin(readonly: aa, bb) copyout(zero: cc)
  !$acc end data

  !$acc data create(aa, bb(:,:)) create(zero: cc(:,:))
  !$acc end data

  !$acc data no_create(aa) present(bb, cc)
  !$acc end data

  !$acc data deviceptr(aa) attach(bb, cc)
  !$acc end data

  !$acc data copy(aa, bb) default(none)
  !$acc end data

  !$acc data copy(aa, bb) default(present)
  !$acc end data

  !ERROR: At most one DEFAULT clause can appear on the DATA directive
  !$acc data copy(aa, bb) default(none) default(present)
  !$acc end data

  !ERROR: At most one IF clause can appear on the DATA directive
  !$acc data copy(aa) if(.true.) if(ifCondition)
  !$acc end data

  !$acc data copyin(i)
  !ERROR: Unmatched PARALLEL directive
  !$acc end parallel

  !ERROR: At least one of DEVICE, HOST, SELF clause must appear on the UPDATE directive
  !$acc update

  !$acc update self(a, f) host(g) device(h)

  !$acc update host(aa) async(1)

  !$acc update device(bb) async(async1)

  !ERROR: At most one ASYNC clause can appear on the UPDATE directive
  !$acc update host(aa, bb) async(1) async(2)

  !$acc update self(bb, cc(:)) wait(1)

  !ERROR: SELF clause on the UPDATE directive must have a var-list
  !$acc update self

  !$acc update device(aa, bb, cc) wait(wait1)

  !$acc update host(aa) host(bb) device(cc) wait(1,2)

  !$acc update device(aa, cc) wait(wait1, wait2)

  !$acc update device(aa) device_type(*) async

  !$acc update host(bb) device_type(*) wait

  !$acc update self(cc) device_type(1,2) async device_type(3) wait

  !ERROR: At most one IF clause can appear on the UPDATE directive
  !$acc update device(aa) if(.true.) if(ifCondition)

  !ERROR: At most one IF_PRESENT clause can appear on the UPDATE directive
  !$acc update device(bb) if_present if_present

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the UPDATE directive
  !$acc update device(i) device_type(*) if(.TRUE.)

  !$acc parallel
  !ERROR: INDEPENDENT and SEQ clauses are mutually exclusive and may not appear on the same LOOP directive
  !$acc loop seq independent
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: SEQ and AUTO clauses are mutually exclusive and may not appear on the same LOOP directive
  !$acc loop auto seq
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop tile(2)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel loop tile(2)
  do i = 1, N
    a(i) = 3.14
  end do

  !$acc parallel loop self
  do i = 1, N
    a(i) = 3.14
  end do

  !ERROR: SELF clause on the PARALLEL LOOP directive only accepts optional scalar logical expression
  !$acc parallel loop self(bb, cc(:))
  do i = 1, N
    a(i) = 3.14
  end do

  !$acc parallel loop self(.true.)
  do i = 1, N
    a(i) = 3.14
  end do

  !$acc parallel loop self(ifCondition)
  do i = 1, N
    a(i) = 3.14
  end do

  !$acc parallel loop tile(2, 2)
  do i = 1, N
    do j = 1, N
      aa(i, j) = 3.14
    end do
  end do

  !$acc parallel device_type(*) num_gangs(2)
  !$acc loop
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop seq
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop independent
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop auto
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: At most one VECTOR clause can appear on the LOOP directive
  !$acc loop vector vector(128)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop vector
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop vector(10)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop vector(vector_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop vector(length: vector_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: At most one WORKER clause can appear on the LOOP directive
  !$acc loop worker worker(10)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop worker
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop worker(10)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop worker(worker_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop worker(num: worker_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: At most one GANG clause can appear on the LOOP directive
  !$acc loop gang gang(gang_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop gang(gang_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop gang(num: gang_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop gang(gang_size, static:*)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop gang(num: gang_size, static:*)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop gang(num: gang_size, static: gang_size)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop private(b, a(:))
  do i = 1, N
    a(i) = b
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop tile(*)
  do i = 1, N
    a(i) = b
  end do
  !$acc end parallel

  !$acc parallel
  !$acc loop tile(2, 2)
  do i = 1, N
    do j = 1, N
      a(i) = b
    end do
  end do
  !$acc end parallel

  !$acc parallel async
  !$acc end parallel

  !$acc parallel async(1)
  !$acc end parallel

  !$acc parallel async(async1)
  !$acc end parallel

  !$acc parallel wait
  !$acc end parallel

  !$acc parallel wait(1)
  !$acc end parallel

  !$acc parallel wait(wait1)
  !$acc end parallel

  !$acc parallel wait(1,2)
  !$acc end parallel

  !$acc parallel wait(wait1, wait2)
  !$acc end parallel

  !$acc parallel num_gangs(8)
  !$acc end parallel

  !$acc parallel num_workers(8)
  !$acc end parallel

  !$acc parallel vector_length(128)
  !$acc end parallel

  !$acc parallel if(.true.)
  !$acc end parallel

  !$acc parallel if(ifCondition)
  !$acc end parallel

  !$acc parallel self
  !$acc end parallel

  !$acc parallel self(.true.)
  !$acc end parallel

  !$acc parallel self(ifCondition)
  !$acc end parallel

  !$acc parallel copy(aa) copyin(bb) copyout(cc)
  !$acc end parallel

  !$acc parallel copy(aa, bb) copyout(zero: cc)
  !$acc end parallel

  !$acc parallel present(aa, bb) create(cc)
  !$acc end parallel

  !$acc parallel copyin(readonly: aa, bb) create(zero: cc)
  !$acc end parallel

  !$acc parallel deviceptr(aa, bb) no_create(cc)
  !$acc end parallel

  !$acc parallel attach(aa, bb, cc)
  !$acc end parallel

  !$acc parallel private(aa) firstprivate(bb, cc)
  !$acc end parallel

  !$acc parallel default(none)
  !$acc end parallel

  !$acc parallel default(present)
  !$acc end parallel

  !$acc parallel device_type(*)
  !$acc end parallel

  !$acc parallel device_type(1)
  !$acc end parallel

  !$acc parallel device_type(1, 3)
  !$acc end parallel

  !ERROR: Clause PRIVATE is not allowed after clause DEVICE_TYPE on the PARALLEL directive
  !ERROR: Clause FIRSTPRIVATE is not allowed after clause DEVICE_TYPE on the PARALLEL directive
  !$acc parallel device_type(*) private(aa) firstprivate(bb)
  !$acc end parallel

  !$acc parallel device_type(*) async
  !$acc end parallel

  !$acc parallel device_type(*) wait
  !$acc end parallel

  !$acc parallel device_type(*) num_gangs(8)
  !$acc end parallel

  !$acc parallel device_type(1) async device_type(2) wait
  !$acc end parallel

  !$acc parallel
  !ERROR: The parameter of the COLLAPSE clause must be a constant positive integer expression
  !$acc loop collapse(-1)
  do i = 1, N
    do j = 1, N
      a(i) = 3.14 + j
    end do
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: Clause PRIVATE is not allowed after clause DEVICE_TYPE on the LOOP directive
  !$acc loop device_type(*) private(i)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: Clause GANG is not allowed if clause SEQ appears on the LOOP directive
  !$acc loop gang seq
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: Clause WORKER is not allowed if clause SEQ appears on the LOOP directive
  !$acc loop worker seq
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !$acc parallel
  !ERROR: Clause VECTOR is not allowed if clause SEQ appears on the LOOP directive
  !$acc loop vector seq
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the PARALLEL directive
  !$acc parallel device_type(*) if(.TRUE.)
  !$acc loop
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the PARALLEL LOOP directive
  !$acc parallel loop device_type(*) if(.TRUE.)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end parallel loop

  !$acc serial
  !$acc end serial

  !$acc serial async
  !$acc end serial

  !$acc serial async(1)
  !$acc end serial

  !ERROR: At most one ASYNC clause can appear on the SERIAL directive
  !$acc serial async(1) async(2)
  !$acc end serial

  !$acc serial async(async1)
  !$acc end serial

  !$acc serial wait
  !$acc end serial

  !$acc serial wait(1)
  !$acc end serial

  !$acc serial wait(wait1)
  !$acc end serial

  !$acc serial wait(1,2)
  !$acc end serial

  !$acc serial wait(wait1, wait2)
  !$acc end serial

  !$acc serial wait(wait1) wait(wait2)
  !$acc end serial

  !ERROR: NUM_GANGS clause is not allowed on the SERIAL directive
  !$acc serial num_gangs(8)
  !$acc end serial

  !ERROR: NUM_WORKERS clause is not allowed on the SERIAL directive
  !$acc serial num_workers(8)
  !$acc end serial

  !ERROR: VECTOR_LENGTH clause is not allowed on the SERIAL directive
  !$acc serial vector_length(128)
  !$acc end serial

  !$acc serial if(.true.)
  !$acc end serial

  !ERROR: At most one IF clause can appear on the SERIAL directive
  !$acc serial if(.true.) if(ifCondition)
  !$acc end serial

  !$acc serial if(ifCondition)
  !$acc end serial

  !$acc serial self
  !$acc end serial

  !$acc serial self(.true.)
  !$acc end serial

  !$acc serial self(ifCondition)
  !$acc end serial

  !$acc serial loop reduction(+: reduction_r)
  do i = 1, N
    reduction_r = a(i) + i
  end do

  !$acc serial loop reduction(*: reduction_r)
  do i = 1, N
    reduction_r = reduction_r * (a(i) + i)
  end do

  !$acc serial loop reduction(min: reduction_r)
  do i = 1, N
    reduction_r = min(reduction_r, a(i) * i)
  end do

  !$acc serial loop reduction(max: reduction_r)
  do i = 1, N
    reduction_r = max(reduction_r, a(i) * i)
  end do

  !$acc serial loop reduction(iand: b)
  do i = 1, N
    b = iand(b, c(i))
  end do

  !$acc serial loop reduction(ior: b)
  do i = 1, N
    b = ior(b, c(i))
  end do

  !$acc serial loop reduction(ieor: b)
  do i = 1, N
    b = ieor(b, c(i))
  end do

  !$acc serial loop reduction(.and.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .and. e(i)
  end do

  !$acc serial loop reduction(.or.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .or. e(i)
  end do

  !$acc serial loop reduction(.eqv.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .eqv. e(i)
  end do

  !$acc serial loop reduction(.neqv.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .neqv. e(i)
  end do

  !$acc serial reduction(.neqv.: reduction_l)
  !$acc loop reduction(.neqv.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .neqv. e(i)
  end do
  !$acc end serial

  !$acc serial copy(aa) copyin(bb) copyout(cc)
  !$acc end serial

  !$acc serial copy(aa, bb) copyout(zero: cc)
  !$acc end serial

  !$acc serial present(aa, bb) create(cc)
  !$acc end serial

  !$acc serial copyin(readonly: aa, bb) create(zero: cc)
  !$acc end serial

  !$acc serial deviceptr(aa, bb) no_create(cc)
  !$acc end serial

  !$acc serial attach(aa, bb, cc)
  !$acc end serial

  !$acc serial firstprivate(bb, cc)
  !$acc end serial

  !$acc serial private(aa)
  !$acc end serial

  !$acc serial default(none)
  !$acc end serial

  !$acc serial default(present)
  !$acc end serial

  !ERROR: At most one DEFAULT clause can appear on the SERIAL directive
  !$acc serial default(present) default(none)
  !$acc end serial

  !$acc serial device_type(*) async wait
  !$acc end serial

  !$acc serial device_type(*) async
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end serial

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the SERIAL directive
  !$acc serial device_type(*) if(.TRUE.)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end serial

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the SERIAL LOOP directive
  !$acc serial loop device_type(*) if(.TRUE.)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end serial loop

  !$acc parallel loop
  do i = 1, N
    a(i) = 3.14
  end do
  !ERROR: Unmatched END KERNELS LOOP directive
  !$acc end kernels loop

  !$acc kernels loop
  do i = 1, N
    a(i) = 3.14
  end do
  !ERROR: Unmatched END SERIAL LOOP directive
  !$acc end serial loop

  !$acc serial loop
  do i = 1, N
    a(i) = 3.14
  end do
  !ERROR: Unmatched END PARALLEL LOOP directive
  !$acc end parallel loop

  !$acc parallel loop reduction(+: reduction_r)
  do i = 1, N
    reduction_r = a(i) + i
  end do

  !$acc parallel loop reduction(*: reduction_r)
  do i = 1, N
    reduction_r = reduction_r * (a(i) + i)
  end do

  !$acc parallel loop reduction(min: reduction_r)
  do i = 1, N
    reduction_r = min(reduction_r, a(i) * i)
  end do

  !$acc parallel loop reduction(max: reduction_r)
  do i = 1, N
    reduction_r = max(reduction_r, a(i) * i)
  end do

  !$acc parallel loop reduction(iand: b)
  do i = 1, N
    b = iand(b, c(i))
  end do

  !$acc parallel loop reduction(ior: b)
  do i = 1, N
    b = ior(b, c(i))
  end do

  !$acc parallel loop reduction(ieor: b)
  do i = 1, N
    b = ieor(b, c(i))
  end do

  !$acc parallel loop reduction(.and.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .and. e(i)
  end do

  !$acc parallel loop reduction(.or.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .or. e(i)
  end do

  !$acc parallel loop reduction(.eqv.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .eqv. e(i)
  end do

  !$acc parallel loop reduction(.neqv.: reduction_l)
  do i = 1, N
    reduction_l = d(i) .neqv. e(i)
  end do

  !$acc kernels async
  !$acc end kernels

  !$acc kernels async(1)
  !$acc end kernels

  !$acc kernels async(async1)
  !$acc end kernels

  !$acc kernels wait(wait1)
  !$acc end kernels

  !$acc kernels wait(wait1, wait2)
  !$acc end kernels

  !$acc kernels wait(1, 2) async(3)
  !$acc end kernels

  !$acc kernels wait(queues: 1, 2) async(3)
  !$acc end kernels

  !$acc kernels wait(1) wait(2) async(3)
  !$acc end kernels

  !$acc kernels wait(devnum: 1: 1, 2) async(3)
  !$acc end kernels

  !$acc kernels wait(devnum: 1: queues: 1, 2) async(3)
  !$acc end kernels

  !$acc kernels num_gangs(8)
  !$acc end kernels

  !$acc kernels num_workers(8)
  !$acc end kernels

  !$acc kernels vector_length(128)
  !$acc end kernels

  !$acc kernels if(.true.)
  !$acc end kernels

  !$acc kernels if(ifCondition)
  !$acc end kernels

  !ERROR: At most one IF clause can appear on the KERNELS directive
  !$acc kernels if(.true.) if(ifCondition)
  !$acc end kernels

  !$acc kernels self
  !$acc end kernels

  !$acc kernels self(.true.)
  !$acc end kernels

  !$acc kernels self(ifCondition)
  !$acc end kernels

  !$acc kernels copy(aa) copyin(bb) copyout(cc)
  !$acc end kernels

  !$acc kernels copy(aa, bb) copyout(zero: cc)
  !$acc end kernels

  !$acc kernels present(aa, bb) create(cc)
  !$acc end kernels

  !$acc kernels copyin(readonly: aa, bb) create(zero: cc)
  !$acc end kernels

  !$acc kernels deviceptr(aa, bb) no_create(cc)
  !$acc end kernels

  !$acc kernels attach(aa, bb, cc)
  !$acc end kernels

  !ERROR: PRIVATE clause is not allowed on the KERNELS directive
  !$acc kernels private(aa, bb, cc)
  !$acc end kernels

  !$acc kernels default(none)
  !$acc end kernels

  !$acc kernels default(present)
  !$acc end kernels

  !ERROR: At most one DEFAULT clause can appear on the KERNELS directive
  !$acc kernels default(none) default(present)
  !$acc end kernels

  !$acc kernels device_type(*)
  !$acc end kernels

  !$acc kernels device_type(1)
  !$acc end kernels

  !$acc kernels device_type(1, 3)
  !$acc end kernels

  !$acc kernels device_type(*) async wait num_gangs(8) num_workers(8) vector_length(128)
  !$acc end kernels

  !$acc kernels device_type(*) async
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end kernels

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the KERNELS directive
  !$acc kernels device_type(*) if(.TRUE.)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end kernels

  !ERROR: Clause IF is not allowed after clause DEVICE_TYPE on the KERNELS LOOP directive
  !$acc kernels loop device_type(*) if(.TRUE.)
  do i = 1, N
    a(i) = 3.14
  end do
  !$acc end kernels loop

  !$acc wait

  !$acc wait async

  !$acc wait(1)
  !$acc wait(1, 2)

  !$acc wait(queues: 1)
  !$acc wait(queues: 1, 2)

  !$acc wait(devnum: 1: 3)
  !$acc wait(devnum: 1: 3, 4)

  !$acc wait(devnum: 1: queues: 3)
  !$acc wait(devnum: 1: queues: 3, 4)

  !$acc wait(1) if(.true.)

  !ERROR: At most one IF clause can appear on the WAIT directive
  !$acc wait(1) if(.true.) if(.false.)

  !$acc wait(1) if(.true.) async

  !$acc wait(1) if(.true.) async(1)

  !ERROR: At most one ASYNC clause can appear on the WAIT directive
  !$acc wait(1) if(.true.) async(1) async

  !$acc parallel
  !$acc atomic update
  c(i) = c(i) + 1

  !$acc atomic update
  c(i) = c(i) + 1
  !$acc end atomic

  !$acc atomic write
  c(i) = 10

  !$acc atomic write
  c(i) = 10
  !$acc end atomic

  !$acc atomic read
  i = c(i)

  !$acc atomic read
  i = c(i)
  !$acc end atomic

  !$acc atomic capture
  c(i) = i
  i = i + 1
  !$acc end atomic
  !$acc end parallel
  t%arr(i) = 2.0

  !$acc cache(a(i))
  !$acc cache(a(1:2,3:4))
  !$acc cache(a)
  !$acc cache(readonly: a, aa)
  !$acc cache(readonly: a(i), aa(i, i))
  !$acc cache(t%arr)
  !$acc cache(ta(1:2)%arr)
  !$acc cache(ta(1:2)%arr(1:4))

  !ERROR: Only array element or subarray are allowed in CACHE directive
  !$acc cache(ta(1:2)%s)

  !ERROR: Only array element or subarray are allowed in CACHE directive
  !$acc cache(i)

  !ERROR: Only array element or subarray are allowed in CACHE directive
  !$acc cache(t%s)

  !ERROR: Only array element or subarray are allowed in CACHE directive
  !$acc cache(/i/)

end program openacc_clause_validity

max_stars_repo_name: glyons-intel/llvm

max_issues_repo_name: glyons-intel/llvm

max_forks_repo_name: glyons-intel/llvm


text:
      subroutine newt(nit0,nitmx,nit1,nit2,nit3,conv,omega,
     &  ichk0,donon,cutoff,i_phimap,i_epsmap,i_debmap,i_qmap,icount1b,
     &  i_iepsmp,i_idebmap,i_gchrg,i_gchrgp,i_medeps,natom,nobject,
     &  nmedia)

c------------------------------------------------------------------
	include 'qdiffpar5.h'
	include 'qlog.h'
c------------------------------------------------------------------
	pointer (i_gchrg,gchrg),(i_gchrgp,gchrgp)
      pointer (i_phimap,phimap),(i_medeps,medeps)
	pointer (i_iepsmp,iepsmp),(i_idebmap,idebmap)
	pointer (i_epsmap,epsmap)
	pointer (i_debmap,debmap),(i_qmap,qmap)

c       W da pointer ma passati
        real phimap(igrid,igrid,igrid)
        real medeps(0:nmedia)
        integer iepsmp(igrid,igrid,igrid,3)
        logical*1 idebmap(igrid,igrid,igrid)
        dimension gchrg(icount1b)
        integer gchrgp(3,icount1b)

c       W passati con pointer aggiunto da me
        real*4 epsmap(igrid,igrid,igrid,3)
        real*4 debmap(igrid,igrid,igrid)
        real*4 qmap(igrid,igrid,igrid)

        logical donon

c       W locali
        real*4 deemap(igrid,igrid,igrid)
        real*4 delmap(igrid,igrid,igrid)
        real*4 eff(igrid,igrid,igrid)
        real*8 normf,normfp,normf1,lam,lamtol
c==================================================================
c inexact Newton method for solving the nonlinear Poisson-Boltzmann
c eqn.
c based on Mike Holst's algorithm, implemented by Jonathan Hecht
c refined by kim sharp, oct 1993
c
c==================================================================

c
c the NLPB eqn is:
c
c  div(eps.grad phi) - K**2.sinh(phi) + 4pirho(fixed) = 0
c
c the finite difference form of the PB eqn is
c
c  sum(epsj.phij) - phi0.sum(epsj) - (K.h)**2.sinh(phi0) + 4pi.q0/h
c    = 0
c
c where 0 refers the the values at tha grid point, the sum is over
c the j
c index for the 6 neighbors, h is the spacing between grid points
c phi is stored in array phimap, eps is stored in array epsmap,
c (K.h)**2 is stored in array debmap, and 4piq/h is stored in qmap
c
c The basic outline of the algorithm is:
c
c From current phi estimate, getf calculates the residual, fmap,
c the amount
c that the lhs of the finite difference PB eqn. differes from 0 at
c each
c lattice point.  The size of normf (the sum of the squares of the
c residual)
c is the criterion for convergence use din the outer, Newton method
c loop
c
c The residuals are then used as the source term in the differential
c form
c of the PB equation used to estimate the correction to the
c potential, dphi:
c
c  sum(epsj.dphij) - dphi0.sum(epsj) - (K.h)**2.dphi0.cosh(phi0) +
c    fmap0 = 0
c
c where sinh(phi) has been replaced by its differential,
c dph.cosh(phi),
c and the charge term is replaced by the residual.  Delmap contains
c the correction to the potential, dph.  The residual of THIS
c equation is calculated in multin, stored in resid, and the
c estimated dphi is
c refined using the standard gauss-seidel relaxtaion scheme, where
c
c dphi0(new) = dphi0 + omega*resid/(sum(epsj) + (K.h)**2.cosh(phi0))
c
c until the normal of this residual, rnorm, reaches the required
c value.
c The relaxation parameter omega, and multigridding if necessary
c aids convergence. omega=1.6 appears to be a good value.
c The potential is updated in the main newton loop
c using phi = phi+lambda*dphi.  Initially we try lambda = 1, and
c recalculate
c the normal, normf. If this is not decreased, then we try an ever
c smaller
c lambda (decreasing it by a factor of a two each time) until the
c normal
c is decreased.  The nonlinear sinh and cosh terms are updates, the
c differential equation is set up, and dphi solved for again, etc...
c until convergence in the outer loop is obtained.  Note that once we
c find the normal decreasing with a lambda of 1, it will always
c converge with lambda=1 and so we omit the lambda search thereafter
c
c====================================================================

c set up arrays
c check for pts assigned charge that are in salt: blows nonlinaer
c equation up!
c

c Added by Joe ------------------------------------
        pi=acos(-1.0)
        debfct = rionst * 8.0 * pi * 6.022136736e-4 / (scale**2)

        do k=1,igrid
          do j=1,igrid
            do i=1,igrid
              if(idebmap(i,j,k)) then
                debmap(i,j,k)=debfct
              else
                debmap(i,j,k)=0.0
              end if

              do m=1,3
c                if(iepsmp(i,j,k,m).gt.1) then
c                  itmp=iepsmp(i,j,k,m)/(natom+nobject+2)
c                  epsmap(i,j,k,m)=medeps(itmp)
c                else
c                  epsmap(i,j,k,m)=epsout
c                end if
                 itmp=iepsmp(i,j,k,m)/(natom+nobject+2)
                 epsmap(i,j,k,m)=medeps(itmp)
              end do
            end do
          end do
        end do

        do i=1,icount1b
          qmap(gchrgp(1,i),gchrgp(2,i),gchrgp(3,i))
     &      =4*pi*scale*gchrg(i)
        end do

c -------------------------------------------------

        start = cputime(0.0)
        print *,'using inexact Newton method to solve NLPB eqn...'
        print *,'dphi relaxation parameter: ',omega
        print *,'dphi multigrid iterations: ',nit0,nit1,nit2,nit3
        print *,' '
        nd1q1 = 0
        do k = 2,igrid-1
          do j = 2,igrid-1
            do i = 2,igrid-1
            deemap(i,j,k) = epsmap(i,j,k,1) + epsmap(i,j,k,2)
     &      + epsmap(i,j,k,3) + epsmap(i-1,j,k,1)
     &      + epsmap(i,j-1,k,2) + epsmap(i,j,k-1,3)
            delmap(i,j,k) = 0.0
                phimap(i,j,k) = 0.
            end do
          end do
        end do

* We need to reset the dielectric
* constants at grid points where the partial charges are not zero

        do k = 1,igrid
          do j = 1,igrid
            do i = 1,igrid
                if((abs(qmap(i,j,k)).gt.0.).and.
     &            (debmap(i,j,k).gt.0.))then
                  nd1q1 = nd1q1+1
                  debmap(i,j,k) = 0.
                  idebmap(i,j,k) = .false.
                end if
            end do
          end do
        end do
        if(nd1q1.gt.0)then
          print *,'WARNING: # of charged points in salt: ',nd1q1
          print *,
     &      'these were set out of salt; scale should be increased'
        end if
c
c intialize 'source' term of PB eqn., which is residual
c and calculate its normal
c
        ntot = 0
        lamtol = 1./2.**15
        print *,'before getf'
        call getf(phimap,epsmap,qmap,debmap,deemap,
     &    eff,normf,donon,cutoff)
        print *,'After getf'
c
c begin outer, newton iteration loop until convergence criteria is
c met
c or maxmum iterations
c
        do while((normf.gt.conv).and.(ntot.lt.nitmx))
          ntot = ntot + 1
          print *,'NEWTON ITERATION#,norm: ',ntot,normf
c
c get new correction, delmap, to potential
c
          ngrid1 = (igrid+1)/2
          ngrid2 = (ngrid1+1)/2
          ngrid3 = (ngrid2+1)/2
          call multin(i_phimap,i_debmap,i_epsmap,deemap,
     &      nit0,nit1,nit2,nit3,omega,delmap,eff,
     &      normf,ichk0,donon,cutoff,nitmx,ngrid1,ngrid2,ngrid3)

c
c apply total correction to potential
c
          do iz=2,igrid-1
            do iy=2,igrid-1
              do ix=2,igrid-1
                phimap(ix,iy,iz)=phimap(ix,iy,iz)+delmap(ix,iy,iz)
              end do
            end do
        end do
c
c periodic boundary conditions
c
          call doper(phimap,igrid,iper)
c
c find new normal and check to see if converging
c
          call getf(phimap,epsmap,qmap,debmap,deemap,
     &      eff,normfp,donon,cutoff)
          if(normfp.ge.normf)then
c
c not converging, so subtract of half of correction until it does
c
            lam = 0.5
            do while((normfp.ge.normf).and.(lam.ge.lamtol))
              do iz=2,igrid-1
                do iy=2,igrid-1
                  do ix=2,igrid-1
                    phimap(ix,iy,iz)=phimap(ix,iy,iz)-
     &                lam*delmap(ix,iy,iz)
                  end do
                end do
            end do
              call doper(phimap,igrid,iper)
c
c get the new residual, its normal and check lambda
c
              call getf(phimap,epsmap,qmap,debmap,deemap,
     &          eff,normfp,donon,cutoff)
              print *,'trial lambda, normal: ',lam,normfp
              lam = lam/2.
            end do
            if(lam.lt.lamtol)print *,'WARNING: did not converge
     &        on lambda'
          end if
          normf = normfp
        end do
        fin = cputime(start)
        print *,'final normal: ',normf
        if(ntot.ge.nitmx)print *,'WARNING: max # of newton
     &    iterations: ',ntot
        print *,'total time for newton iterations: ',fin

c        nphi=0
c
c        do k=1,igrid
c          do j=1,igrid
c            do i=1,igrid
c              if(phimap(i,j,k).eq.0.) nphi=nphi+1
c            end do
c          end do
c        end do
c        print *,'total number of times when phiimap is 0 ',nphi

        return
        end

        subroutine getf(phimap,epsmap,qmap,debmap,deemap,
     &    fmap,normf,donon,cutoff)
c
c this gets the new source term for the current newton iteration
c
c-------------------------------------------------
	include 'qdiffpar5.h'
	include 'qlog.h'
c-------------------------------------------------

      real*4 phimap(igrid,igrid,igrid)
      real*4 epsmap(igrid,igrid,igrid,3)
      real*4 qmap(igrid,igrid,igrid)
      real*4 debmap(igrid,igrid,igrid)
      real*4 deemap(igrid,igrid,igrid)

        real*4 fmap(igrid,igrid,igrid),pold,pold2
        real*8 normf
        logical donon
c-------------------------------------------------
        normf = 0.
        if(donon)then
        do iz = 2, igrid-1
          do iy = 2, igrid-1
            do ix = 2, igrid-1
                  pold = phimap(ix,iy,iz)
                  fmap(ix,iy,iz) =
     &          phimap(ix,iy,iz)*deemap(ix,iy,iz) - qmap(ix,iy,iz)
     &          - phimap(ix-1,iy,iz)*epsmap(ix-1,iy,iz,1)
     &        - phimap(ix+1,iy,iz)*epsmap(ix,iy,iz,1)
     &        - phimap(ix,iy-1,iz)*epsmap(ix,iy-1,iz,2)
     &        - phimap(ix,iy+1,iz)*epsmap(ix,iy,iz,2)
     &        - phimap(ix,iy,iz-1)*epsmap(ix,iy,iz-1,3)
     &        - phimap(ix,iy,iz+1)*epsmap(ix,iy,iz,3)
                if(debmap(ix,iy,iz).ne.0.) then
                pold= max(phimap(ix,iy,iz),-cutoff)
                  pold= min(pold,cutoff)
                    fmap(ix,iy,iz) = fmap(ix,iy,iz) +
     &                debmap(ix,iy,iz)*sinh(pold)
                end if
                  normf = normf + fmap(ix,iy,iz)**2
              end do
            end do
          end do
        else
        do iz = 2, igrid-1
          do iy = 2, igrid-1
            do ix = 2, igrid-1
                  pold = phimap(ix,iy,iz)
                  fmap(ix,iy,iz) =
     &          phimap(ix,iy,iz)*(deemap(ix,iy,iz)+debmap(ix,iy,iz))
     &        - qmap(ix,iy,iz)
     &          - phimap(ix-1,iy,iz)*epsmap(ix-1,iy,iz,1)
     &        - phimap(ix+1,iy,iz)*epsmap(ix,iy,iz,1)
     &        - phimap(ix,iy-1,iz)*epsmap(ix,iy-1,iz,2)
     &        - phimap(ix,iy+1,iz)*epsmap(ix,iy,iz,2)
     &        - phimap(ix,iy,iz-1)*epsmap(ix,iy,iz-1,3)
     &        - phimap(ix,iy,iz+1)*epsmap(ix,iy,iz,3)
                  normf = normf + fmap(ix,iy,iz)**2
              end do
            end do
          end do
        end if
        return
        end

        subroutine getcoh(phimap,debmap,coh,donon,cutoff)
c
c calculates the hyperbolic cosines for this iteration
c-------------------------------------------------
	include 'qdiffpar5.h'
	include 'qlog.h'
c-------------------------------------------------

        real*4 phimap(igrid,igrid,igrid)
        real*4 debmap(igrid,igrid,igrid)

        real*4 coh(igrid,igrid,igrid)
        logical donon
c-------------------------------------------------------
        if(donon)then
        do iz = 2, igrid-1
          do iy = 2, igrid-1
            do ix = 2, igrid-1
                if (debmap(ix,iy,iz).eq.0.) then
                  coh(ix,iy,iz) = 0.
                else
                temp= max(phimap(ix,iy,iz),-cutoff)
                  temp= min(temp,cutoff)
                  coh(ix,iy,iz) = cosh(temp)*debmap(ix,iy,iz)
              end if
              end do
            end do
          end do
        else
        do iz = 2, igrid-1
          do iy = 2, igrid-1
            do ix = 2, igrid-1
                  coh(ix,iy,iz) = debmap(ix,iy,iz)
              end do
            end do
          end do
        end if
        return
        end

        subroutine multin(i_phimap,i_debmap,i_epsmap,deemap,
     &    nit0,nit1s,nit2,nit3,omega,delmap,eff,normf,ichk0,
     &    donon,cutoff,nitmx,ngrid1,ngrid2,ngrid3)
c
c-------------------------------------------------
	include 'qdiffpar5.h'
	include 'qlog.h'
c-------------------------------------------------
      pointer (i_phimap,phimap),(i_epsmap,epsmap)
	pointer (i_debmap,debmap)
c        parameter (nitmx=300)
c
c uppper level (0) arrays
c

c	parameter (ngrid1 = (ngrid+1)/2)
c	parameter (ngrid2 = (ngrid1+1)/2)
c	parameter (ngrid3 = (ngrid2+1)/2)

      real*4 deemap(igrid,igrid,igrid)

      real*4 phimap(igrid,igrid,igrid)
      real*4 debmap(igrid,igrid,igrid)
      real*4 epsmap(igrid,igrid,igrid,3)

        real*4 delmap(igrid,igrid,igrid)
        real*4 eff(igrid,igrid,igrid)
        real*4 rsumeps(igrid,igrid,igrid)
        real*4 sumeps(igrid,igrid,igrid)
        real*4 resid(igrid,igrid,igrid)
        real*4 ddelmap(igrid,igrid,igrid)
        real*4 dresid(igrid,igrid,igrid)
        real*4 coh(igrid,igrid,igrid)

c level 1 arrays

        real*4 epsmap1(ngrid1,ngrid1,ngrid1,3)
        real*4 eff1(ngrid1,ngrid1,ngrid1)
        real*4 rsumeps1(ngrid1,ngrid1,ngrid1)
        real*4 sumeps1(ngrid1,ngrid1,ngrid1)
        real*4 resid1(ngrid1,ngrid1,ngrid1)
        real*4 ddelmap1(ngrid1,ngrid1,ngrid1)
        real*4 debmap1(ngrid1,ngrid1,ngrid1)
        real*4 delmap1(ngrid1,ngrid1,ngrid1)
        real*4 coh1(ngrid1,ngrid1,ngrid1)

c level 2 arrays

        real*4 epsmap2(ngrid2,ngrid2,ngrid2,3)
        real*4 eff2(ngrid2,ngrid2,ngrid2)
        real*4 rsumeps2(ngrid2,ngrid2,ngrid2)
        real*4 sumeps2(ngrid2,ngrid2,ngrid2)
        real*4 resid2(ngrid2,ngrid2,ngrid2)
        real*4 ddelmap2(ngrid2,ngrid2,ngrid2)
        real*4 debmap2(ngrid2,ngrid2,ngrid2)
        real*4 delmap2(ngrid2,ngrid2,ngrid2)
        real*4 coh2(ngrid2,ngrid2,ngrid2)

c level 3 arrays

        real*4 epsmap3(ngrid3,ngrid3,ngrid3,3)
        real*4 eff3(ngrid3,ngrid3,ngrid3)
        real*4 rsumeps3(ngrid3,ngrid3,ngrid3)
        real*4 sumeps3(ngrid3,ngrid3,ngrid3)
        real*4 resid3(ngrid3,ngrid3,ngrid3)
        real*4 debmap3(ngrid3,ngrid3,ngrid3)
        real*4 delmap3(ngrid3,ngrid3,ngrid3)
        real*4 coh3(ngrid3,ngrid3,ngrid3)

        character*8 hour
        real*8 normf,rnorm,rnorm1,rnorm2,rnorm3
        logical donon
        data zero / 0. /

c-------------------------------------------------------
        start = cputime(0.0)
c-------------------------------------------------------
c
c start
c
        nit1 = abs(nit1s)
        igrid1 = 0
        igrid2 = 0
        igrid3 = 0
c        print *,' '
        if(nit1.gt.0)then
          if(mod(igrid-1,16).eq.0)then
            igrid1 = (igrid+1)/2
            igrid2 = (igrid1+1)/2
            igrid3 = (igrid2+1)/2
c            print *,'using 4 level multigrid'
          else if(mod(igrid-1,8).eq.0)then
            nit3 = 0
            igrid1 = (igrid+1)/2
            igrid2 = (igrid1+1)/2
            igrid3 = 0
c            print *,'using 3 level multigrid'
          else if(mod(igrid-1,4).eq.0)then
            nit2 = 0
            nit3 = 0
            igrid1 = (igrid+1)/2
            igrid2 = 0
            igrid3 = 0
c            print *,'using 2 level multigrid'
          else
            nit1 = 0
            nit2 = 0
            nit3 = 0
            igrid1 = 0
            igrid2 = 0
            igrid3 = 0
c            print *,'multigridding switched off'
          end if
        else
c          print *,'multigridding switched off'
        end if
c        print *,'grid dimensions: ',ngrid,ngrid1,ngrid2,ngrid3
c        print *,'grid sizes     : ',igrid,igrid1,igrid2,igrid3
c        print *,'level 0,1,2,3 iterations: ',nit0,nit1,nit2,nit3
c===============================
c
c set up coefficient arrays
c
c===============================
c level 0
c===============================
c         get cosh * debmap array
        call getcoh(phimap,debmap,coh,donon,cutoff)
        call setarr(ddelmap,igrid,igrid,zero)
        call setarr(delmap,igrid,igrid,zero)
        do k = 2,igrid-1
          do j = 2,igrid-1
            do i = 2,igrid-1
                sumeps(i,j,k) = deemap(i,j,k)+coh(i,j,k)
                rsumeps(i,j,k) = 1./sumeps(i,j,k)
            end do
          end do
        end do
c===============================
c level 1
c===============================
c
c coarse epsmap by harmonic averaging
c
        call hrmav(epsmap,igrid,igrid,epsmap1,ngrid1,igrid1)
c
c coarse debmap directly
c
        ibound = 1
        call transfr(debmap,igrid,igrid,debmap1,ngrid1,igrid1,ibound)
        call transfr(coh,igrid,igrid,coh1,ngrid1,igrid1,ibound)
c
c secondary coefficients
c
        call setarr(ddelmap1,ngrid1,igrid1,zero)
        do k = 2,igrid1-1
          do j = 2,igrid1-1
            do i = 2,igrid1-1
                delmap1(i,j,k) = 0.
                debmap1(i,j,k) = debmap1(i,j,k)*4.
                coh1(i,j,k) = coh1(i,j,k)*4.
                sumeps1(i,j,k) = epsmap1(i,j,k,1) + epsmap1(i,j,k,2)
     &                       + epsmap1(i,j,k,3) + epsmap1(i-1,j,k,1)
     &                       + epsmap1(i,j-1,k,2) + epsmap1(i,j,k-1,3)
     &                       + coh1(i,j,k)
                rsumeps1(i,j,k) = 1./sumeps1(i,j,k)
            end do
          end do
        end do
c===============================
c level 2
c===============================
        call hrmav(epsmap1,ngrid1,igrid1,epsmap2,ngrid2,igrid2)
        ibound = 1
        call transfr(debmap1,ngrid1,igrid1,debmap2,ngrid2,
     &    igrid2,ibound)
        call transfr(coh1,ngrid1,igrid1,coh2,ngrid2,igrid2,ibound)
        do k = 2,igrid2-1
          do j = 2,igrid2-1
            do i = 2,igrid2-1
                delmap2(i,j,k) = 0.
                debmap2(i,j,k) = debmap2(i,j,k)*4.
                coh2(i,j,k) = coh2(i,j,k)*4.
                sumeps2(i,j,k) = epsmap2(i,j,k,1) + epsmap2(i,j,k,2)
     &                       + epsmap2(i,j,k,3) + epsmap2(i-1,j,k,1)
     &                       + epsmap2(i,j-1,k,2) + epsmap2(i,j,k-1,3)
     &                       + coh2(i,j,k)
                rsumeps2(i,j,k) = 1./sumeps2(i,j,k)
            end do
          end do
        end do
c===============================
c level 3
c===============================
        call hrmav(epsmap2,ngrid2,igrid2,epsmap3,ngrid3,igrid3)
        ibound = 1
        call transfr(debmap2,ngrid2,igrid2,debmap3,ngrid3,
     &    igrid3,ibound)
        call transfr(coh2,ngrid2,igrid2,coh3,ngrid3,igrid3,ibound)
        do k = 2,igrid3-1
          do j = 2,igrid3-1
            do i = 2,igrid3-1
                delmap3(i,j,k) = 0.
                debmap3(i,j,k) = debmap3(i,j,k)*4.
                coh3(i,j,k) = coh3(i,j,k)*4.
                sumeps3(i,j,k) = epsmap3(i,j,k,1) + epsmap3(i,j,k,2)
     &                       + epsmap3(i,j,k,3) + epsmap3(i-1,j,k,1)
     &                       + epsmap3(i,j-1,k,2) + epsmap3(i,j,k-1,3)
     &                       + coh3(i,j,k)
                rsumeps3(i,j,k) = 1./sumeps3(i,j,k)
            end do
          end do
        end do
c---------------------------------------------------------------
        om0 = omega
c
c start level 1,2 relaxation parameters at 1.8, max val is 1.9
c
C        om1 = 1.6
C        om2 = 1.6
C        om3 = 1.6
        om1 = 1.0
        om2 = 1.0
        om3 = 1.0
c        print *,'level 0,1,2,3 omegas: ',om0,om1,om2,om3
c
        ntot = 0
        icheck=0
         print *,'LEVEL  #it    Norm of dphi residual'
300        continue ! big loop over grids
c
c=============================================
c smoooth for nit0 iterations at LEVEL 0
c 1st time thru do at least n2 iterations
c=============================================
        do n = 1,nit0
          ntot=ntot+1
          if(ntot.gt.nitmx)then
            print *,'WARNING: exceeding maximum # of iterations: ',
     &        nitmx
            ntot = nitmx
            goto 301
          end if
c
c calculate top level residuals and update using chequerboard
c gauss-seidel
c do odd parity points first iodd = 0,
c do even parity points last iodd = 1, so that
c points on coarser grid (which are odd) are left with the residuals
c
          do iodd = 0,1

            call getres(delmap,epsmap,sumeps,eff,resid,igrid,
     &       igrid,iodd,rnorm)
            call gausei(delmap,resid,rsumeps,igrid,igrid,iodd,om0)
          end do
          write(6,200)' 0      ',ntot,rnorm
c          if (rnorm.le.normf) goto 301
200          format(a,i3,f13.6)
c
c periodic boundary conditions
c
        end do
c
c update residual
c
        do iodd = 0,1
            call getres(delmap,epsmap,sumeps,eff,resid,igrid,
     &        igrid,iodd,rnorm)
        end do
        write(6,200)' 0      ',ntot,rnorm
        icheck=icheck+1
        if(icheck.eq.ichk0)then
          icheck = 0
          if (rnorm.le.normf) goto 301
        end if
c
c do level 1 and below only if multigrid iterations
c
        if(nit1.eq.0)goto 300
c==================================================
c transfer residuals to LEVEL 1 and solve for error
c==================================================
        ibound = 0
        call transfr(resid,igrid,igrid,eff1,ngrid1,igrid1,ibound)
        call setarr(delmap1,ngrid1,igrid1,zero)
c==========================================
c presmooth at level 1
c==========================================
        do n = 1,nit1
          do iodd = 0,1
            call getres(delmap1,epsmap1,sumeps1,eff1,resid1,ngrid1,
     &    igrid1,iodd,rnorm1)
            call gausei(delmap1,resid1,rsumeps1,ngrid1,igrid1,
     &        iodd,om1)
          end do
        end do
        do iodd = 0,1
          call getres(delmap1,epsmap1,sumeps1,eff1,resid1,
     &  ngrid1,igrid1,iodd,rnorm1)
        end do
        write(6,200)'    pre1      ',n,rnorm1
        if(nit2.eq.0)goto 400
c==================================================
c transfer residuals to LEVEL 2 and solve for error in error
c==================================================
        ibound = 0
        call transfr(resid1,ngrid1,igrid1,eff2,ngrid2,igrid2,ibound)
        call setarr(delmap2,ngrid2,igrid2,zero)
c==========================================
c presmooth at level 2
c==========================================
        do n = 1,nit2
          do iodd = 0,1
            call getres(delmap2,epsmap2,sumeps2,eff2,resid2,
     &    ngrid2,igrid2,iodd,rnorm2)
            call gausei(delmap2,resid2,rsumeps2,ngrid2,igrid2,
     &        iodd,om2)
          end do
        end do
        do iodd = 0,1
          call getres(delmap2,epsmap2,sumeps2,eff2,resid2,
     &  ngrid2,igrid2,iodd,rnorm2)
        end do
        write(6,200)'      pre2      ',n,rnorm2
        if(nit3.eq.0)goto 500
c==================================================
c transfer residuals to LEVEL 3 and solve for error in error
c==================================================
        ibound = 0
        call transfr(resid2,ngrid2,igrid2,eff3,ngrid3,igrid3,ibound)
        call setarr(delmap3,ngrid3,igrid3,zero)
c==========================================
c smooth at level 3
c==========================================
        do n = 1,nit3
          do iodd = 0,1
            call getres(delmap3,epsmap3,sumeps3,eff3,resid3,
     &    ngrid3,igrid3,iodd,rnorm3)
            call gausei(delmap3,resid3,rsumeps3,ngrid3,igrid3,
     &        iodd,om3)
          end do
        end do
        do iodd = 0,1
          call getres(delmap3,epsmap3,sumeps3,eff3,
     &  resid3,ngrid3,igrid3,iodd,rnorm3)
        end do
        write(6,200)'       post3      ',n,rnorm3
c===================================================
c transfer errors to level 2 and correct the error
c===================================================
        call intcrc(delmap3,ngrid3,igrid3,ddelmap2,epsmap2,
     &debmap2,ngrid2,igrid2)
        call addarr(delmap2,ddelmap2,ngrid2,igrid2)
c==========================================
c postsmooth at level 2
c==========================================
        do n = 1,nit2
          do iodd = 0,1
            call getres(delmap2,epsmap2,sumeps2,eff2,
     &    resid2,ngrid2,igrid2,iodd,rnorm2)
            call gausei(delmap2,resid2,rsumeps2,ngrid2,igrid2,
     &        iodd,om2)
          end do
        end do
        do iodd = 0,1
          call getres(delmap2,epsmap2,sumeps2,eff2,
     &  resid2,ngrid2,igrid2,iodd,rnorm2)
        end do
        write(6,200)'     post2      ',n,rnorm2
500   continue
c===================================================
c transfer errors to level 1 and correct the error
c===================================================
        call intcrc(delmap2,ngrid2,igrid2,ddelmap1,
     &epsmap1,debmap1,ngrid1,igrid1)
        call addarr(delmap1,ddelmap1,ngrid1,igrid1)
c==============================
c postsmooth at level 1
c==============================
        do n = 1,nit1
          do iodd = 0,1
            call getres(delmap1,epsmap1,sumeps1,eff1,resid1,
     &    ngrid1,igrid1,iodd,rnorm1)
            call gausei(delmap1,resid1,rsumeps1,ngrid1,igrid1,
     &        iodd,om1)
          end do
        end do
        do iodd = 0,1
          call getres(delmap1,epsmap1,sumeps1,eff1,resid1,
     &  ngrid1,igrid1,iodd,rnorm1)
        end do
        write(6,200)'   post1      ',n,rnorm1
400   continue
c===================================================
c transfer errors to level 0 and correct the potential
c===================================================
        call intcrc(delmap1,ngrid1,igrid1,ddelmap,epsmap,
     &debmap,igrid,igrid)
        call addarr(delmap,ddelmap,igrid,igrid)
c=================================================================
c end of big loop
c=================================================================
        goto 300
301        continue   !here if converged at level 0

        finish = cputime(start)
        call time(hour)
c        write(6,*)'finished ',ntot,' iterations in ',finish
c        print *,' '
        end
c--------------------------------------------------------

        subroutine getres(delmap,epsmap,sumeps,eff,resid,ngrid,
     &    igrid,iodd,rnorm)
c
c calculate residual: res = f - L(func)
c do odd parity points first (iodd = 0->1, so that
c points on coarser grid are left with residuals
c
        real*4 resid(ngrid,ngrid,ngrid)
        real*4 eff(ngrid,ngrid,ngrid)
        real*4 sumeps(ngrid,ngrid,ngrid)
        real*4 delmap(ngrid,ngrid,ngrid)
        real*4 epsmap(ngrid,ngrid,ngrid,3)

        real*8 rnorm
c------------------------------------------------
c        print *,'iodd: ',iodd
        if(iodd.eq.0)rnorm = 0.
        do k = 2,igrid-1
          do j = 2,igrid-1
c iodd = 0: do odd parity points
            ist = mod(iodd+k+j+1,2)+2
            do i = ist,igrid-1,2
                resid(i,j,k) = eff(i,j,k) + sumeps(i,j,k)*
     &  delmap(i,j,k)
     &        - delmap(i-1,j,k)*epsmap(i-1,j,k,1)
     &        - delmap(i+1,j,k)*epsmap(i,j,k,1)
     &        - delmap(i,j-1,k)*epsmap(i,j-1,k,2)
     &        - delmap(i,j+1,k)*epsmap(i,j,k,2)
     &        - delmap(i,j,k-1)*epsmap(i,j,k-1,3)
     &        - delmap(i,j,k+1)*epsmap(i,j,k,3)
                rnorm = rnorm + resid(i,j,k)**2
            end do
          end do
        end do
        return
        end

        subroutine gausei(delmap,resid,rsumeps,ngrid,igrid,iodd,
     &    omega)
c
c update potential estimate using Gauss-Seidel method
c with chequerboard order
c
        real*4 resid(ngrid,ngrid,ngrid)
        real*4 rsumeps(ngrid,ngrid,ngrid)
        real*4 delmap(ngrid,ngrid,ngrid)
c------------------------------------
        do k = 2,igrid-1
          do j = 2,igrid-1
c iodd = 0: do odd parity points
            ist = mod(iodd+k+j+1,2)+2
            do i = ist,igrid-1,2
                delmap(i,j,k) = delmap(i,j,k) - omega*resid(i,j,k)
     &            *rsumeps(i,j,k)
            end do
          end do
        end do
        return
        end

        subroutine doper(phimap,igrid,iper)
c
c if periodic boundary condition option
c force periodicity using wrap around update of boundary values:
c 2nd slice-->last
c last-1 slice-->first
c
        real*4 phimap(igrid,igrid,igrid)
        logical iper(3)
c----------------------------------------
c
c z periodicity
c
        if(iper(3)) then
          do iy = 2,igrid-1
            do ix = 2,igrid-1
                 phimap(ix,iy,igrid) = phimap(ix,iy,2)
                 phimap(ix,iy,1) = phimap(ix,iy,igrid-1)
            end do
          end do
        end if
c
c y periodicity
c
        if(iper(2)) then
          do iy = 2,igrid-1
            do ix = 2,igrid-1
                 phimap(ix,igrid,iy) = phimap(ix,2,iy)
                 phimap(ix,1,iy) = phimap(ix,igrid-1,iy)
            end do
          end do
        end if
c
c x periodicity
c
        if(iper(1)) then
          do iy = 2,igrid-1
            do ix = 2,igrid-1
                 phimap(igrid,ix,iy) = phimap(2,ix,iy)
                 phimap(1,ix,iy) = phimap(igrid-1,ix,iy)
            end do
          end do
        end if
        return
        end

        subroutine transfr(array0,ngrid,igrid,array1,ngrid1,
     &    igrid1,ibound)
c
c transfer points from fine grid to coarse one. if ibound = 0
c then set boundary points to be 0, else transfer these too
c
        real*4 array0(ngrid,ngrid,ngrid)
        real*4 array1(ngrid1,ngrid1,ngrid1)
c------------------------------------------------
        do k = 2,igrid1-1
          k0 = 2*k-1
          do j = 2,igrid1-1
            j0 = 2*j-1
            do i = 2,igrid1-1
                i0 = 2*i-1
c straight transfer
                array1(i,j,k) = array0(i0,j0,k0)
c average
c                array1(i,j,k) = 0.5*array0(i0,j0,k0) + (
c     & array0(i0+1,j0,k0)+array0(i0-1,j0,k0) +
c     & array0(i0,j0+1,k0)+array0(i0,j0-1,k0) +
c     & array0(i0,j0,k0+1)+array0(i0,j0,k0-1))/12.
            end do
          end do
        end do
        if(ibound.eq.0)then
          do k = 1,igrid1,igrid1-1
            do j = 1,igrid1
                do i = 1,igrid1
                  array1(i,j,k) =  0.
                  array1(j,k,i) =  0.
                  array1(k,i,j) =  0.
                end do
            end do
          end do
        else
          do k = 1,igrid1,igrid1-1
            k0 = 2*k-1
            do j = 1,igrid1
                j0 = 2*j-1
                do i = 1,igrid1
                  i0 = 2*i-1
                  array1(i,j,k) =  array0(i0,j0,k0)
                  array1(j,k,i) =  array0(j0,k0,i0)
                  array1(k,i,j) =  array0(k0,i0,j0)
                end do
            end do
          end do
        end if
        return
        end
c
        subroutine hrmav(epsmap,ngrid,igrid,epsmap1,ngrid1,igrid1)
c
c coarse eps is harmonic average of fine eps
c
        real*4 epsmap(ngrid,ngrid,ngrid,3)
        real*4 epsmap1(ngrid1,ngrid1,ngrid1,3)
c------------------------------------------------
        do k = 1,igrid1-1
          k0 = 2*k-1
          do j = 1,igrid1-1
            j0 = 2*j-1
            do i = 1,igrid1-1
                i0 = 2*i-1
                eps1 = epsmap(i0,j0,k0,1)
                eps2 = epsmap(i0+1,j0,k0,1)
                epsmap1(i,j,k,1) = 2.*eps1*eps2/(eps1+eps2)
c debug arithmetic av.
c                epsmap1(i,j,k,1) = (eps1+eps2)/2.
                eps1 = epsmap(i0,j0,k0,2)
                eps2 = epsmap(i0,j0+1,k0,2)
                epsmap1(i,j,k,2) = 2.*eps1*eps2/(eps1+eps2)
c debug arithmetic av.
c                epsmap1(i,j,k,2) = (eps1+eps2)/2.
                eps1 = epsmap(i0,j0,k0,3)
                eps2 = epsmap(i0,j0,k0+1,3)
                epsmap1(i,j,k,3) = 2.*eps1*eps2/(eps1+eps2)
c debug arithmetic av.
c                epsmap1(i,j,k,3) = (eps1+eps2)/2.
            end do
          end do
        end do
        return
        end

        subroutine setarr(array,ngrid,igrid,val)
c
c initialize array to value
c
        real*4 array(ngrid,ngrid,ngrid)
        do k = 1,igrid
          do j = 1,igrid
            do i = 1,igrid
                array(i,j,k) = val
          end do
          end do
        end do
        return
        end

        subroutine addarr(array1,array2,ngrid,igrid)
        real*4 array1(ngrid,ngrid,ngrid)
        real*4 array2(ngrid,ngrid,ngrid)
        do k = 2,igrid-1
          do j = 2,igrid-1
            do i = 2,igrid-1
                array1(i,j,k) = array1(i,j,k)+array2(i,j,k)
            end do
          end do
        end do
        return
        end
c
        subroutine intcrc(delmap1,ngrid1,igrid1,ddelmap,epsmap,
     &debmap,ngrid,igrid)
c
c interpolate errors from coarse to fine grid and correct potential
c
        real*4 delmap1(ngrid1,ngrid1,ngrid1)
        real*4 ddelmap(ngrid,ngrid,ngrid)
        real*4 epsmap(ngrid,ngrid,ngrid,3)
        real*4 debmap(ngrid,ngrid,ngrid)
c------------------------------------------------
c
c common points
c
        do k = 2,igrid1-1
          k0 = 2*k-1
          do j = 2,igrid1-1
            j0 = 2*j-1
            do i = 2,igrid1-1
                i0 = 2*i-1
                ddelmap(i0,j0,k0) = delmap1(i,j,k)
            end do
          end do
        end do
c
c on line points 1st
c
        do k0 = 3,igrid-1,2
          do j0 = 3,igrid-1,2
            do i0 = 2,igrid-1,2
c x dir
                eps1 = epsmap(i0-1,j0,k0,1)
                eps2 = epsmap(i0,j0,k0,1)
                ddelmap(i0,j0,k0) =
     &        (ddelmap(i0-1,j0,k0)*eps1
     &      + ddelmap(i0+1,j0,k0)*eps2)/(eps1+eps2
     &      +debmap(i0,j0,k0))
c y dir
                eps1 = epsmap(j0,i0-1,k0,2)
                eps2 = epsmap(j0,i0,k0,2)
                ddelmap(j0,i0,k0) =
     &        (ddelmap(j0,i0-1,k0)*eps1
     &      + ddelmap(j0,i0+1,k0)*eps2)/(eps1+eps2
     &      +debmap(j0,i0,k0))
c z dir
                eps1 = epsmap(j0,k0,i0-1,3)
                eps2 = epsmap(j0,k0,i0,3)
                ddelmap(j0,k0,i0) =
     &        (ddelmap(j0,k0,i0-1)*eps1
     &      + ddelmap(j0,k0,i0+1)*eps2)/(eps1+eps2
     &      +debmap(j0,k0,i0))
            end do
          end do
        end do
c
c on plane points
c
        do k0 = 3,igrid-1,2
          do j0 = 2,igrid-1,2
            do i0 = 2,igrid-1,2
c xy plane
                eps1 = epsmap(i0-1,j0,k0,1)
                eps2 = epsmap(i0,j0,k0,1)
                eps3 = epsmap(i0,j0-1,k0,2)
                eps4 = epsmap(i0,j0,k0,2)
                ddelmap(i0,j0,k0) =
     &        (ddelmap(i0-1,j0,k0)*eps1 + ddelmap(i0+1,j0,k0)*eps2 +
     &        ddelmap(i0,j0-1,k0)*eps3 + ddelmap(i0,j0+1,k0)*eps4)/
     &        (eps1+eps2+eps3+eps4+debmap(i0,j0,k0))
c yz plane
                eps1 = epsmap(k0,i0-1,j0,2)
                eps2 = epsmap(k0,i0,j0,2)
                eps3 = epsmap(k0,i0,j0-1,3)
                eps4 = epsmap(k0,i0,j0,3)
                ddelmap(k0,i0,j0) =
     &        (ddelmap(k0,i0-1,j0)*eps1 + ddelmap(k0,i0+1,j0)*eps2 +
     &        ddelmap(k0,i0,j0-1)*eps3 + ddelmap(k0,i0,j0+1)*eps4)/
     &        (eps1+eps2+eps3+eps4+debmap(k0,i0,j0))
c zx plane
                eps1 = epsmap(j0,k0,i0-1,3)
                eps2 = epsmap(j0,k0,i0,3)
                eps3 = epsmap(j0-1,k0,i0,1)
                eps4 = epsmap(j0,k0,i0,1)
                ddelmap(j0,k0,i0) =
     &        (ddelmap(j0,k0,i0-1)*eps1 + ddelmap(j0,k0,i0+1)*eps2 +
     &        ddelmap(j0-1,k0,i0)*eps3 + ddelmap(j0+1,k0,i0)*eps4)/
     &        (eps1+eps2+eps3+eps4+debmap(j0,k0,i0))
            end do
          end do
        end do
c
c center points
c
        do k0 = 2,igrid-1,2
          do j0 = 2,igrid-1,2
            do i0 = 2,igrid-1,2
                eps1 = epsmap(i0-1,j0,k0,1)
                eps2 = epsmap(i0,j0,k0,1)
                eps3 = epsmap(i0,j0-1,k0,2)
                eps4 = epsmap(i0,j0,k0,2)
                eps5 = epsmap(i0,j0,k0-1,3)
                eps6 = epsmap(i0,j0,k0,3)
                ddelmap(i0,j0,k0) =
     &        (ddelmap(i0-1,j0,k0)*eps1 + ddelmap(i0+1,j0,k0)*eps2 +
     &        ddelmap(i0,j0-1,k0)*eps3 + ddelmap(i0,j0+1,k0)*eps4 +
     &        ddelmap(i0,j0,k0-1)*eps5 + ddelmap(i0,j0,k0+1)*eps6)/
     &        (eps1+eps2+eps3+eps4+eps5+eps6+debmap(i0,j0,k0))
            end do
          end do
        end do
        return
        end

max_stars_repo_name: caixiuhong/Stable-MCCE

max_issues_repo_name: caixiuhong/Stable-MCCE

max_forks_repo_name: caixiuhong/Stable-MCCE


text:
! Compile and run: make

program phys395_hw2_curve_fitting
  implicit none

  integer, parameter :: ifh = 1
  integer, parameter :: max_lines = 1048576
  integer :: i, stat
  real, dimension(:), allocatable :: x, y
  real, dimension(max_lines) :: x_, y_

  ! Read data
  open(unit=ifh, file="data.dat", action="read", status="old", &
    access="sequential", form="formatted")
  do i = 1, max_lines
    read(ifh, *, iostat=stat) x_(i), y_(i)
    if (stat < 0) exit
  end do
  close(ifh)

  ! Allocate properly-sized arrays
  i = i - 1
  allocate(x(i), y(i))
  x = x_(1:i)
  y = y_(1:i)

  print *
  call write_csv(x, y, "results_svd.csv", solver="svd")
  call write_csv(x, y, "results_lss.csv", solver="lss")

  print "(a)", "The condition number for the SVD matrix is larger than "
  print "(a)", "the condition number for the Linear Least Squares (LSS) matrix."
  print "(a)", "Otherwise, the best fit parameters appear to be the same."
  print *

  deallocate(x, y)

contains

  !! Write csv for specified solver (svd/lss)
  subroutine write_csv(x, y, filename, solver)
    real, dimension(:) :: x, y
    character(len=*) :: filename, solver
    integer, parameter :: ofh = 2
    integer :: i
    real, dimension(size(x)) :: y_fit3, y_fit7

    y_fit3 = fit(x, y, n=3, solver=solver)
    y_fit7 = fit(x, y, n=7, solver=solver)

    open(unit=ofh, file=filename, action="write", status="replace")
    write(ofh, *) "x, y, $f_3(x)$, $f_7(x)$"
    do i = 1, size(x)
      write(ofh, *) x(i), ", ", y(i), ", ", y_fit3(i), ", ", y_fit7(i)
    end do
    close(ofh)
  end subroutine

  !! Fit data to n+1-dimensional basis via specified solver
  function fit(x, y, n, solver) result(ys_fit)
    integer :: a, i, j, k, n
    real, dimension(:) :: x, y
    real, dimension(size(x)) :: ys_fit
    real, dimension(n + 1) :: coeffs, p, s
    real, dimension(size(x), n + 1) :: bxa
    real, dimension(n + 1, n + 1) :: B, U, Vh
    real :: cond_num, chi_actual, chi_expect
    character(len=64) :: fmt_str
    character(len=*) :: solver

    ! Construct bxa to minimize matrix equation | y - bxa x |^2
    ! Construct B and p to solve the matrix equation Bc = p
    forall (a=0:n, i=1:size(bxa, 1)) bxa(i, a + 1) = basis(a, x(i))
    forall (j=1:n+1, k=1:n+1) B(j, k) = sum(bxa(:, j) * bxa(:, k))
    p = matmul(transpose(bxa), y)

    call svd(n + 1, B, U, s, Vh)

    select case (solver)
    case ("svd"); coeffs = solve_svd(p, U, s, Vh, 1.0e-6)
    case ("lss"); call solve_lss(bxa, y, -1.0, coeffs, s)
    end select

    ys_fit = matmul(bxa, coeffs)
    chi_actual = sum((y - ys_fit)**2)
    chi_expect = size(x) - (n + 1)
    cond_num = s(1) / s(n + 1)

    write(fmt_str, "(a, i10, a)") "(", n + 1, "f7.2)"
    print "(a, a)", "solver = ", solver
    print "(a, i1)", "n = ", n
    print *
    write(*, "(a21)", advance="no") "Best fit params: "
    print fmt_str, coeffs
    print "(a21, f10.2)", "Condition number: ",    cond_num
    print "(a21, f10.2)", "Chi-square (actual): ", chi_actual
    print "(a21, f10.2)", "Chi-square (expect): ", chi_expect
    print "(/, a, /)", "----"
  end function

  !! Solve for x minimizing |Ax - B|^2
  !! Results are returned in x, s
  subroutine solve_lss(A, B, rcond, x, s)
    real :: A(:, :), B(:)
    real :: A_(size(A, 1), size(A, 2))
    real :: B_(size(B, 1))
    real :: x(size(A, 2))
    real :: s(min(size(A, 1), size(A, 2)))
    real :: work(6 * max(size(A, 1), size(A, 2)))
    real :: rcond
    integer :: m, n, rank, stat

    m = size(A, 1)
    n = size(A, 2)
    A_ = A
    B_ = B
    call dgelss(m, n, 1, A_, m, B_, m, s, rcond, rank, work, size(work), stat)
    x = B_(1:size(x))
    if (stat /= 0) call abort
  end subroutine

  !! Solve the matrix equation (U s Vh) x = b
  function solve_svd(b, U, s, Vh, eps) result(x)
    real :: b(:), s(:), U(:, :), Vh(:, :), x(size(Vh, 2)), eps

    x = matmul(transpose(U), b)
    where (s > eps * s(1)); x = x / s
    elsewhere;              x = 0.0
    end where
    x = matmul(transpose(Vh), x)
  end function

  !! Decompose A into U, s, Vh
  subroutine svd(n, A, U, s, Vh)
    integer :: n, stat
    real :: A(n, n), A_(n, n), U(n, n), Vh(n, n), s(n), work(6 * n)

    A_ = A
    call dgesvd('A', 'A', n, n, A_, n, s, U, n, Vh, n, work, 6 * n, stat)
    if (stat /= 0) call abort
  end subroutine

  !! An even basis... for a more civilized age
  elemental function basis(a, x)
    integer, intent(in) :: a
    real, intent(in) :: x
    real :: basis
    real, parameter :: pi = 3.14159265358979323846264338327950288419716939937510

    basis = cos(2 * pi * a * x)
  end function

end program phys395_hw2_curve_fitting

max_stars_repo_name: YodaEmbedding/experiments

max_issues_repo_name: YodaEmbedding/experiments

max_forks_repo_name: YodaEmbedding/experiments


text:
program da_thin
!----------------------------------------------------------------------
! Purpose: Thinning wrfinput by using decimation
!
! Input  : wrfinput_hires -- High resolution wrfinput
!
! Output : wrfinput_lores -- Thinned wrfinput
!
! jliu@ucar.edu, 2011-12-15
!----------------------------------------------------------------------

  use netcdf

  implicit none

  integer :: i, n

  integer :: ncidin, ncidout, varid, varid_out, status
  integer :: nDims, nVars, nGlobalAtts, numsAtts, nTimes
  integer :: dLen, attLen, xtype, dID, unlimDimID, TID
  integer :: divided_exactly, dimid

  integer :: dsizes(4), start(4), stride(4)

  integer, dimension(nf90_max_var_dims) :: vDimIDs

  integer, dimension(:),    allocatable :: vdimsizes

  real :: fVal

  real,                 dimension(:,:,:,:),  allocatable :: fVar
  integer,              dimension(:,:,:,:),  allocatable :: iVar
  character (len = 19), dimension(:),        allocatable :: times

  character (len = 14 )           :: coordinates
  character (len = NF90_MAX_NAME) :: vNam, dNam, attNam

  integer             :: decimation_factor = 3
  integer             ::            offset = 2
  character (len=255) ::            filin  = "wrfinput_hires"
  character (len=255) ::            filout = "wrfinput_lores"
  character (len=255) ::               arg = ""
  character (len=255) ::           appname = ""
  character(len=8)    ::            i_char =""

  integer iargc

  call getarg(0, appname)
  n=index(appname, '/', BACK=.true.)
  appname = trim(appname(n+1:))

  DO i = 1, iargc(), 2
    call getarg(i, arg)
    select case ( trim(arg) )
      case ("-i")
        call getarg(i+1, arg)
        filin=trim(arg)
      case ("-o")
        call getarg(i+1, arg)
        filout=trim(arg)
      case ("-thin")
        call getarg(i+1, arg)
        read(arg, '(i3)') decimation_factor
      case default
        Write(*,*) "Usage : "//trim(appname)//" [-i inputfile] [-o outputfile] [-thin decimation_factor] [-h]"
        Write(*,*) "  -i     Optional, input filename,    default - wrfinput_hires"
        Write(*,*) "  -o     Optional, output filename,   default - wrfinput_lores"
        Write(*,*) "  -thin  Optional, decimation factor, default - 3"
        Write(*,*) "  -h     Shwo this usage"
        call exit(0)
    end select
  END DO

  if ( mod(decimation_factor,2) == 0 ) then
    Write(*,*) "\nError : decimation factor must be odd number\n"
    call exit(-1)
  endif

  status = nf90_open(filin, NF90_NOWRITE, ncidin)
  if ( status /= nf90_noerr ) then
    Write (*,*) "File open error. Please link the input file to "//trim(filin)
    call exit(-1)
  endif

  status = nf90_inq_dimid(ncidin, "west_east_stag", dimid)
  status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
  divided_exactly = mod((dLen-1),decimation_factor)

  status = nf90_inq_dimid(ncidin, "south_north_stag", dimid)
  status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
  divided_exactly = divided_exactly + mod((dLen-1),decimation_factor)

  if ( divided_exactly /= 0 ) then
    Write (*,fmt='(a,i2,a)') "Failed to thinning. Grids need to match : ( n - 1 ) mod ",decimation_factor," = 0"
    call exit(-1)
  endif

  status = nf90_create(filout, NF90_CLOBBER, ncidout)
  if ( status /= nf90_noerr) call nf90_handle_err(status)

  status = nf90_inquire(ncidin, nDims, nVars, nGlobalAtts, unlimDimID)
  if ( status /= nf90_noerr ) call nf90_handle_err(status)

  write (i_char, '(i8)') decimation_factor

  Write (*,*) "       Input file : "//trim(filin)
  Write (*,*) "      Output file : "//trim(filout)
  Write (*,*) "decimation factor : "//adjustl(i_char)

  do i=1, nGlobalAtts
    status = nf90_inq_attname(ncidin, NF90_GLOBAL, i, attNam)
    select case (trim(attNam))
      case ( "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION", &
             "WEST-EAST_PATCH_END_UNSTAG", "WEST-EAST_PATCH_END_STAG", &
             "SOUTH-NORTH_PATCH_END_UNSTAG", "SOUTH-NORTH_PATCH_END_STAG" )
        status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
        status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, int(( fVal - 1 ) / decimation_factor + 1) )
      case ("DX","DY", "DT" )
        status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
        status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, fVal * decimation_factor )
      case default
        status = nf90_copy_att(ncidin, NF90_GLOBAL, attNam, ncidout, NF90_GLOBAL)
      end select
  end do

  allocate (vdimsizes(nDims), stat=status)

  do i=1, nDims

    status = nf90_inquire_dimension(ncidin, i, name=dNam, len = dLen)

    vdimsizes(i) = dLen
    select case (trim(dNam))
      case ("south_north_stag", "west_east_stag")
        vdimsizes(i) = (dLen - 1 ) / decimation_factor + 1
      case ("west_east", "south_north")
        vdimsizes(i) = dLen / decimation_factor
      case ("Time")
        allocate(times(dLen), stat=status)
        vdimsizes(i) = NF90_UNLIMITED
        nTimes = dLen
        TID = i
    end select

    status = nf90_def_dim(ncidout, dNam, vdimsizes(i), dID)

  end do

  vdimsizes(TID) = nTimes

  do varid=1, nVars
    status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs,natts=numsAtts)
    status = nf90_def_var(ncidout, trim(vNam), xtype, vDimIDs(1:nDims), varid_out)
    if(status /= nf90_NoErr) call nf90_handle_err(status)
    do i=1, numsAtts
      status = nf90_inq_attname(ncidin, varid, i, attNam)
      status = nf90_copy_att(ncidin, varid, trim(attNam), ncidout, varid_out)
      if(status /= nf90_NoErr) call nf90_handle_err(status)
    end do
  end do

  status = nf90_enddef(ncidout)

  offset = (decimation_factor + 1) / 2

  do varid=1, nVars

    status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs)

    dsizes = 1
    do i = 1 , nDims
      dsizes(i) = vdimsizes(vDimIDs(i))
    end do

    status = nf90_inquire_attribute(ncidin,varid,"coordinates")

    if ( status == nf90_noerr ) then

      Write(*,*) "Thinning for "//trim(vNam)

      coordinates=char(0)
      status = nf90_get_att(ncidin, varid, "coordinates" , coordinates)
      !print *, coordinates

      stride=(/decimation_factor,decimation_factor,1,1/)

      n = index(coordinates, char(0)) - 1
      if ( n < 0 ) n = len(coordinates)

      select case (trim(coordinates(1:n)))
        case ("XLONG_V XLAT_V")
          start=(/offset,1,1,1/)
        case ("XLONG_U XLAT_U")
          start=(/1,offset,1,1/)
        case ("XLONG XLAT")
          start=(/offset,offset,1,1/)
        case ("XLONG XLAT XTI")
          start=(/offset,offset,1,1/)
        case default
          print *, "Unkown coordinates : "//coordinates
          call exit(-1)
      end select

    else

      stride = 1
      start  = 1

      if ( trim(vNam) == 'XLONG' .or. trim(vNam) == 'XLAT' ) then
          stride = (/decimation_factor,decimation_factor,1,1/)
          start  = (/offset,offset,1,1/)
      endif

    endif

    select case (xtype)
      case (nf90_float)
        allocate(fVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
        status = nf90_get_var(ncidin, varid, fVar, start=start, stride=stride)
        if ( status == nf90_noerr ) then
          if ( vNam == "RDX" .or. vNam == "RDY" ) then
            status = nf90_put_var(ncidout, varid, fVar / decimation_factor)
          else
            status = nf90_put_var(ncidout, varid, fVar)
          endif
          if ( status /= nf90_noerr) call nf90_handle_err(status)
        else
          call nf90_handle_err(status)
        endif
        deallocate(fVar, stat=status)
      case (nf90_int)
        allocate(iVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
        status = nf90_get_var(ncidin, varid, iVar, start=start, stride=stride)
        if ( status == nf90_noerr ) then
          status = nf90_put_var(ncidout, varid, iVar)
          if ( status /= nf90_noerr) call nf90_handle_err(status)
        else
          call nf90_handle_err(status)
        endif
        deallocate(iVar, stat=status)
      case (nf90_char)
        if ( trim(vNam) == "Times") then
          status = nf90_get_var(ncidin, varid, times)
          if ( status == nf90_noerr ) then
            status = nf90_put_var(ncidout, varid, times)
            if ( status /= nf90_noerr) call nf90_handle_err(status)
          else
            call nf90_handle_err(status)
          endif
          deallocate(times, stat=status)
        else
          print *, "Unkown character variable :"//trim(vNam)
          call exit(-1)
        endif
      case default
        print *, "Unkown xtype : ", xtype
        call exit(-1)
    end select
  end do

  status = nf90_close(ncidin)
  status = nf90_close(ncidout)

  Write(*,*) "Completed thinning successfully"

contains

  subroutine nf90_handle_err(status)
    integer, intent ( in) :: status

    if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
      call exit(-1)
    end if
  end subroutine nf90_handle_err

end program da_thin

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: matzegoebel/WRF-fluxavg

max_forks_repo_name: matzegoebel/WRF-fluxavg


text:
C
C $Id: tdgetr.f,v 1.5 2008-07-27 00:17:32 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE TDGETR (PNAM,RVAL)
C
      CHARACTER*(*) PNAM
C
C The subroutine TDGETR may be used to get TDPACK parameters which have
C values of type REAL.
C
C The variables in the following common block define the mapping from
C 3-space to 2-space.
C
        COMMON /TDCOM1/ IH,IT,XM,YM,ZM,XO,YO,ZO,XT,YT,ZT,OE,XE,YE,ZE
        COMMON /TDCOM1/ A1,B1,C1,D1,E1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
        COMMON /TDCOM1/ IS,FV,VL,VR,VB,VT,WL,WR,WB,WT
        SAVE   /TDCOM1/
C
C The variables in the following common block define the parallelogram.
C
        COMMON /TDCOM2/ XACP,YACP,ZACP,XCDX,YCDX,ZCDX,XCDY,YCDY,ZCDY
        SAVE   /TDCOM2/
C
C The variables in the following common block together determine how
C big characters written by TDLBLA will be.
C
        COMMON /TDCOM4/ CSM1,CSM2
        SAVE   /TDCOM4/
C
C The variables in the following common block define the shading type
C and the position of the light source.
C
        COMMON /TDCOM6/ ISHD,XPLS,YPLS,ZPLS
        SAVE   /TDCOM6/
C
C Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
        CALL TDBLDA
C
C Check for an uncleared prior error.
C
        IF (ICFELL('TDGETR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Get the selected parameter.
C
        IF      (PNAM(1:3).EQ.'CS1'.OR.PNAM(1:3).EQ.'cs1') THEN
          RVAL=CSM1
        ELSE IF (PNAM(1:3).EQ.'CS2'.OR.PNAM(1:3).EQ.'cs2') THEN
          RVAL=CSM2
        ELSE IF (PNAM(1:3).EQ.'FOV'.OR.PNAM(1:3).EQ.'fov') THEN
          RVAL=FV
        ELSE IF (PNAM(1:3).EQ.'HND'.OR.PNAM(1:3).EQ.'hnd') THEN
          RVAL=REAL(IH)
        ELSE IF (PNAM(1:3).EQ.'LSU'.OR.PNAM(1:3).EQ.'lsu') THEN
          RVAL=XPLS
        ELSE IF (PNAM(1:3).EQ.'LSV'.OR.PNAM(1:3).EQ.'lsv') THEN
          RVAL=YPLS
        ELSE IF (PNAM(1:3).EQ.'LSW'.OR.PNAM(1:3).EQ.'lsw') THEN
          RVAL=ZPLS
        ELSE IF (PNAM(1:3).EQ.'SET'.OR.PNAM(1:3).EQ.'set') THEN
          RVAL=REAL(IS)
        ELSE IF (PNAM(1:3).EQ.'SHD'.OR.PNAM(1:3).EQ.'shd') THEN
          RVAL=REAL(ISHD)
        ELSE IF (PNAM(1:3).EQ.'STE'.OR.PNAM(1:3).EQ.'ste') THEN
          RVAL=REAL(IT)
        ELSE IF (PNAM(1:3).EQ.'VPB'.OR.PNAM(1:3).EQ.'vpb') THEN
          RVAL=VB
        ELSE IF (PNAM(1:3).EQ.'VPL'.OR.PNAM(1:3).EQ.'vpl') THEN
          RVAL=VL
        ELSE IF (PNAM(1:3).EQ.'VPR'.OR.PNAM(1:3).EQ.'vpr') THEN
          RVAL=VR
        ELSE IF (PNAM(1:3).EQ.'VPT'.OR.PNAM(1:3).EQ.'vpt') THEN
          RVAL=VT
        ELSE
          CALL SETER ('TDGETR - UNRECOGNIZED PARAMETER NAME',2,1)
          RETURN
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      subroutine spinorcurr(N,p,za,zb,zab,zba)
c---Calculate spinor products
c---extended to deal with negative energies ie with all momenta outgoing
c---Arbitrary conventions of Bern, Dixon, Kosower, Weinzierl,
c---za(i,j)*zb(j,i)=s(i,j)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      double precision p(mxpart,4),rt(mxpart)
      double complex c23(mxpart),f(mxpart),fac,
     & zab(mxpart,4,mxpart),zba(mxpart,4,mxpart)
      integer i,j,N

c---if one of the vectors happens to be zero this routine fails.
      do j=1,N
         za(j,j)=czip
         zb(j,j)=za(j,j)
         zab(j,:,j)=dcmplx(2d0*p(j,:))
         zba(j,:,j)=zab(j,:,j)
C-----positive energy case
         if (p(j,4) .gt. 0d0) then
            rt(j)=dsqrt(p(j,4)+p(j,1))
            c23(j)=dcmplx(p(j,3),-p(j,2))
            f(j)=cone
         else
C-----negative energy case
            rt(j)=dsqrt(-p(j,4)-p(j,1))
            c23(j)=dcmplx(-p(j,3),p(j,2))
            f(j)=im
         endif
      enddo

      do i=2,N
         do j=1,i-1
         s(i,j)=two*(p(i,4)*p(j,4)-p(i,1)*p(j,1)
     &              -p(i,2)*p(j,2)-p(i,3)*p(j,3))
         za(i,j)=f(i)*f(j)
     &   *(c23(i)*dcmplx(rt(j)/rt(i))-c23(j)*dcmplx(rt(i)/rt(j)))

      fac=f(i)*f(j)
C      elseif(h .eq. +1) then
      zba(i,3,j)=fac
     . *(dcmplx(rt(i)/rt(j))*c23(j)+dcmplx(rt(j)/rt(i))*dconjg(c23(i)))
      zba(i,2,j)=fac*im
     . *(dcmplx(rt(i)/rt(j))*c23(j)-dcmplx(rt(j)/rt(i))*dconjg(c23(i)))
      zba(i,1,j)=fac
     . *(dcmplx(rt(i)*rt(j))-dconjg(c23(i))*c23(j)/dcmplx(rt(i)*rt(j)))
      zba(i,4,j)=fac
     . *(dcmplx(rt(i)*rt(j))+dconjg(c23(i))*c23(j)/dcmplx(rt(i)*rt(j)))

C      elseif(h .eq. -1) then
      zab(i,3,j)=+fac
     . *(dcmplx(rt(i)/rt(j))*dconjg(c23(j))+dcmplx(rt(j)/rt(i))*c23(i))
      zab(i,2,j)=fac*im
     . *(dcmplx(rt(j)/rt(i))*c23(i)-dcmplx(rt(i)/rt(j))*dconjg(c23(j)))
      zab(i,1,j)=fac
     . *(dcmplx(rt(i)*rt(j))-c23(i)*dconjg(c23(j))/dcmplx(rt(i)*rt(j)))
      zab(i,4,j)=fac
     . *(dcmplx(rt(i)*rt(j))+c23(i)*dconjg(c23(j))/dcmplx(rt(i)*rt(j)))

      if (abs(s(i,j)).lt.1d-5) then
      zb(i,j)=-(f(i)*f(j))**2*dconjg(za(i,j))
      else
      zb(i,j)=-dcmplx(s(i,j))/za(i,j)
      endif

      za(j,i)=-za(i,j)
      zb(j,i)=-zb(i,j)
      s(j,i)=s(i,j)
      zab(j,:,i)=zba(i,:,j)
      zba(j,:,i)=zab(i,:,j)
      enddo
      enddo
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      program chaine05

C     Fortran does not take trailing SPACEs into account when comparing
C     strings because any string is assumed padded with SPACE up to its
C     declared length

      character*10 s
      character*10 t
      character*10 u
      character*20 v
      character*20 w

      s = "hello"
      t = "hello "

      if(s.eq.t) then
!        FALSE
         print *, "Trailing blanks are ignored"
      else
!        TRUE
         print *, "Trailing blanks are not ignored"
      endif

      s = "hello"
      t = " hello"

      if(s.eq.t) then
!        FALSE
         print *, "Leading blanks are ignored"
      else
!        TRUE
         print *, "Leading blanks are not ignored"
      endif

      s = "hello"
      t = "hello "
      u = "world"
      v = s // u
      w = t // u

      if(v.eq.w) then
!        FALSE
         print *, "Trailing blanks are ignored when concatenating"
      else
!        TRUE because s and t are declared with the same length and hence padding
         print *, "Trailing blanks are not ignored when concatenating"
      endif

      print *, s, t, u, v, w

      end


text:
Tutormatcher was started in 2003 in Davis, CA under the name Aggie Tutors by Craig Shacklett, a UC Davis Sophomore. Craig had tutored a variety of students and subjects while attending high school in the bay area and had earned a reputation from local teachers as someone that could dramatically improve students’ grades and train them to become exceptional students.  Aggie Tutors still operates in the Davis region with support from Tutormatcher.

Parents loved having a young tutor/mentor that their children could look to as a role model. Referrals were pouring in and soon Craig had to start turning away potential clients. That’s when he had the idea, “if I hire friendly, likable people and train them to tutor in the same techniques that I used, there would be no limit to how many students we could help.”

Using a stringent selection process, Craig put together a staff of incredibly talented college students and trained them to be exceptional tutors. Using only wordofmouth advertising, Craig’s new tutoring service expanded to one of the most wellknown and respected in the area.
We feel that we have perfected the formula for helping our clients become the successful and ambitious students that their parents know they can be. That is why we are bringing this opportunity to your area. Take advantage of your local resources and help your child today!

Other Information:
Hire UC Davis students with cumulative GPAs of no less than 3.5
Match Tutors to Staff via Free initial interview with client
Tutor Clients of all ages, Kindergarten to College

Rates
25$30$/Hour for oneonone in home
Group Rate varies from 2510$/hour per person
Discounts Available during Summer
College Student Discount Available for those who qualify
Referral Discounts

   Need help with classwork?  Davis Wiki has a list of tutoring tutoring resources!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C     =================
      SUBROUTINE GETLUMI
C     =================
C     This routine reconstructs the energy and position of deposits
C     in the Lumi calorimeter (electron and photon), the BPC,
C     the 8 meter tagger and the 44 meter tagger.
C
      IMPLICIT NONE

#include "partap.inc"
#include "lmresu.inc"
#include "lmposd.inc"
#include "lmeb.inc"
c#include "t8_data.inc"
#include "common.inc"
C
      INTEGER  NO_USE, I,J,K, Ierr, Irun, JBYT
      EXTERNAL JBYT
      REAL     XPOS, ERR_X, YPOS, ERR_Y
      DATA     NO_USE / -999.99 /

C =========================
C === LUMI calorimeters ===
C =========================
      TEMPLUMG = -500
      TEMPLUME = -500
      ENE44M = 1000
      POSX_LG = NO_USE
      POSY_LG = NO_USE
      POSX_LE = NO_USE
      POSY_LE = NO_USE

      ENE_LG  = NO_USE
      ENE_LE  = NO_USE

C =======================
C === 44 meter tagger ===
C =======================
      TAG44_E  = 0
c      TAG44_X1 = 0
c      TAG44_X2 = 0
c      TAG44_X3 = 0
      CALL FETTAB(LMEB, ID, 1)
      TAG44_E =255-jbyt(LMEB_Rexy(8),25,8)

C -------------------
C --- LUMI ENERGY ---
C -------------------
      IF (COUTAB(LMRESU).GT.0) THEN
          CALL FETTAB(LMRESU,ID,1)
          ENE_LG = LMRESU_Enrg
          ENE_LE = LMRESU_Enre
c       else
c          ene_lg = -999
c          ene_le = -999
      ENDIF

C ----------------
C --- POSITION ---
C ----------------
	xpos = -1000.
	ypos = -1000.
      CALL LMPOSE(xpos,err_x,ypos,err_y,Ierr)
      IF (Ierr .GT. 0) THEN
         POSX_LE = XPOS
         POSY_LE = YPOS
      ENDIF

	xpos = -1000.
	ypos = -1000.
      CALL LMPOSG(xpos,err_x,ypos,err_y,Ierr)
      IF (Ierr .GT. 0) THEN
         POSX_LG = XPOS
         POSY_LG = YPOS
      ENDIF

C ======================
C === 8 meter tagger ===
C ======================
c      TAG8_X = NO_USE
c      TAG8_E = NO_USE
c      TAG8_P = NO_USE
c      TAG8_V = NO_USE

c      Irun = RUN_NUM
c      CALL GT_TAG8(irun,ierr)  !!! routine by Ulrike Wollmer

c      IF (Ierr.GE.0) THEN

c        TAG8_X = t8_x
c        TAG8_E = t8_e
c        TAG8_P = t8_pres
c        TAG8_V = t8_veto

c      ENDIF

C =======================
C === 44 meter tagger ===
C =======================
c      TAG44_E  = 0
c      TAG44_X1 = 0
c      TAG44_X2 = 0
c      TAG44_X3 = 0

c      IF (COUTAB(LMEB).GT.0) THEN

c      CALL FETTAB(LMEB, ID, 1)
c      TAG44_E =255-jbyt(LMEB_Rexy(8),25,8) ! ene deposit in cal

c      TAG44_X1=255-jbyt(LMEB_Rexy(8),17,8) ! ene deposit in scinti.
C                                            finger close to the bpipe

c      TAG44_X2=255-jbyt(LMEB_Rexy(8),9,8)  ! ene deposit in scinti.
C                                            finger ... next one ..

c      TAG44_X3=255-jbyt(LMEB_Rexy(8),1,8)  ! ene deposit in scinti
C                                            finger far away from bpipe
c      ENDIF

c      call fettab(LMEB,ID,1)
      RETURN
      END

ext: fpp


text:
C
C $Id: mplndr.f,v 1.14 2008-07-27 00:17:07 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MPLNDR (FLNM,ILVL)
        CHARACTER*(*) FLNM
        INTEGER       ILVL
        IF (ICFELL('MPLNDR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
        CALL MDLNDR (FLNM,ILVL)
        IF (ICFELL('MPLNDR',2).NE.0) RETURN
        RETURN
      END