text:
      subroutine getVrhoL(AL,vion_pL,vionT_pL,xatom,
     &   ntype,iatom,rhocr_pL,totNel,ido_rho,islda)

******************************************
cc     Written by Lin-Wang Wang, March 30, 2001.
*************************************************************************
**  copyright (c) 2003, The Regents of the University of California,
**  through Lawrence Berkeley National Laboratory (subject to receipt of any
**  required approvals from the U.S. Dept. of Energy).  All rights reserved.
*************************************************************************

******************************************

ccccc if ido_rho.eq.1, generate rho_n from atom,
ccccc if ido_rho.eq.0, do not touch rho_n, rho_n is inputted.

      use fft_data
      use load_data
      use data

      implicit double precision (a-h,o-z)

      include 'param.escan_real'
      include 'mpif.h'

      real*8 xatom(3,matom),AL(3,3)

      real*8 vion_pL(mr_nL)
      real*8 vionT_pL(mr_nL),rhocr_pL(mr_nL)

      real*8 qi2(mnq),vq(mnq,mtype),rhoq(mnq,mtype)
      real*8 vqT(mnq,mtype),rhocq(mnq,mtype)
      real*8 occ_t(mtype)
c
      integer iiatom(mtype),iatom(matom),icore(mtype),numref(matom),
     &  ityatom(matom)

      integer icoul
      real*8 xcoul(3)

      character*20 file_tmp

      complex*16 cc

      real*8, allocatable, dimension(:)  :: workr_nL
      real*8, allocatable, dimension(:)  :: workr_nL2,vion_pL2

***************************************************
****  xatom(1),xatom(2),xatom(3) are the coord in unit of AL(3,3)
****  supercell edges
***************************************************

      common /comNL2/occ_t,iiatom,icore,numref,ityatom
      common /comVrho/qi2,vq,rhoq,vqT,rhocq
      common /comcoul/icoul,xcoul

*******************************************************
**** generate the Kleiman-Bylander reference wavefunction
*******************************************************
      allocate(workr_nL(mr_nL))

      icc=0
      do itype=1,ntype
      icc=icc+icore(itype)
      enddo

      vins=1.d0/vol

      nh1L=(n1L+1)/2+1
      nr_nL=n1L*n2L*n3L/nnodes
      nrL=n1L*n2L*n3L

      ng2_nL=ngtotnod2L(inode)

      vion_pL = 0.0d0
      vionT_pL = 0.0d0
      rhocr_pL = 0.0d0

      if(ido_rho.eq.1) rho_nL = 0.0d0

      do 10 i=1,ng2_nL
      y=0.d0

      do 9  itype=1,ntype
        cc=dcmplx(0.d0,0.d0)
        do ia=1,natom

	if(iatom(ia).eq.iiatom(itype)) then
        x1=xatom(1,ia)
        y1=xatom(2,ia)
        z1=xatom(3,ia)

        x11=AL(1,1)*x1+AL(1,2)*y1+AL(1,3)*z1
        y11=AL(2,1)*x1+AL(2,2)*y1+AL(2,3)*z1
        z11=AL(3,1)*x1+AL(3,2)*y1+AL(3,3)*z1

        ph=gkx2_nL(i)*x11+gky2_nL(i)*y11+gkz2_nL(i)*z11
        cc=cc+cdexp(dcmplx(0.d0,ph))
	endif

        enddo

      q=dsqrt(gkx2_nL(i)**2+gky2_nL(i)**2+gkz2_nL(i)**2)

      iq=1+q*(mnq-1.d0)/qi2(mnq)

      x=(q-qi2(iq))/(qi2(iq+1)-qi2(iq))    ! assuming equal distance grid

cccccc for even smoother treatment, we should use four points interpolation
cccccc But that might not be necessary if we use mqline=4000

      f1=1-x-0.5d0*x*(1-x)
      f2=x+x*(1-x)
      f3=-0.5d0*x*(1-x)       ! using quadratic interpolation

      if(icoul.eq.0) then
      y=vq(iq,itype)*qi2(iq)**2*f1+vq(iq+1,itype)*qi2(iq+1)**2*f2+
     &   vq(iq+2,itype)*qi2(iq+2)**2*f3
      endif

      y2=vqT(iq,itype)*qi2(iq)**2*f1+vqT(iq+1,itype)*qi2(iq+1)**2*f2
     &   +vqT(iq+2,itype)*qi2(iq+2)**2*f3

      if(q.lt.1.D-6) then
      y=0.d0
      y2=vqT(4,itype)
      else
      y=y/q**2
      y2=y2/q**2
      endif

      y1=rhoq(iq,itype)*f1+rhoq(iq+1,itype)*f2
     &  +rhoq(iq+2,itype)*f3        ! rhoq is the full charge, not just psi^2

      i2 = i*2

      if(icoul.eq.0) then
      vion_pL(i2-1)=vion_pL(i2-1)+y*dreal(cc)*vins
      vion_pL(i2)=vion_pL(i2)+y*dimag(cc)*vins
      endif

      vionT_pL(i2-1)=vionT_pL(i2-1)+y2*dreal(cc)*vins
      vionT_pL(i2)=vionT_pL(i2)+y2*dimag(cc)*vins

      if(ido_rho.eq.1) then
      rho_nL(i2-1,1)=rho_nL(i2-1,1)+y1*dreal(cc)*vins
      rho_nL(i2,1)=rho_nL(i2,1)+y1*dimag(cc)*vins
      endif

      if(icore(itype).ne.0) then
      y1=rhocq(iq,itype)*f1+rhocq(iq+1,itype)*f2+
     &      rhocq(iq+2,itype)*f3
      rhocr_pL(i2-1)=rhocr_pL(i2-1)+y1*dreal(cc)*vins
      rhocr_pL(i2)=rhocr_pL(i2)+y1*dimag(cc)*vins
      endif

  9   continue
 10   continue

      scale=1.d0

      if(ido_rho.eq.1) then
      call d3fft_real2L(rho_nL(1,1),workr_nL,-1,0)
      if(islda.eq.1) rho_nL(:,1) = workr_nL     ! rho_nL is the full charge
      if(islda.eq.2) then
      rho_nL(:,1)=workr_nL*0.6d0
      rho_nL(:,2)=workr_nL*0.4d0
      endif
      endif

      if(icoul.eq.0) then
      call d3fft_real2L(vion_pL,workr_nL,-1,0)
      vion_pL = workr_nL
      endif

*************************************

      call d3fft_real2L(vionT_pL,workr_nL,-1,0)
      vionT_pL = workr_nL

*************************************
      if(icc.gt.0) then
      call d3fft_real2L(rhocr_pL,workr_nL,-1,0)

      s=0.d0
      s0=0.d0
      do i=1,nr_nL
      rhocr_pL(i) = workr_nL(i)
      s0=s0+rhocr_pL(i)
      if(rhocr_pL(i).lt.0.d0) rhocr_pL(i)=0.d0
      s=s+rhocr_pL(i)
      enddo
      call global_sumr(s)
      call global_sumr(s0)
      s=s*vol/nrL
      s0=s0*vol/nrL
      if(inode_tot.eq.1) then
      write(6,*) "The total core charge for core correction",s,s0
      write(22,*) "The total core charge for core correction",s,s0
      endif
      s=s0/s
      do i=1,nr_nL
      rhocr_pL(i)=rhocr_pL(i)*s
      enddo

      else
      rhocr_pL=0.d0
      endif
*************************************
      s=0.d0
      do iislda=1,islda
      do i=1,nr_nL
      s=s+rho_nL(i,iislda)
      enddo
      enddo

      call global_sumr(s)

      s=s*vol/nrL

      if(inode_tot.eq.1) then
      write(6,*) "input charge sum=",s
      write(6,*) "correct that to totNel", totNel
      endif

      s=totNel/s               ! rho_nL is the full charge, not just \sum psi^2(r)
      do iislda=1,islda
      do i=1,nr_nL
      rho_nL(i,iislda)=rho_nL(i,iislda)*s        ! do not use abs(),
      enddo
      enddo

      deallocate(workr_nL)

ccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccc

      if(icoul.eq.0) return
cccc calculate vion_pL for icoul.ne.0
cccc for icoul.ne.0, vion_pL must be calculated based on n1L2,n2L2,n3L2,AL2

      vins2=1.d0/vol2
      ng2_nL2=ngtotnod2L2(inode)

      allocate(workr_nL2(mr_nL2))
      allocate(vion_pL2(mr_nL2))

      vion_pL2=0.d0

      do 110 i=1,ng2_nL2

      do 19  itype=1,ntype
        cc=dcmplx(0.d0,0.d0)
        do ia=1,natom

	if(iatom(ia).eq.iiatom(itype)) then   ! even for icoul=11,12,13, shifting x1,y1,z1 is okay
        x1=xatom(1,ia)
         if(x1.gt.1.d0) x1=x1-1.d0
         if(x1.lt.0.d0) x1=x1+1.d0
         if(x1.gt.xcoul(1)) x1=x1-1.d0
        y1=xatom(2,ia)                      ! actually x2, not really y
         if(y1.gt.1.d0) y1=y1-1.d0
         if(y1.lt.0.d0) y1=y1+1.d0
         if(y1.gt.xcoul(2)) y1=y1-1.d0
        z1=xatom(3,ia)                      ! actually x3, not really z
         if(z1.gt.1.d0) z1=z1-1.d0
         if(z1.lt.0.d0) z1=z1+1.d0
         if(z1.gt.xcoul(3)) z1=z1-1.d0

         if(icoul.eq.1) then
         x1=x1/2
         y1=y1/2
         z1=z1/2
         endif
         if(icoul.eq.11) then
         x1=x1/2
         endif
         if(icoul.eq.12) then
         y1=y1/2
         endif
         if(icoul.eq.13) then
         z1=z1/2
         endif

        x11=AL2(1,1)*x1+AL2(1,2)*y1+AL2(1,3)*z1
        y11=AL2(2,1)*x1+AL2(2,2)*y1+AL2(2,3)*z1
        z11=AL2(3,1)*x1+AL2(3,2)*y1+AL2(3,3)*z1

        ph=gkx2_nL2(i)*x11+gky2_nL2(i)*y11+gkz2_nL2(i)*z11
        cc=cc+cdexp(dcmplx(0.d0,ph))
	endif

        enddo

      i2 = i*2

      vion_pL2(i2-1)=vion_pL2(i2-1)+
     &        vcoul_nL2(i2-1,itype)*dreal(cc)-
     &        vcoul_nL2(i2,itype)*dimag(cc)
      vion_pL2(i2)=vion_pL2(i2)+
     &        vcoul_nL2(i2-1,itype)*dimag(cc)+
     &        vcoul_nL2(i2,itype)*dreal(cc)

 19   continue

110   continue

      call d3fft_real2L2(vion_pL2,workr_nL2,-1,0)
      vion_pL2=workr_nL2

      deallocate(workr_nL2)

      call convert_2LtoL(vion_pL,vion_pL2,2)

      deallocate(vion_pL2)

      return
cccccccccccccccccccccccccccccccccccccccc

      end

ext: f

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
! Issue # 96: https://github.com/OpenFortranProject/open-fortran-parser/issues/96
!
! prefix-spec was not updated for F2008 syntax. For example
!
!
impure subroutine impure_proc
end subroutine

ELEMENTAL subroutine elemental_proc
end subroutine

MODULE subroutine module_proc
end subroutine

PURE subroutine pure_proc
end subroutine

RECURSIVE subroutine recursive_proc
end subroutine

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
	SUBROUTINE JTB_PDMP ( lunout, seqnm, parms, np, iret )
C************************************************************************
C* JTB_PDMP								*
C*									*
C* This subroutine writes out information about each parameter in	*
C* parms.								*
C*									*
C* JTB_PDMP ( LUNOUT, SEQNM, PARMS, NP, IRET )				*
C*									*
C* Input parameters:							*
C*	LUNOUT		INTEGER		Open output unit number		*
C*	SEQNM		CHAR*		Sequence mnemonic		*
C*	PARMS (NP)	CHAR*8		Parameter names			*
C*	NP		INTEGER		# of parameter names		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-10 = wrtie error		*
C**									*
C* Log:									*
C* K. Brill/EMC		 9/98						*
C* T. Lee/GSC		12/98	Fixed FORMAT statement, X -> 1X		*
C************************************************************************
	INCLUDE		'jtb.cmn'
C*
	CHARACTER*(*)	seqnm, parms (*)
C*
	CHARACTER*40	descr
	CHARACTER*8	prm
	CHARACTER*12	unit
	LOGICAL		fnddsc, fndpck
C-----------------------------------------------------------------------
	iret = 0
	WRITE (lunout, 1001, IOSTAT=ios) seqnm
1001	FORMAT ( /, ' FOR SEQUENCE ', A, ':', / )
	IF ( ios .ne. 0 ) THEN
	    iret = -10
	    RETURN
	END IF
	IF ( np .eq. 0 ) RETURN
	WRITE (lunout, 1002)
1002	FORMAT ( 3X, 'NAME', 12X, 'PARAMETER DESCRIPTION', 11X, 'SCL',
     +          4X, 'REF', 3X, 'BIT', 4X, 'UNITS' )
1003	FORMAT ( 1X, A8, 1X, A40, 1X, I3, 1X, I8, 1X, I3, 1X, A12 )
	DO ip = 1, np
	    fnddsc = .false.
	    jdesc = 0
	    CALL ST_LSTR ( parms (ip), l1, ier )
	    DO WHILE ( .not. fnddsc .and. jdesc .lt. nmbprm )
		jdesc = jdesc + 1
		CALL ST_LSTR ( tbbprm (jdesc), l2, ier )
		fnddsc = ( parms (ip) (1:l1) .eq.
     +			   tbbprm (jdesc) (1:l2) )
	    END DO
	    IF ( .not. fnddsc ) THEN
		descr = ' NOT FOUND'
	    ELSE
		descr = ' NO DESCRIPTION AVAILABLE'
	    END IF
	    fndpck = fnddsc
	    jpack = jdesc
	    IF ( .not. fndpck ) THEN
		iscl = 0
		iref = 0
		ibit = 0
		unit = ' '
	    ELSE
		iscl = ipkscl (jpack)
		iref = ipkref (jpack)
		ibit = ipkbit (jpack)
		unit = pkunit (jpack)
	    END IF
	    prm = parms (ip)
	    CALL ST_LSTR ( prm, lng, ier )
	    IF ( lng .lt. 8 ) THEN
		DO iq = lng+1, 8
		    prm (iq:iq) = ' '
		END DO
	    END IF
	    CALL ST_LSTR ( descr, lng, ier )
	    IF ( lng .lt. 40 ) THEN
		DO iq = lng+1, 40
		    descr (iq:iq) = ' '
		END DO
	    END IF
	    CALL ST_LSTR ( unit, lng, ier )
	    IF ( lng .lt. 12 ) THEN
		DO iq = lng+1, 12
		    unit (iq:iq) = ' '
		END DO
	    END IF
	    WRITE (lunout, 1003) prm (1:8), descr (1:40),
     +				 iscl, iref, ibit, unit (1:12)
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!----------------------------------------------------------------------
! NEMO system team, System and Interface for oceanic RElocable Nesting
!----------------------------------------------------------------------
!
! DESCRIPTION:
!> @file
!> @brief
!> this program creates fine grid coordinate file.
!>
!> @details
!> @section sec1 method
!> variables from the input coordinates coarse/source grid file, are extracted
!> and interpolated to create a fine/taget grid coordinates file.

!> @note
!>    interpolation method could be different for each variable.
!>
!> \image html  header_coord_40.png
!>  \image latex header_coord_40.png
!>
!>
!> @section sec2 how to
!> USAGE: create_coord create_coord.nam [-v] [-h]

!>    - positional arguments:

!>       - create_coord.nam

!>          namelist of create_coord
!>          @note
!>             a template of the namelist could be created running (in templates directory):
!>             @code{.sh}
!>                python create_templates.py create_coord
!>             @endcode
!>
!>    - optional arguments:

!>       - -h, --help

!>          show this help message (and exit)

!>       - -v, --version

!>          show Siren's version   (and exit)
!>
!> @section sec_coord create_coord.nam
!>    create_coord.nam contains 6 sub-namelists:

!>       - **namlog** to set logger parameters
!>       - **namcfg** to set configuration file parameters
!>       - **namsrc** to set source/coarse grid parameters
!>       - **namvar** to set variable parameters
!>       - **namnst** to set sub domain and nesting paramters
!>       - **namout** to set output parameters
!>
!>    here after, each sub-namelist parameters is detailed.
!>    @note
!>       default values are specified between brackets
!>
!> @subsection sublog namlog
!>    the logger sub-namelist parameters are :
!>
!>    - **cn_logfile** [@a create_coord.log]

!>       logger filename
!>
!>    - **cn_verbosity** [@a warning]

!>       verbosity level, choose between :
!>          - trace
!>          - debug
!>          - info
!>          - warning
!>          - error
!>          - fatal
!>          - none
!>
!>    - **in_maxerror** [@a 5]

!>       maximum number of error allowed
!>
!> @subsection subcfg namcfg
!>    the configuration sub-namelist parameters are :
!>
!>    - **cn_varcfg** [@a ./cfg/variable.cfg]

!>       path to the variable configuration file.

!>       the variable configuration file defines standard name,
!>       default interpolation method, axis,...
!>       to be used for some known variables.

!>
!>    - **cn_dimcfg** [@a ./cfg/dimension.cfg]

!>       path to the dimension configuration file.

!>       the dimension configuration file defines dimensions allowed.

!>
!>    - **cn_dumcfg** [@a ./cfg/dummy.cfg]

!>       path to the useless (dummy) configuration file.

!>       the dummy configuration file defines useless
!>       dimension or variable. these dimension(s) or variable(s) will not be
!>       processed.

!>
!> @subsection subsrc namsrc
!>    the source/coarse grid sub-namelist parameters are :
!>
!>    - **cn_coord0** [@a ]

!>       path to the coordinate file
!>
!>    - **in_perio0** [@a ]

!>       NEMO periodicity index

!>       the NEMO periodicity could be choose between 0 to 6:
!>

!>
in_perio=0

!>
standard regional model

!>
in_perio=1

!>
east-west cyclic model

!>
in_perio=2

!>
model with symmetric boundary condition across the equator

!>
in_perio=3

!>
regional model with North fold boundary and T-point pivot

!>
in_perio=4

!>
global model with a T-point pivot.

!>          example: ORCA2, ORCA025, ORCA12

!>
in_perio=5

!>
regional model with North fold boundary and F-point pivot

!>
in_perio=6

!>
global model with a F-point pivot

!>          example: ORCA05

!>

!>

!>       @sa For more information see @ref md_src_docsrc_6_perio
!>       and Model Boundary Condition paragraph in the
!>       [NEMO documentation](https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/manual/pdf/NEMO_manual.pdf)
!>
!> @subsection subvar namvar
!>    the variable sub-namelist parameters are :
!>
!>    - **cn_varinfo** [@a ]

!>       list of variable and extra information about request(s) to be used

!>
!>       each elements of *cn_varinfo* is a string character (separated by ',').

!>       it is composed of the variable name follow by ':',
!>       then request(s) to be used on this variable.

!>       request could be:
!>          - int = interpolation method
!>          - ext = extrapolation method
!>
!>             requests must be separated by ';'.

!>             order of requests does not matter.

!>
!>       informations about available method could be find in @ref interp,
!>       @ref extrap and @ref filter modules.

!>       Example:
!>          - 'glamt: int=linear; ext=dist_weight', 'e1t: int=cubic/rhoi'
!>
!>       @note
!>          If you do not specify a method which is required,
!>          default one is apply.
!>
!> @subsection subnst namnst
!>    the nesting sub-namelist parameters are :
!>
!>    - **in_imin0** [@a ]

!>       i-direction lower left  point indice of coarse grid subdomain to be used
!>    - **in_imax0** [@a ]

!>       i-direction upper right point indice of coarse grid subdomain to be used
!>    - **in_jmin0** [@a ]

!>       j-direction lower left  point indice of coarse grid subdomain to be used
!>    - **in_jmax0** [@a ]

!>       j-direction upper right point indice of coarse grid subdomain to be used
!>
or

!>    - **rn_lonmin0** [@a ]

!>       lower left  longitude of coarse grid subdomain to be used
!>    - **rn_lonmax0** [@a ]

!>       upper right longitude of coarse grid subdomain to be used
!>    - **rn_latmin0** [@a ]

!>       lower left  latitude  of coarse grid subdomain to be used
!>    - **rn_latmax0** [@a ]

!>       upper right latitude  of coarse grid subdomain to be used
!>    @note you could define sub domain with
!>       - coarse/source grid indices
!>
or

!>       - coordinates.

!>    if coordinates are defined (-180 < lon < 360 and -90 < lat < 90),
!>    SIREN does not take into account indices.
!>
!>    - **in_rhoi**  [@a 1]

!>       refinement factor in i-direction
!>
!>    - **in_rhoj**  [@a 1]

!>       refinement factor in j-direction
!>
!>       \image html  grid_zoom_60.png
!>        \image latex grid_zoom_40.png
!>
!>
!> @subsection subout namout
!>    the output sub-namelist parameter is :
!>
!>    - **cn_fileout** [@a coord_fine.nc]

!>       output bathymetry filename
!>
!>

!> @author J.Paul
!>
!> @date November, 2013 - Initial Version
!> @date September, 2014
!> - add header for user
!> - compute offset considering grid point
!> - add global attributes in output file
!> @date September, 2015
!> - manage useless (dummy) variable, attributes, and dimension
!> @date September, 2016
!> - allow to use coordinate to define subdomain
!> @date October, 2016
!> - dimension to be used select from configuration file
!> @date January, 2019
!> - add url path to global attributes of output file(s)
!> @date February, 2019
!> - rename sub namelist namcrs to namsrc
!> - create and clean file structure to avoid memory leaks
!> @date Ocober, 2019
!> - add help and version optional arguments
!>
!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
!----------------------------------------------------------------------
PROGRAM create_coord

   USE global                          ! global variable
   USE kind                            ! F90 kind parameter
   USE logger                          ! log file manager
   USE fct                             ! basic useful function
   USE date                            ! date manager
   USE att                             ! attribute manager
   USE dim                             ! dimension manager
   USE var                             ! variable manager
   USE file                            ! file manager
   USE iom                             ! I/O manager
   USE grid                            ! grid manager
   USE extrap                          ! extrapolation manager
   USE interp                          ! interpolation manager
   USE filter                          ! filter manager
   USE mpp                             ! MPP manager
   USE dom                             ! domain manager
   USE iom_mpp                         ! MPP I/O manager
   USE iom_dom                         ! DOM I/O manager

   IMPLICIT NONE

   ! parameters
   CHARACTER(LEN=lc), PARAMETER  :: cp_myname = "create_coord"

   ! local variable
   CHARACTER(LEN=lc)                                    :: cl_arg
   CHARACTER(LEN=lc)                                    :: cl_namelist
   CHARACTER(LEN=lc)                                    :: cl_date
   CHARACTER(LEN=lc)                                    :: cl_url
   CHARACTER(LEN=lc)                                    :: cl_errormsg

   INTEGER(i4)                                          :: il_narg
   INTEGER(i4)                                          :: il_status
   INTEGER(i4)                                          :: il_fileid
   INTEGER(i4)                                          :: il_attid
   INTEGER(i4)                                          :: il_ind
   INTEGER(i4)                                          :: il_nvar
   INTEGER(i4)                                          :: il_ew
   INTEGER(i4)                                          :: il_imin0
   INTEGER(i4)                                          :: il_imax0
   INTEGER(i4)                                          :: il_jmin0
   INTEGER(i4)                                          :: il_jmax0

   INTEGER(i4)      , DIMENSION(ip_maxdim)              :: il_rho
   INTEGER(i4)      , DIMENSION(2)                      :: il_index
   INTEGER(i4)      , DIMENSION(2,2,ip_npoint)          :: il_offset

   LOGICAL                                              :: ll_exist

   TYPE(TATT)                                           :: tl_att

   TYPE(TDOM)                                           :: tl_dom

   TYPE(TVAR)       , DIMENSION(:)        , ALLOCATABLE :: tl_var

   TYPE(TDIM)       , DIMENSION(ip_maxdim)              :: tl_dim

   TYPE(TFILE)                                          :: tl_file

   TYPE(TMPP)                                           :: tl_coord0
   TYPE(TFILE)                                          :: tl_fileout

   ! loop indices
   INTEGER(i4) :: ji
   INTEGER(i4) :: jj

   ! namelist variable
   ! namlog
   CHARACTER(LEN=lc) :: cn_logfile  = 'create_coord.log'
   CHARACTER(LEN=lc) :: cn_verbosity= 'warning'
   INTEGER(i4)       :: in_maxerror = 5

   ! namcfg
   CHARACTER(LEN=lc) :: cn_varcfg   = './cfg/variable.cfg'
   CHARACTER(LEN=lc) :: cn_dimcfg   = './cfg/dimension.cfg'
   CHARACTER(LEN=lc) :: cn_dumcfg   = './cfg/dummy.cfg'

   ! namsrc
   CHARACTER(LEN=lc) :: cn_coord0   = ''
   INTEGER(i4)       :: in_perio0   = -1

   ! namvar
   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = ''

   !namnst
   REAL(sp)          :: rn_lonmin0  = -360.
   REAL(sp)          :: rn_lonmax0  = -360.
   REAL(sp)          :: rn_latmin0  = -360.
   REAL(sp)          :: rn_latmax0  = -360.
   INTEGER(i4)       :: in_imin0 = 0
   INTEGER(i4)       :: in_imax0 = 0
   INTEGER(i4)       :: in_jmin0 = 0
   INTEGER(i4)       :: in_jmax0 = 0
   INTEGER(i4)       :: in_rhoi  = 1
   INTEGER(i4)       :: in_rhoj  = 1

   !namout
   CHARACTER(LEN=lc) :: cn_fileout  = 'coord_fine.nc'
   !-------------------------------------------------------------------

   NAMELIST /namlog/ &  !  logger namelist
   &  cn_logfile,    &  !< logger file name
   &  cn_verbosity,  &  !< logger verbosity
   &  in_maxerror       !< logger maximum error

   NAMELIST /namcfg/ &  !< configuration namelist
   &  cn_varcfg,     &  !< variable configuration file
   &  cn_dimcfg,     &  !< dimension configuration file
   &  cn_dumcfg         !< dummy configuration file

   NAMELIST /namsrc/ &  !< source/coarse grid namelist
   &  cn_coord0 ,    &  !< coordinate file
   &  in_perio0         !< periodicity index

   NAMELIST /namvar/ &  !< variable namelist
   &  cn_varinfo        !< list of variable and extra information about
                        !< interpolation, extrapolation or filter method to be used.
                        !< (ex: 'votemper:linear,hann,dist_weight','vosaline:cubic' )

   NAMELIST /namnst/ &  !< nesting namelist
   &  rn_lonmin0,    &  !< lower left  coarse grid longitude
   &  rn_lonmax0,    &  !< upper right coarse grid longitude
   &  rn_latmin0,    &  !< lower left  coarse grid latitude
   &  rn_latmax0,    &  !< upper right coarse grid latitude
   &  in_imin0,      &  !< i-direction lower left  point indice
   &  in_imax0,      &  !< i-direction upper right point indice
   &  in_jmin0,      &  !< j-direction lower left  point indice
   &  in_jmax0,      &  !< j-direction upper right point indice
   &  in_rhoi,       &  !< refinement factor in i-direction
   &  in_rhoj           !< refinement factor in j-direction

   NAMELIST /namout/ &  !< output namelist
   &  cn_fileout        !< fine grid coordinate file
   !-------------------------------------------------------------------

   !
   ! Initialisation
   ! --------------
   !
   il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec

   ! Traitement des arguments fournis
   ! --------------------------------
   IF( il_narg /= 1 )THEN
      WRITE(cl_errormsg,*) ' ERROR : one argument is needed '
      CALL fct_help(cp_myname,cl_errormsg)
      CALL EXIT(1)
   ELSE

      CALL GET_COMMAND_ARGUMENT(1,cl_arg) !f03 intrinsec
      SELECT CASE (cl_arg)
         CASE ('-v', '--version')

            CALL fct_version(cp_myname)
            CALL EXIT(0)

         CASE ('-h', '--help')

            CALL fct_help(cp_myname)
            CALL EXIT(0)

         CASE DEFAULT

            cl_namelist=cl_arg

            ! read namelist
            INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
            IF( ll_exist )THEN

               il_fileid=fct_getunit()

               OPEN( il_fileid, FILE=TRIM(cl_namelist),  &
               &                FORM='FORMATTED',        &
               &                ACCESS='SEQUENTIAL',     &
               &                STATUS='OLD',            &
               &                ACTION='READ',           &
               &                IOSTAT=il_status)
               CALL fct_err(il_status)
               IF( il_status /= 0 )THEN
                  WRITE(cl_errormsg,*) " ERROR : error opening "//TRIM(cl_namelist)
                  CALL fct_help(cp_myname,cl_errormsg)
                  CALL EXIT(1)
               ENDIF

               READ( il_fileid, NML = namlog )

               ! define logger file
               CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
               CALL logger_header()

               READ( il_fileid, NML = namcfg )
               ! get variable extra information on configuration file
               CALL var_def_extra(TRIM(cn_varcfg))

               ! get dimension allowed
               CALL dim_def_extra(TRIM(cn_dimcfg))

               ! get dummy variable
               CALL var_get_dummy(TRIM(cn_dumcfg))
               ! get dummy dimension
               CALL dim_get_dummy(TRIM(cn_dumcfg))
               ! get dummy attribute
               CALL att_get_dummy(TRIM(cn_dumcfg))

               READ( il_fileid, NML = namsrc )
               READ( il_fileid, NML = namvar )
               ! add user change in extra information
               CALL var_chg_extra( cn_varinfo )

               READ( il_fileid, NML = namnst )
               READ( il_fileid, NML = namout )

               CLOSE( il_fileid, IOSTAT=il_status )
               CALL fct_err(il_status)
               IF( il_status /= 0 )THEN
                  CALL logger_error("CREATE COORD: closing "//TRIM(cl_namelist))
               ENDIF

            ELSE

               WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cl_namelist)
               CALL fct_help(cp_myname,cl_errormsg)
               CALL EXIT(1)

            ENDIF

      END SELECT
   ENDIF

   ! open files
   IF( cn_coord0 /= '' )THEN
      tl_file=file_init(TRIM(cn_coord0))
      tl_coord0=mpp_init( tl_file, id_perio=in_perio0)
      ! clean
      CALL file_clean(tl_file)
      CALL grid_get_info(tl_coord0)
   ELSE
      CALL logger_fatal("CREATE COORD: no coarse grid coordinate found. "//&
      &     "check namelist")
   ENDIF

   ! check
   ! check output file do not already exist
   INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist)
   IF( ll_exist )THEN
      CALL logger_fatal("CREATE COORD: output file "//TRIM(cn_fileout)//&
      &  " already exist.")
   ENDIF

   ! check nesting parameters
   il_index(:)=0
   IF( rn_lonmin0 >= -180. .AND. rn_lonmin0 <= 360 .AND. &
     & rn_latmin0 >= -90.  .AND. rn_latmin0 <= 90. )THEN

      il_index(:)=grid_get_closest(tl_coord0, &
         &                         REAL(rn_lonmin0,dp), REAL(rn_latmin0,dp), &
         &                         cd_pos='ll')
      il_imin0=il_index(1)
      il_jmin0=il_index(2)
   ELSE
      il_imin0=in_imin0
      il_jmin0=in_jmin0
   ENDIF

   il_index(:)=0
   IF( rn_lonmax0 >= -180. .AND. rn_lonmax0 <= 360 .AND. &
     & rn_latmax0 >= -90.  .AND. rn_latmax0 <= 90. )THEN

      il_index(:)=grid_get_closest(tl_coord0, &
         &                         REAL(rn_lonmax0,dp), REAL(rn_latmax0,dp), &
         &                         cd_pos='ur')
      il_imax0=il_index(1)
      il_jmax0=il_index(2)
   ELSE
      il_imax0=in_imax0
      il_jmax0=in_jmax0
   ENDIF

   ! forced indices for east west cyclic domain
   IF( rn_lonmin0 == rn_lonmax0 .AND. &
     & rn_lonmin0 /= -360. )THEN
      il_imin0=0
      il_imax0=0
   ENDIF

   IF( il_imin0 < 0 .OR. il_imax0 < 0 .OR. il_jmin0 < 0 .OR. il_jmax0 < 0)THEN
      CALL logger_fatal("CREATE COORD: invalid points indices."//&
      &  " check namelist "//TRIM(cl_namelist))
   ENDIF

   il_rho(:)=1
   IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN
      CALL logger_error("CREATE COORD: invalid refinement factor."//&
      &  " check namelist "//TRIM(cl_namelist))
   ELSE
      il_rho(jp_I)=in_rhoi
      il_rho(jp_J)=in_rhoj

      il_offset(:,:,:)=create_coord_get_offset(il_rho(:))
   ENDIF

   ! check domain validity
   CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0 )

   ! compute domain
   tl_dom=dom_init( tl_coord0,         &
   &                il_imin0, il_imax0,&
   &                il_jmin0, il_jmax0 )

   ! add extra band (if need be) to compute interpolation
   CALL dom_add_extra(tl_dom)

   ! open mpp files
   CALL iom_dom_open(tl_coord0, tl_dom)

   il_nvar=tl_coord0%t_proc(1)%i_nvar
   ALLOCATE( tl_var(il_nvar) )
   DO ji=1,il_nvar

      tl_var(ji)=iom_dom_read_var(tl_coord0, &
      &                          TRIM(tl_coord0%t_proc(1)%t_var(ji)%c_name),&
      &                          tl_dom)

      SELECT CASE(TRIM(tl_var(ji)%c_point))
         CASE('T')
            jj=jp_T
         CASE('U')
            jj=jp_U
         CASE('V')
            jj=jp_V
         CASE('F')
            jj=jp_F
      END SELECT

      ! interpolate variables
      CALL create_coord_interp( tl_var(ji), il_rho(:), &
      &                         il_offset(:,:,jj) )

      ! remove extraband added to domain
      CALL dom_del_extra( tl_var(ji), tl_dom, il_rho(:), .true. )

      ! filter
      CALL filter_fill_value(tl_var(ji))

   ENDDO

   ! clean
   CALL dom_clean_extra( tl_dom )

   ! close mpp files
   CALL iom_dom_close(tl_coord0)

   ! clean
   CALL mpp_clean(tl_coord0)

   ! create file
   tl_fileout=file_init(TRIM(cn_fileout))

   ! add dimension
   ! save biggest dimension
   tl_dim(:)=var_max_dim(tl_var(:))

   DO ji=1,ip_maxdim
      IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji))
   ENDDO

   ! add variables
   DO ji=il_nvar,1,-1
      CALL file_add_var(tl_fileout, tl_var(ji))
      CALL var_clean(tl_var(ji))
   ENDDO

   ! add some attribute
   tl_att=att_init("Created_by","SIREN create_coord")
   CALL file_add_att(tl_fileout, tl_att)

   !add source url
   cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
   tl_att=att_init("SIREN_url",cl_url)
   CALL file_add_att(tl_fileout, tl_att)

   ! add date of creation
   cl_date=date_print(date_now())
   tl_att=att_init("Creation_date",cl_date)
   CALL file_add_att(tl_fileout, tl_att)

   tl_att=att_init("src_file",TRIM(fct_basename(cn_coord0)))
   CALL file_add_att(tl_fileout, tl_att)

   tl_att=att_init("src_i_indices",(/tl_dom%i_imin,tl_dom%i_imax/))
   CALL file_add_att(tl_fileout, tl_att)
   tl_att=att_init("src_j_indices",(/tl_dom%i_jmin,tl_dom%i_jmax/))
   CALL file_add_att(tl_fileout, tl_att)
   IF( .NOT. ALL(il_rho(:)==1) )THEN
      tl_att=att_init("refinment_factor",(/il_rho(jp_I),il_rho(jp_J)/))
      CALL file_add_att(tl_fileout, tl_att)
   ENDIF

   ! add attribute periodicity
   il_attid=0
   IF( ASSOCIATED(tl_fileout%t_att) )THEN
      il_attid=att_get_index(tl_fileout%t_att(:),'periodicity')
   ENDIF
   IF( tl_dom%i_perio >= 0 .AND. il_attid == 0 )THEN
      tl_att=att_init('periodicity',tl_dom%i_perio)
      CALL file_add_att(tl_fileout,tl_att)
   ENDIF

   ! add attribute east west overlap
   il_attid=0
   IF( ASSOCIATED(tl_fileout%t_att) )THEN
      il_attid=att_get_index(tl_fileout%t_att(:),'ew_overlap')
   ENDIF
   IF( il_attid == 0 )THEN
      il_ind=var_get_index(tl_fileout%t_var(:),'longitude')
      IF( il_ind == 0 )THEN
         il_ind=var_get_index(tl_fileout%t_var(:),'longitude_T')
      ENDIF
      il_ew=grid_get_ew_overlap(tl_fileout%t_var(il_ind))
      IF( il_ew >= 0 )THEN
         tl_att=att_init('ew_overlap',il_ew)
         CALL file_add_att(tl_fileout,tl_att)
      ENDIF
   ENDIF

   ! create file
   CALL iom_create(tl_fileout)

   ! write file
   CALL iom_write_file(tl_fileout)

   ! close file
   CALL iom_close(tl_fileout)

   ! clean
   CALL att_clean(tl_att)
   CALL var_clean(tl_var(:))
   DEALLOCATE( tl_var)

   CALL file_clean(tl_fileout)
   CALL var_clean_extra()

   ! close log file
   CALL logger_footer()
   CALL logger_close()

CONTAINS
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   FUNCTION create_coord_get_offset(id_rho) &
         & RESULT (if_offset)
   !-------------------------------------------------------------------
   !> @brief
   !> This function compute offset over Arakawa grid points,
   !> given refinement factor.
   !>
   !> @author J.Paul
   !> @date August, 2014 - Initial Version
   !>
   !> @param[in] id_rho array of refinement factor
   !> @return array of offset
   !-------------------------------------------------------------------

      IMPLICIT NONE

      ! Argument
      INTEGER(i4), DIMENSION(:), INTENT(IN) :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2,ip_npoint) :: if_offset

      ! local variable
      ! loop indices
      !----------------------------------------------------------------

      ! case 'T'
      if_offset(jp_I,:,jp_T)=FLOOR(REAL(id_rho(jp_I)-1,dp)*0.5)
      if_offset(jp_J,:,jp_T)=FLOOR(REAL(id_rho(jp_J)-1,dp)*0.5)
      ! case 'U'
      if_offset(jp_I,1,jp_U)=0
      if_offset(jp_I,2,jp_U)=id_rho(jp_I)-1
      if_offset(jp_J,:,jp_U)=FLOOR(REAL(id_rho(jp_J)-1,dp)*0.5)
      ! case 'V'
      if_offset(jp_I,:,jp_V)=FLOOR(REAL(id_rho(jp_I)-1,dp)*0.5)
      if_offset(jp_J,1,jp_V)=0
      if_offset(jp_J,2,jp_V)=id_rho(jp_J)-1
      ! case 'F'
      if_offset(jp_I,1,jp_F)=0
      if_offset(jp_I,2,jp_F)=id_rho(jp_I)-1
      if_offset(jp_J,1,jp_F)=0
      if_offset(jp_J,2,jp_F)=id_rho(jp_J)-1

   END FUNCTION create_coord_get_offset
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   SUBROUTINE create_coord_interp(td_var, id_rho, id_offset, &
         &                        id_iext, id_jext)
   !-------------------------------------------------------------------
   !> @brief
   !> This subroutine interpolate variable, given refinment factor.
   !>
   !> @details
   !>  Optionaly, you could specify number of points
   !>    to be extrapolated in i- and j-direction.

   !>  variable mask is first computed (using _FillValue) and interpolated.

   !>  variable is then extrapolated, and interpolated.

   !>  Finally interpolated mask is applied on refined variable.
   !>
   !> @author J.Paul
   !> @date November, 2013 - Initial Version
   !>
   !> @param[inout] td_var variable strcuture
   !> @param[in] id_rho    array of refinement factor
   !> @param[in] id_offset offset between fine grid and coarse grid
   !> @param[in] id_iext   number of points to be extrapolated in i-direction
   !> @param[in] id_jext   number of points to be extrapolated in j-direction
   !>
   !> @todo check if mask is really needed
   !-------------------------------------------------------------------

      IMPLICIT NONE

      ! Argument
      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
      INTEGER(i4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
      INTEGER(i4), DIMENSION(:,:), INTENT(IN)    :: id_offset
      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_iext
      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_jext

      ! local variable
      TYPE(TVAR)  :: tl_mask

      INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask

      INTEGER(i4) :: il_iext
      INTEGER(i4) :: il_jext

      ! loop indices
      !----------------------------------------------------------------

      IF( ANY(SHAPE(id_offset(:,:)) /= 2) )THEN
         CALL logger_error("CREATE COORD INTERP: invalid dimension of "//&
         &                 "offset array")
      ENDIF

      !WARNING: two extrabands are required for cubic interpolation
      il_iext=2
      IF( PRESENT(id_iext) ) il_iext=id_iext

      il_jext=2
      IF( PRESENT(id_jext) ) il_jext=id_jext

      IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
         CALL logger_warn("CREATE COORD INTERP: at least extrapolation "//&
         &  "on two points are required with cubic interpolation ")
         il_iext=2
      ENDIF

      IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN
         CALL logger_warn("CREATE COORD INTERP: at least extrapolation "//&
         &  "on two points are required with cubic interpolation ")
         il_jext=2
      ENDIF

      IF( ANY(id_rho(:)>1) )THEN
         ! work on mask
         ! create mask
         ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
         &                td_var%t_dim(2)%i_len, &
         &                td_var%t_dim(3)%i_len, &
         &                td_var%t_dim(4)%i_len) )

         bl_mask(:,:,:,:)=1
         WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0

         SELECT CASE(TRIM(td_var%c_point))
         CASE DEFAULT ! 'T'
            tl_mask=var_init('tmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
            &                id_ew=td_var%i_ew )
         CASE('U')
            tl_mask=var_init('umask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
            &                id_ew=td_var%i_ew )
         CASE('V')
            tl_mask=var_init('vmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
            &                id_ew=td_var%i_ew )
         CASE('F')
            tl_mask=var_init('fmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),&
            &                id_ew=td_var%i_ew )
         END SELECT

         DEALLOCATE(bl_mask)

         ! interpolate mask
         CALL interp_fill_value( tl_mask, id_rho(:), &
         &                       id_offset=id_offset(:,:) )

         ! work on variable
         ! add extraband
         CALL extrap_add_extrabands(td_var, il_iext, il_jext)

         ! extrapolate variable
         CALL extrap_fill_value( td_var )

         ! interpolate variable
         CALL interp_fill_value( td_var, id_rho(:), &
         &                       id_offset=id_offset(:,:))

         ! remove extraband
         CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J))

         ! keep original mask
         WHERE( tl_mask%d_value(:,:,:,:) == 0 )
            td_var%d_value(:,:,:,:)=td_var%d_fill
         END WHERE
      ENDIF

      ! clean variable structure
      CALL var_clean(tl_mask)

   END SUBROUTINE create_coord_interp
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
END PROGRAM create_coord


text:
 SUBROUTINE dsurf1(nsurf,ndime,ndofn,npoin,iwrit,coord,loadv,heads)
!     Apply load distributed over surface
 USE lispa0, ONLY : lures
 USE loa_db
 IMPLICIT NONE
 INTEGER (kind=4), INTENT (IN) :: nsurf,ndime,ndofn,npoin,iwrit
 REAL (kind=8), INTENT(IN) :: coord(ndime,npoin)
 REAL (kind=8), INTENT(IN OUT) :: loadv(ndofn,npoin)
 TYPE (srf_nod),POINTER :: heads

 INTEGER (kind=4), PARAMETER :: nnmax = 6
 INTEGER (kind=4) :: i,n,l,nn,isurf,nvert,&
                     noprs(nnmax),g,ll(3)
 REAL (kind=8) :: x(ndime,nnmax),direc(ndime),&
                  dx(ndime),dy(ndime),dz(ndime),&
                  rload(ndime,nnmax),px(ndime),p(nnmax),&
                  auxv,area,fc
 REAL (kind=8) weigh(3,3),posgp(9),shape(21),deriv(21),gldir(ndime)
 INTEGER (kind=4) :: ii(4)= (/6,4,5,6/), kk(4)= (/4,5,6,5/)
 TYPE (srf_nod),POINTER :: srf

 CHARACTER :: systm
 INTEGER (kind=4) chnode

 IF(iwrit == 1) WRITE(lures,"(//6x,'No. of Loaded Surfaces =',i5,&
&                            /)",ERR=9999) nsurf

 !*** distributed surface loads section

 IF(iwrit == 1) WRITE(lures,"(/6x,'list of loaded Faces  and ',&
&  'Applied loads'/)",ERR=9999)

 srf => heads
 DO isurf = 1,nsurf

   systm = srf%systm3
   nn = srf%nn3
   noprs(1:nn) = srf%lnod3(1:nn)
   nvert = nn
   IF( nn > 4) nvert = nn/2
   DO i=1,nn
     noprs(i) = chnode(noprs(i))
   END DO
   direc = 0d0
   direc(1:ndime) = srf%direc3(1:ndime)
   CALL vecuni(ndime,direc,area)
   IF (area == 0d0) direc(3) = 1d+0
   p(1:nn) = srf%press3(1:nn)

   !       prints actual DATA
   IF(iwrit == 1) THEN
     WRITE(lures,"(' FACE',i3,' Sys= ',a1,' NN=',i2,' DIR'&
&         ,3f6.3,' Nodes',6i6)",ERR=9999) isurf,systm,nn,direc,noprs(1:NN)
     WRITE(lures,"(' Nodal p ',6g12.3)",ERR=9999) p(1:nn)
   END IF

   rload = 0d0
   x(1:ndime,1:nn) = coord(1:ndime,noprs(1:nn))

   SELECT CASE (nn)
   CASE (3)
     dx = x(1:ndime,2) - x(1:ndime,1)
     dy = x(1:ndime,3) - x(1:ndime,1)
     CALL vecpro(dx,dy,dz)
     CALL vecuni(ndime,dz,area)
     area = area/24d0    !area/12

     !         direction of the load

     IF ( systm == 'L') THEN
       !         SELECT local x=t1 in the global xy plane
       dx = (/ -dz(2), dz(1) , 0d0 /)
       !         of course local y = t2 = t3 x t1
       dy=(/-dz(1)*dz(3),-dz(2)*dz(3),(dz(1)*dz(1)+dz(2)*dz(2))/)
       IF(ABS(dy(3)) < 1.0d-5) THEN
         dx = (/  1d0, 0d0, 0d0 /)
         dy = (/  0d0, 1d0, 0d0 /)
       ELSE
         !           normalizes t1 & t2
         CALL vecuni(ndime,dx,auxv)
         CALL vecuni(ndime,dy,auxv)
       END IF
       !  transform direc from local to global
       gldir = direc(1)*dx+direc(2)*dy+direc(3)*dz

     ELSE

       gldir = direc

     END IF
     px(1) = 2*p(1)+  p(2)+  p(3)
     px(2) =   p(1)+2*p(2)+  p(3)
     px(3) =   p(1)+  p(2)+2*p(3)

     DO g=1,nvert
       auxv = area*px(g)
       rload(1:ndime,g) = rload(1:ndime,g) + gldir*auxv
     END DO

   CASE (4)
     CALL gaussq(2 ,posgp ,weigh )
     DO g = 1,4
       CALL shape3(deriv,shape,posgp(MOD(g-1,2)+1),posgp((g+1)/2),     4)
       dx(1:ndime) = MATMUL(x(1:ndime,1:4),deriv(1:4))
       dy(1:ndime) = MATMUL(x(1:ndime,1:4),deriv(5:8))
       CALL vecpro(dx,dy,dz)
       CALL vecuni(ndime,dz,area)

       ! direction of the load
       IF (systm == 'L') THEN
         ! transform direc from local to global
         !         SELECT local x=t1 in the global xy plane
         dx = (/ -dz(2), dz(1) , 0d0 /)
         !         of course local y = t2 = t3 x t1
         dy = (/-dz(1)*dz(3), -dz(2)*dz(3),(dz(1)*dz(1)+dz(2)*dz(2))/)
         IF(ABS(dy(3)) < 1.0d-5) THEN
           dx = (/  1d0, 0d0, 0d0 /)
           dy = (/  0d0, 1d0, 0d0 /)
         ELSE
           !           normalizes t1 & t2
           CALL vecuni(ndime,dx,auxv)
           CALL vecuni(ndime,dy,auxv)
         END IF
         gldir = direc(1)*dx+direc(2)*dy+direc(3)*dz

       ELSE

         gldir = direc

       END IF

       auxv = DOT_PRODUCT(shape(1:4),p(1:4))*area
       DO n = 1,nn
         fc = auxv*shape(n)
         rload(1:ndime,n) = rload(1:ndime,n) + fc*gldir
       END DO

     END DO
   CASE (6)
     ! for each sub-triangle
     DO n=1,4
       ll(1) = n
       ll(2) = ii(n)
       ll(3) = kk(n)
       dx = x(1:ndime,kk(n)) - x(1:ndime,n)
       dy = x(1:ndime,ii(n)) - x(1:ndime,n)
       CALL vecpro(dx,dy,dz)
       CALL vecuni(ndime,dz,area)
       area = area/24d0

       !           direction of the load

       IF ( systm == 'L') THEN
         !         SELECT local x=t1 in the global xy plane
         dx = (/ -dz(2), dz(1) , 0d0 /)
         !         of course local y = t2 = t3 x t1
         dy = (/-dz(1)*dz(3), -dz(2)*dz(3),(dz(1)*dz(1)+dz(2)*dz(2))/)
         IF(ABS(dy(3)) < 1.0d-5) THEN
           dx = (/  1d0, 0d0, 0d0 /)
           dy = (/  0d0, 1d0, 0d0 /)
         ELSE
           !           normalizes t1 & t2
           CALL vecuni(ndime,dx,auxv)
           CALL vecuni(ndime,dy,auxv)
         END IF
         !  transform direc from local to global
         gldir = direc(1)*dx+direc(2)*dy+direc(3)*dz

       ELSE

         gldir = direc

       END IF
       px(1) = 2*p(ll(1))+  p(ll(2))+  p(ll(3))
       px(2) =   p(ll(1))+2*p(ll(2))+  p(ll(3))
       px(3) =   p(ll(1))+  p(ll(2))+2*p(ll(3))

       DO g=1,nvert
         auxv = area*px(g)
         rload(1:ndime,ll(g))= rload(1:ndime,ll(g)) + gldir*auxv
       END DO
     END DO
   CASE DEFAULT
     CALL runend('DSURF1:ERROR NN can be only 3 4 6 ')
   END SELECT

   DO n=1,nn
     l=noprs(n)
     loadv(1:ndime,l) = loadv(1:ndime,l) + rload(1:ndime,n)
   END DO
   srf => srf%next
 END DO

 RETURN
 9999 CALL runen2('')
 END SUBROUTINE dsurf1

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


! { dg-do compile }
! { dg-options "-fcoarray=lib" }
!
! Support for stat= in caf reference
!
program whitespace
  implicit none

  integer :: me[*],tmp,stat,i

  me = this_image()
  stat = 0
  i = 1

  sync all(stat = stat)

  if(stat /= 0) write(*,*) 'failure during sync'

  stat = 0

  if(me == 1) then
     tmp = me[num_images(),stat = stat]
     if(stat /= 0) write(*,*) 'failure in img:',me
  else if(me == 2) then
     tmp = me[i,stat=stat]
     if(stat /= 0) write(*,*) 'failure in img:',me
  endif

end program whitespace


text:
!***********************************************************************
! Determines the form of precipitation and distributes precipitation
! and temperatures to each HRU using a multiple linear regression of
! measured data from a group of measurement stations or from atmospheric
! model simulation
!
!                      Converted by Steve Markstrom
!                      Wed Feb 10 15:16:04 MST 1999
!              revised Wed Mar 17 15:48:52 MST 1999
!              revised Mon Aug 30 16:47:07 MDT 1999
!              revised Mon Aug 30 16:47:07 MDT 1999
!              revised Wed Mar  8 09:06:18 MST 2000
!              revised Thu Feb  3 10:00:00 MST 2005
! temp_nsta - number of temperature stations used
! temp_nuse(temp_nsta) - indicies of temperature stations used
! rain_nsta - number of precipitation stations used
! rain_nuse (rain_nsta) - indicies of rain stations used
!***********************************************************************
      MODULE PRMS_XYZ_DIST
      IMPLICIT NONE
!   Local Variables
      INTEGER, PARAMETER :: MAXLAPSE = 3
      CHARACTER(LEN=8), SAVE :: MODNAME
      INTEGER, SAVE :: Nlapse, Temp_nsta, Rain_nsta
      INTEGER, SAVE, ALLOCATABLE :: Rain_nuse(:), Temp_nuse(:)
      DOUBLE PRECISION, SAVE :: Basin_centroid_x, Basin_centroid_y
      DOUBLE PRECISION, SAVE :: Meantmax(12), Meantmin(12)
      DOUBLE PRECISION, SAVE :: Temp_meanx(12), Temp_meany(12)
      DOUBLE PRECISION, SAVE :: Rain_meanx(12), Rain_meany(12)
      DOUBLE PRECISION, SAVE :: Temp_meanz(12), Rain_meanz(12)
      DOUBLE PRECISION, SAVE :: Meanppt(12)
      REAL, SAVE, ALLOCATABLE :: Precip_xyz(:)
      REAL, SAVE, ALLOCATABLE :: Temp_STAelev(:)
! transformed versions of these values
      REAL, SAVE, ALLOCATABLE :: MRUelev(:), Pstaelev(:)
      REAL, SAVE, ALLOCATABLE :: Pstax(:), Pstay(:)
      REAL, SAVE :: Solradelev
!   Declared Variables
      INTEGER, SAVE :: Is_rain_day
      REAL, SAVE, ALLOCATABLE :: Tmax_rain_sta(:), Tmin_rain_sta(:)
!   Declared Parameters
      INTEGER, SAVE :: Conv_flag
      REAL, SAVE, ALLOCATABLE :: Max_lapse(:, :), Min_lapse(:, :)
      REAL, SAVE, ALLOCATABLE :: Ppt_lapse(:, :)
      REAL, SAVE :: Solrad_elev
      REAL, SAVE, ALLOCATABLE :: Adjust_snow(:, :), Adjust_rain(:, :)
      REAL, SAVE :: Tmax_add, Tmax_div, Tmin_add, Tmin_div, Ppt_add
      REAL, SAVE :: X_add, X_div, Y_add, Y_div, Z_add, Z_div, Ppt_div
      INTEGER, SAVE, ALLOCATABLE :: Tsta_nuse(:)
      INTEGER, SAVE, ALLOCATABLE :: Psta_nuse(:), Psta_freq_nuse(:)
      REAL, SAVE, ALLOCATABLE :: MRUx(:), MRUy(:)
!     Temp_STAx = tsta_x, Temp_STAy = tsta_y
      REAL, SAVE, ALLOCATABLE :: Temp_STAx(:), Temp_STAy(:)
!     Rain_STAx = psta_x, Rain_STAy = psta_y
      REAL, SAVE, ALLOCATABLE :: Rain_STAx(:), Rain_STAy(:)
!     TmaxMTH = tsta_month_max, TminMTH = tsta_month_min
!     PptMTH = psta_month_ppt
      REAL, SAVE, ALLOCATABLE :: TmaxMTH(:, :), TminMTH(:, :)
      REAL, SAVE, ALLOCATABLE :: PptMTH(:, :)
      !REAL, SAVE, ALLOCATABLE :: Tmax_allsnow_sta(:, :)
      !REAL, SAVE, ALLOCATABLE :: Tmax_allrain_sta(:, :)
      REAL, SAVE :: Tmax_allsnow_dist, Tmax_allrain_dist(12)
      END MODULE PRMS_XYZ_DIST

!***********************************************************************
!     Main xyz_dist routine
!***********************************************************************
      INTEGER FUNCTION xyz_dist()
      USE PRMS_MODULE, ONLY: Process, Save_vars_to_file,
     &    Init_vars_from_file
      IMPLICIT NONE
! Functions
      INTEGER, EXTERNAL :: xyzdecl, xyzinit, xyzrun, xyzsetdims
      EXTERNAL :: xyz_restart
!***********************************************************************
      xyz_dist = 0

      IF ( Process(:3)=='run' ) THEN
        xyz_dist = xyzrun()
      ELSEIF ( Process(:7)=='setdims' ) THEN
        xyz_dist = xyzsetdims()
      ELSEIF ( Process(:4)=='decl' ) THEN
        xyz_dist = xyzdecl()
      ELSEIF ( Process(:4)=='init' ) THEN
        IF ( Init_vars_from_file==1 )  CALL xyz_restart(1)
        xyz_dist = xyzinit()
      ELSEIF ( Process(:5)=='clean' ) THEN
        IF ( Save_vars_to_file==1 ) CALL xyz_restart(0)
      ENDIF

      END FUNCTION xyz_dist

!***********************************************************************
!     xyzsetdims - declares xyz_dist specific dimensions
!***********************************************************************
      INTEGER FUNCTION xyzsetdims()
      USE PRMS_XYZ_DIST, ONLY: Nlapse
      IMPLICIT NONE
! Functions
      INTEGER, EXTERNAL :: declfix
      EXTERNAL read_error
!***********************************************************************
      xyzsetdims = 0

      IF ( declfix('nlapse', 3, 3,
     +     'Number of lapse rates in X, Y, and Z directions')
     +     /=0 ) CALL read_error(7, 'nlapse')
      Nlapse = 3

      END FUNCTION xyzsetdims
!***********************************************************************

!***********************************************************************
!     xyzdecl - set up parameters for temperature computations
!   Declared Parameters
!     hru_x, hru_y, max_lapse, min_lapse, ppt_lapse, tsta_elev
!     tmax_adj, tmin_adj, tsta_x, tsta_y, psta_elev, psta_x, psta_y
!     tsta_nuse, psta_nuse, psta_freq_nuse, adjust_snow, adjust_rain
!     tsta_month_max, tsta_month_min, psta_month_ppt, rain_code
!     x_add, x_div, y_add, y_div, z_add, z_div, solrad_elev
!     tmin_add, tmin_div, tmax_add, tmax_div, ppt_add, ppt_div
!     tmax_allrain, tmax_allsnow, adjmix_rain, conv_flag
!     hru_elev, hru_area
!***********************************************************************
      INTEGER FUNCTION xyzdecl()
      USE PRMS_XYZ_DIST
      USE PRMS_MODULE, ONLY: Nhru, Ntemp, Nrain
      IMPLICIT NONE
! Functions
      INTEGER, EXTERNAL :: declparam, declvar
      EXTERNAL read_error, print_module
! Local Variables
      CHARACTER(LEN=80), SAVE :: Version_xyz_dist
!***********************************************************************
      xyzdecl = 0

      Version_xyz_dist =
     +'xyz_dist.f 2016-07-22 17:22:00Z'
      CALL print_module(Version_xyz_dist,
     +                  'Temp & Precip Distribution  ', 77)
      MODNAME = 'xyz_dist'

      IF ( declvar(MODNAME, 'is_rain_day', 'one', 1, 'integer',
     +     'Flag to indicate if it is raining anywhere in the basin',
     +     'none',
     +     Is_rain_day)/=0 ) CALL read_error(3, 'is_rain_day')

      ALLOCATE ( Tmax_rain_sta(Nrain) )
      IF ( declvar(MODNAME, 'tmax_rain_sta', 'nrain', Nrain, 'real',
     +     'Maximum temperature distributed to the precipitation'//
     +     ' measurement stations',
     +     'degrees Fahrenheit',
     +     Tmax_rain_sta)/=0 ) CALL read_error(3, 'tmax_rain_sta')

      ALLOCATE ( Tmin_rain_sta(Nrain) )
      IF ( declvar(MODNAME, 'tmin_rain_sta', 'nrain', Nrain, 'real',
     +     'Minimum temperature distributed to the precipitation'//
     +     ' measurement stations',
     +     'degrees Fahrenheit',
     +     Tmin_rain_sta)/=0 ) CALL read_error(3, 'tmin_rain_sta')

! declare parameters
      ALLOCATE ( Adjust_snow(Nrain,12) )
      IF ( declparam(MODNAME, 'adjust_snow', 'nrain,nmonths', 'real',
     +     '-0.4', '-0.5', '0.5',
     +     'Monthly (January to December) snow downscaling adjustment'//
     +     ' factor for each precipitation measurement station',
     +     'Monthly (January to December) snow downscaling adjustment'//
     +     ' factor for each precipitation measurement station',
     +     'decimal fraction')/=0 ) CALL read_error(1, 'adjust_snow')

      ALLOCATE ( Adjust_rain(Nrain,12) )
      IF ( declparam(MODNAME, 'adjust_rain', 'nrain,nmonths', 'real',
     +     '-0.4', '-0.5', '0.5',
     +     'Monthly (January to December) rain downscaling adjustment'//
     +     ' factor for each precipitation measurement station',
     +     'Monthly (January to December) rain downscaling adjustment'//
     +     ' factor for each precipitation measurement station',
     +     'decimal fraction')/=0 ) CALL read_error(1, 'adjust_rain')

      ALLOCATE ( MRUelev(Nhru) )
      ALLOCATE ( MRUx(Nhru) )
      IF ( declparam(MODNAME, 'hru_x', 'nhru', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'X for each HRU (albers)',
     +     'Longitude (X) for each HRU in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'hru_x')

      ALLOCATE ( MRUy(Nhru) )
      IF ( declparam(MODNAME, 'hru_y', 'nhru', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'Y for each HRU (albers)',
     +     'Latitude (Y) for each HRU in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'hru_y')

      ALLOCATE ( Max_lapse(MAXLAPSE,12) )
      IF ( declparam(MODNAME, 'max_lapse', 'nlapse,nmonths', 'real',
     +     '0.0', '-100.0', '100.0',
     +     'Monthly maximum temperature lapse rate for each direction',
     +     'Monthly (January to December) maximum air temperature'//
     +     ' lapse rate for each direction (X, Y, and Z)',
     +     'none')/=0 ) CALL read_error(1, 'max_lapse')

      ALLOCATE ( Min_lapse(MAXLAPSE,12) )
      IF ( declparam(MODNAME, 'min_lapse', 'nlapse,nmonths', 'real',
     +     '0.0', '-100.0', '100.0',
     +     'Monthly minimum temperature lapse rate for each direction',
     +     'Monthly (January to December) minimum air temperature'//
     +     ' lapse rate for each direction (X, Y, and Z)',
     +     'none')/=0 ) CALL read_error(1, 'min_lapse')

      ALLOCATE ( Ppt_lapse(MAXLAPSE,12) )
      IF ( declparam(MODNAME, 'ppt_lapse', 'nlapse,nmonths', 'real',
     +     '0.0', '-10.0', '10.0',
     +     'Precipitation lapse rate',
     +     'Monthly (January to December) precipitation lapse rate'//
     +     ' for each direction (X, Y, and Z)',
     +     'none')/=0 ) CALL read_error(1, 'ppt_lapse')

      ALLOCATE ( Temp_STAx(Ntemp), Temp_STAelev(Ntemp) )
      IF ( declparam(MODNAME, 'tsta_x', 'ntemp', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'X for each temperature station (albers)',
     +     'Longitude (X) for each air-temperature-measurement'//
     +     ' station in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'tsta_x')

      ALLOCATE ( Temp_STAy(Ntemp) )
      IF ( declparam(MODNAME, 'tsta_y', 'ntemp', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'Y for each temperature station (albers)',
     +     'Latitude (Y) for each air-temperature-measurement'//
     +     ' station in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'tsta_y')

      ALLOCATE ( Rain_STAx(Nrain), Pstaelev(Nrain) )
      IF ( declparam(MODNAME, 'psta_x', 'nrain', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'X for each precipitation station (albers)',
     +     'Longitude (X) for each precipitation measurement station'//
     +     ' in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'psta_x')

      ALLOCATE ( Rain_STAy(Nrain) )
      IF ( declparam(MODNAME, 'psta_y', 'nrain', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'Y for each precipitation station (albers)',
     +     'Latitude (Y) for each precipitation measurement station'//
     +     ' in albers projection',
     +     'meters')/=0 ) CALL read_error(1, 'psta_y')

      ALLOCATE ( Tsta_nuse(Ntemp), Temp_nuse(Ntemp) )
      IF ( declparam(MODNAME, 'tsta_nuse', 'ntemp', 'integer',
     +     '1', '0', '1',
     +     '0 = station not used; 1 = station used',
     +     'The subset of temperature measurement stations used in'//
     +     ' the distribution regression (0=station not used;'//
     +     ' 1=station used)',
     +     'none')/=0 ) CALL read_error(1, 'tsta_nuse')

      IF ( declparam(MODNAME, 'solrad_elev', 'one', 'real',
     +     '0.0', '-300.0', '30000.0',
     +     'Elevation of the solrad station used for DD curves',
     +     'Elevation of the solar radiation station used for'//
     +     ' degree-day curves to distribute temperature',
     +     'meters')/=0 ) CALL read_error(1, 'solrad_elev')

      ALLOCATE ( Psta_nuse(Nrain), Rain_nuse(Nrain) )
      IF ( declparam(MODNAME, 'psta_nuse', 'nrain', 'integer',
     +     '1', '0', '1',
     +     'The subset of precipitation stations used in the'//
     +     ' distribution regression (0=station not used;'//
     +     ' 1=station used)',
     +     'The subset of precipitation measurement stations used in'//
     +     ' the distribution regression (0=station not used;'//
     +     ' 1=station used)',
     +     'none')/=0 ) CALL read_error(1, 'psta_nuse')

      ALLOCATE ( Psta_freq_nuse(Nrain) )
      IF ( declparam(MODNAME, 'psta_freq_nuse', 'nrain', 'integer',
     +     '1', '0', '1',
     +     'The subset of precipitation stations used to determine if'//
     +     ' there is distribution in the basin (0=station not used;'//
     +     ' 1=station used)',
     +     'The subset of precipitation measurement stations used to'//
     +     ' determine if there is precipitation in the basin'//
     +     ' (0=station not used; 1=station used)',
     +     'none')/=0 ) CALL read_error(1, 'psta_freq_nuse')

      ALLOCATE ( TmaxMTH(Ntemp, 12) )
      IF ( declparam(MODNAME, 'tsta_month_max', 'ntemp,nmonths',
     +     'real',
     +     '0.0', '-100.0', '100.0',
     +     'Average monthly (January to December) maximum air'//
     +     ' temperature at each station',
     +     'Average monthly (January to December) maximum air'//
     +     ' temperature at each air-temperature-measurement station',
     +     'temp_units')/=0 ) CALL read_error(1, 'tsta_month_max')

      ALLOCATE ( TminMTH(Ntemp, 12) )
      IF ( declparam(MODNAME, 'tsta_month_min', 'ntemp,nmonths',
     +     'real',
     +     '0.0', '-100.0', '100.0',
     +     'Average monthly (January to December) minimum air'//
     +     ' temperature at each station',
     +     'Average monthly (January to December) minimum air'//
     +     ' temperature at each air-temperature-measurement station',
     +     'temp_units')/=0 ) CALL read_error(1, 'tsta_month_min')

      ALLOCATE ( PptMTH(Nrain, 12) )
      IF ( declparam(MODNAME, 'psta_month_ppt', 'nrain,nmonths',
     +     'real',
     +     '0.0', '0.0', '20.0',
     +     'Average monthly precipitation at each station',
     +     'Average monthly (January to December) precipitation at'//
     +     ' each precipitation measurement station',
     +     'precip_units')/=0 ) CALL read_error(1, 'psta_month_ppt')

      IF ( declparam(MODNAME, 'x_add', 'one', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'X additive term for climate station transform',
     +     'Mean value for the climate station transformation'//
     +     ' equation for the longitude (X) coordinate',
     +     'meters')/=0 ) CALL read_error(1, 'x_add')

      IF ( declparam(MODNAME, 'x_div', 'one', 'real',
     +     '1.0', '-1.0E7', '1.0E7',
     +     'X divisor term for climate station transform',
     +     'Standard deviation for the climate station transformation'//
     +     ' equation for the longitude (X) coordinate (not 0.0)',
     +     'meters')/=0 ) CALL read_error(1, 'x_div')

      IF ( declparam(MODNAME, 'y_add', 'one', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'Y additive term for climate station transform',
     +     'Mean value for the climate station transformation'//
     +     ' equation for the latitude (Y) coordinate',
     +     'meters')/=0 ) CALL read_error(1, 'y_add')

      IF ( declparam(MODNAME, 'y_div', 'one', 'real',
     +     '1.0', '-1.0E7', '1.0E7',
     +     'Y divisor term for climate station transform',
     +     'Standard deviation for the climate station transformation'//
     +     ' equation for the latitude (Y) coordinate (not 0.0)',
     +     'meters')/=0 ) CALL read_error(1, 'y_div')

      IF ( declparam(MODNAME, 'z_add', 'one', 'real',
     +     '0.0', '-1.0E7', '1.0E7',
     +     'Z additive term for climate station transform',
     +     'Mean value for the climate station transformation'//
     +     ' equation for the elevation (Z) coordinate',
     +     'meters')/=0 ) CALL read_error(1, 'z_add')

      IF ( declparam(MODNAME, 'z_div', 'one', 'real',
     +     '1.0', '-1.0E7', '1.0E7',
     +     'Z divisor term for climate station transform',
     +     'Standard deviation for the climate station transformation'//
     +     ' equation for the elevation (Z) coordinate (not 0.0)',
     +     'meters')/=0 ) CALL read_error(1, 'z_div')

      IF ( declparam(MODNAME, 'tmax_add', 'one', 'real',
     +     '0.0', '-100.0', '100.0',
     +     'Maximum temperature additive term for climate station'//
     +     ' transform',
     +     'Mean value for the air-temperature-measurement station'//
     +     ' transformation equation for maximum temperature',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmax_add')

      IF ( declparam(MODNAME, 'tmax_div', 'one', 'real',
     +     '1.0', '-100.0', '100.0',
     +     'Maximum temperature divisor term for climate station'//
     +     ' transform',
     +     'Standard deviation for the air-temperature-measurement'//
     +     ' station transformation'//
     +     ' equation for maximum temperature (not 0.0)',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmax_div')

      IF ( declparam(MODNAME, 'tmin_add', 'one', 'real',
     +     '0.0', '-100.0', '100.0',
     +     'Minimum temperature additive term for the air-temperature'//
     +     '-measurement station transformation transform',
     +     'Mean value for the air-temperature-measurement station'//
     +     ' transformation equation for minimum temperature',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmin_add')

      IF ( declparam(MODNAME, 'tmin_div', 'one', 'real',
     +     '1.0', '-100.0', '100.0',
     +     'Minimum temperature divisor term for the air-temperature'//
     +     '-measurement station transformation transform',
     +     'Standard deviation for the air-temperature-measurement'//
     +     ' station transformation'//
     +     ' equation for minimum temperature (not 0.0)',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmin_div')

      IF ( declparam(MODNAME, 'ppt_add', 'one', 'real',
     +     '0.0', '-10.0', '10.0',
     +     'Precipitation additive term for precipitation measurement'//
     +     ' station transform',
     +     'Mean value for the precipitation measurement station'//
     +     ' transformation equation',
     +     'precip_units')/=0 ) CALL read_error(1, 'ppt_add')

      IF ( declparam(MODNAME, 'ppt_div', 'one', 'real',
     +     '1.0', '-10.0', '10.0',
     +     'Precipitation divisor term for precipitation measurement'//
     +     ' station transform',
     +     'Standard deviation for the precipitation measurement'//
     +     ' station transformation equation (not 0.0)',
     +     'precip_units')/=0 ) CALL read_error(1, 'ppt_div')

      IF ( declparam(MODNAME, 'conv_flag', 'one', 'integer',
     +     '0', '0', '2',
     +     'Elevation conversion flag',
     +     'Elevation conversion flag (0=none, 1=feet to meters,'//
     +     ' 2=meters to feet)',
     +     'none')/=0 ) CALL read_error(1, 'conv_flag')

      !ALLOCATE ( Tmax_allrain_sta(Nrain,12) )
      ! IF ( declparam(MODNAME, 'tmax_allrain_sta', 'nrain,nmonths',
      !+     'real', '38.0', '-8.0', '45.0',
      !+     'Precipitation is rain if HRU max temperature >= this value',
      !+     'Monthly (January to December) maximum air temperature'//
      !+     ' when precipitation is assumed to be rain; if HRU air'//
      !+     ' temperature is greater than or equal to this value,'//
      !+     ' precipitation is rain',
      !+     'temp_units')/=0 ) CALL read_error(1, 'tmax_allrain_sta')

      !ALLOCATE ( Tmax_allsnow_sta(Nrain,12) )
      ! IF ( declparam(MODNAME, 'tmax_allsnow_sta', 'nrain,nmonths',
      !+     'real', '32.0', '-10.0', '40.0',
      !+     'Maximum temperature when precipitation is all snow',
      !+     'Maximum air temperature when precipitation is assumed'//
      !+     ' to be snow; if HRU air temperature is less than or'//
      !+     ' equal to this value, precipitation is snow',
      !+     'temp_units')/=0 ) CALL read_error(1, 'tmax_allsnow_sta')

      IF ( declparam(MODNAME, 'tmax_allrain_dist', 'nmonths', 'real',
     +     '38.0', '-8.0', '75.0',
     +     'Precipitation is rain if station mean max temperature >='//
     +     ' this value',
     +     'Monthly (January to December) maximum air temperature'//
     +     ' when precipitation is assumed to be rain; if mean air'//
     +     ' temperature is greater than or equal to this value,'//
     +     ' precipitation is rain',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmax_allrain_dist')

      IF ( declparam(MODNAME, 'tmax_allsnow_dist', 'one', 'real',
     +     '32.0', '-10.0', '40.0',
     +     'Maximum air temperature when precipitation is all snow',
     +     'Maximum air temperature when precipitation is assumed'//
     +     ' to be snow; if mean air temperature is less than or'//
     +     ' equal to this value, precipitation is snow',
     +     'temp_units')/=0 ) CALL read_error(1, 'tmax_allsnow_dist')

      ALLOCATE ( Precip_xyz(Nrain) )

      END FUNCTION xyzdecl

!***********************************************************************
!     xyzinit - Initialize xyz_dist module - get parameter values,
!***********************************************************************
      INTEGER FUNCTION xyzinit()
      USE PRMS_XYZ_DIST
      USE PRMS_MODULE, ONLY: Nhru, Inputerror_flag, Ntemp, Nrain,
     +    Init_vars_from_file
      USE PRMS_BASIN, ONLY: Hru_area, Basin_area_inv,
     +    Hru_elev, Active_hrus, Hru_route_order, FEET2METERS
      USE PRMS_CLIMATEVARS, ONLY: Psta_elev, Tsta_elev
      IMPLICIT NONE
! Functions
      INTRINSIC DBLE
      INTEGER, EXTERNAL :: getparam
      EXTERNAL mean_by_month, read_error
! Local Variables
      INTEGER :: i, m, ii, ierr
!***********************************************************************
      xyzinit = 0

! Initialize declared variables
      IF ( Init_vars_from_file==0 ) THEN
        Is_rain_day = 0
        Tmax_rain_sta = 0.0
        Tmin_rain_sta = 0.0
      ENDIF

      IF ( getparam(MODNAME, 'adjust_rain', Nrain*12, 'real',
     +     Adjust_rain)/=0 ) CALL read_error(2, 'adjust_rain')

      IF ( getparam(MODNAME, 'adjust_snow', Nrain*12, 'real',
     +     Adjust_snow)/=0 ) CALL read_error(2, 'adjust_snow')

      IF ( getparam (MODNAME, 'solrad_elev', 1, 'real', Solrad_elev)
     +     /=0 ) CALL read_error(2, 'solrad_elev')

      IF ( getparam(MODNAME, 'hru_x', Nhru, 'real', MRUx)
     +     /=0 ) CALL read_error(2, 'hru_x')

      IF ( getparam(MODNAME, 'hru_y', Nhru, 'real', MRUy)
     +     /=0 ) CALL read_error(2, 'hru_y')

      IF ( getparam(MODNAME, 'max_lapse', MAXLAPSE*12, 'real',
     +     Max_lapse)/=0 ) CALL read_error(2, 'max_lapse')

      IF ( getparam(MODNAME, 'min_lapse', MAXLAPSE*12, 'real',
     +     Min_lapse)/=0 ) CALL read_error(2, 'min_lapse')

      IF ( getparam(MODNAME, 'ppt_lapse', MAXLAPSE*12, 'real',
     +     Ppt_lapse)/=0 ) CALL read_error(2, 'ppt_lapse')

      IF ( getparam(MODNAME, 'tsta_x', Ntemp, 'real', Temp_STAx)
     +     /=0 ) CALL read_error(2, 'tsta_x')

      IF ( getparam(MODNAME, 'tsta_y', Ntemp, 'real', Temp_STAy)
     +     /=0 ) CALL read_error(2, 'tsta_y')

      IF ( getparam(MODNAME, 'psta_x', Nrain, 'real', Rain_STAx)
     +     /=0 ) CALL read_error(2, 'psta_x')

      IF ( getparam(MODNAME, 'psta_y', Nrain, 'real', Rain_STAy)
     +     /=0 ) CALL read_error(2, 'psta_y')

      IF ( getparam(MODNAME, 'tsta_nuse', Ntemp, 'integer',
     +     Tsta_nuse)/=0 ) CALL read_error(2, 'tsta_nuse')

      IF ( getparam(MODNAME, 'psta_nuse', Nrain, 'integer',
     +     Psta_nuse)/=0 ) CALL read_error(2, 'psta_nuse')

      IF ( getparam(MODNAME, 'psta_freq_nuse', Nrain, 'integer',
     +     Psta_freq_nuse)/=0 ) CALL read_error(2, 'psta_freq_nuse')

      IF ( getparam(MODNAME, 'tsta_month_min', Ntemp*12, 'real',
     +     TminMTH)/=0 ) CALL read_error(2, 'tsta_month_min')

      IF ( getparam(MODNAME, 'tsta_month_max', Ntemp*12, 'real',
     +     TmaxMTH)/=0 ) CALL read_error(2, 'tsta_month_max')

      IF ( getparam(MODNAME, 'psta_month_ppt', Nrain*12, 'real',
     +     PptMTH)/=0 ) CALL read_error(2, 'psta_month_ppt')

      IF ( getparam(MODNAME, 'z_add', 1, 'real', Z_add)
     +     /=0 ) CALL read_error(2, 'z_add')

      IF ( getparam(MODNAME, 'z_div', 1, 'real', Z_div)
     +     /=0 ) CALL read_error(2, 'z_div')

      IF ( getparam(MODNAME, 'x_add', 1, 'real', X_add)
     +     /=0 ) CALL read_error(2, 'x_add')

      IF ( getparam(MODNAME, 'x_div', 1, 'real', X_div)
     +     /=0 ) CALL read_error(2, 'x_div')

      IF ( getparam(MODNAME, 'y_add', 1, 'real', Y_add)
     +     /=0 ) CALL read_error(2, 'y_add')

      IF ( getparam(MODNAME, 'y_div', 1, 'real', Y_div)
     +     /=0 ) CALL read_error(2, 'y_div')

      IF ( getparam(MODNAME, 'tmax_add', 1, 'real', Tmax_add)
     +     /=0 ) CALL read_error(2, 'tmax_add')

      IF ( getparam(MODNAME, 'tmax_div', 1, 'real', Tmax_div)
     +     /=0 ) CALL read_error(2, 'tmax_div')

      IF ( getparam(MODNAME, 'tmin_add', 1, 'real', Tmin_add)
     +     /=0 ) CALL read_error(2, 'tmin_add')

      IF ( getparam(MODNAME, 'tmin_div', 1, 'real', Tmin_div)
     +     /=0 ) CALL read_error(2, 'tmin_div')

      IF ( getparam(MODNAME, 'ppt_add', 1, 'real', Ppt_add)
     +     /=0 ) CALL read_error(2, 'ppt_add')

      IF ( getparam(MODNAME, 'ppt_div', 1, 'real', Ppt_div)
     +     /=0 ) CALL read_error(2, 'ppt_div')

      IF ( getparam(MODNAME, 'conv_flag', 1, 'integer', Conv_flag)
     +     /=0 ) CALL read_error(2, 'conv_flag')

      ! IF ( getparam(MODNAME, 'tmax_allrain_sta', Nrain*12, 'real',
      !+     Tmax_allrain_sta)/=0 ) CALL read_error(2, 'tmax_allrain_sta')
      ! IF ( getparam(MODNAME, 'tmax_allsnow_sta', Nrain*12, 'real',
      !+     Tmax_allsnow_sta)/=0 ) CALL read_error(2, 'tmax_allsnow_sta')

      IF ( getparam(MODNAME, 'tmax_allrain_dist', 12, 'real',
     +   Tmax_allrain_dist)/=0 ) CALL read_error(2, 'tmax_allrain_dist')

      IF ( getparam(MODNAME, 'tmax_allsnow_dist', 1, 'real',
     +   Tmax_allsnow_dist)/=0 ) CALL read_error(2, 'tmax_allsnow_dist')
!
! Compute basin centroid
!
      Basin_centroid_x = 0.0D0
      Basin_centroid_y = 0.0D0
      DO ii = 1, Active_hrus
        i = Hru_route_order(ii)
        Basin_centroid_x = Basin_centroid_x + DBLE(Hru_area(i)*MRUx(i))
        Basin_centroid_y = Basin_centroid_y + DBLE(Hru_area(i)*MRUy(i))
      ENDDO
      Basin_centroid_x = Basin_centroid_x*Basin_area_inv
      Basin_centroid_y = Basin_centroid_y*Basin_area_inv
!
! convert elevations from feet to meters
!
      IF ( Conv_flag==1 ) THEN
        DO i = 1, Nhru
          MRUelev(i) = Hru_elev(i)*FEET2METERS
        ENDDO
        DO i = 1, Ntemp
          Temp_STAelev(i) = Tsta_elev(i)*FEET2METERS
        ENDDO
        DO i = 1, Nrain
          Pstaelev(i) = Psta_elev(i)*FEET2METERS
        ENDDO
        Solradelev = Solrad_elev*FEET2METERS
      ELSE
        MRUelev = Hru_elev
        Temp_STAelev = Tsta_elev
        Pstaelev = Psta_elev
        Solradelev = Solrad_elev
      ENDIF
!
! transform Z, X and Y
!
! transform values only once
      Basin_centroid_x = (Basin_centroid_x+DBLE(X_add))/DBLE(X_div)
      Basin_centroid_y = (Basin_centroid_y+DBLE(Y_add))/DBLE(Y_div)
      Solradelev = (Solradelev+Z_add)/Z_div

      DO ii = 1, Active_hrus
        i = Hru_route_order(ii)
        MRUelev(i) = (MRUelev(i)+Z_add)/Z_div
        MRUx(i) = (MRUx(i)+X_add)/X_div
        MRUy(i) = (MRUy(i)+Y_add)/Y_div
      ENDDO

      ierr = 0
      Temp_nsta = 0
      Temp_nuse = 0
      DO i = 1, Ntemp
!       Temp_STAelev(i) = (Temp_STAelev(i)+Z_add)/Z_div
!       temp_STAx(i) = (temp_STAx(i)+X_add)/X_div
!       temp_STAy(i) = (temp_STAy(i)+Y_add)/Y_div
        IF ( Tsta_nuse(i)==1 ) THEN
          Temp_nsta = Temp_nsta + 1
          Temp_nuse(Temp_nsta) = i
        ENDIF
      ENDDO
      IF ( Temp_nsta<2 ) THEN
        PRINT *, 'ERROR, need to select at least 2 temperature stations'
        PRINT *, '       using tsta_nuse for xyz_dist'
        ierr = 1
      ENDIF

      ALLOCATE ( Pstax(Nrain), Pstay(Nrain) )
      Rain_nsta = 0
      Rain_nuse = 0
      DO i = 1, Nrain
        Pstaelev(i) = (Pstaelev(i)+Z_add)/Z_div
        Pstax(i) = (Rain_STAx(i)+X_add)/X_div
        Pstay(i) = (Rain_STAy(i)+Y_add)/Y_div
        IF ( Psta_nuse(i)==1 ) THEN
          Rain_nsta = Rain_nsta + 1
          Rain_nuse(Rain_nsta) = i
        ENDIF
      ENDDO
      IF ( Rain_nsta<2 ) THEN
        PRINT*,'ERROR, need to select at least 2 precipitation stations'
        PRINT *, '       using psta_nuse for xyz_dist'
        ierr = 1
      ENDIF
      IF ( ierr==1 ) THEN
        Inputerror_flag = 1
        RETURN
      ENDIF
!
! calculate the station mean by month
!
      DO m = 1, 12
        CALL mean_by_month(PptMTH(1, m), TminMTH(1, m), TmaxMTH(1, m),
     +                     Meanppt(m), Meantmax(m), Meantmin(m),
     +                     Rain_meanx(m), Rain_meany(m), Rain_meanz(m),
     +                     Temp_meanx(m), Temp_meany(m), Temp_meanz(m))
      ENDDO

      END FUNCTION xyzinit

!***********************************************************************
!     xyzrun - Temperature calculation
!               calculates daily max and min temperature
!               using data from available stations
!               Outputs a daily max and min Temperature by HRU elevation
!***********************************************************************
      INTEGER FUNCTION xyzrun()
      USE PRMS_XYZ_DIST
      USE PRMS_SET_TIME, ONLY: Nowmonth
      USE PRMS_OBS, ONLY: Rain_code
      IMPLICIT NONE
! Functions
      EXTERNAL xyz_temp_run, xyz_rain_run
!***********************************************************************
      xyzrun = 0

      CALL xyz_temp_run(Max_lapse(1, Nowmonth), Min_lapse(1, Nowmonth),
     +                  Meantmax(Nowmonth), Meantmin(Nowmonth),
     +                  Temp_meanx(Nowmonth), Temp_meany(Nowmonth),
     +                  Temp_meanz(Nowmonth))

      CALL xyz_rain_run(Ppt_lapse(1, Nowmonth), Rain_meanx(Nowmonth),
     +                  Rain_meany(Nowmonth), Rain_meanz(Nowmonth),
     +                  Meanppt(Nowmonth), Rain_code(Nowmonth))

      END FUNCTION xyzrun

!***********************************************************************
!     xyz_temp_run - Temperature calculation
!               calculates daily max and min temperature
!               using data from available stations
!               Outputs a daily max and min Temperature by HRU elevation
!***********************************************************************
      SUBROUTINE xyz_temp_run(Max_lapse, Min_lapse, Meantmax, Meantmin,
     +                        Temp_meanx, Temp_meany, Temp_meanz)
      USE PRMS_XYZ_DIST, ONLY: MRUx, MRUy, Tmax_rain_sta, Solradelev,
     +    Tmin_rain_sta, Temp_nuse, Tmin_add, Tmin_div, Tmax_add,
     +    Tmax_div, Temp_nsta, X_div, Y_div, Z_div, X_add, Y_add, Z_add,
     +    Temp_STAx, Temp_STAy, Basin_centroid_y, Basin_centroid_x,
     +    MAXLAPSE, Pstaelev, Pstax, Pstay, MRUelev, Temp_STAelev
      USE PRMS_MODULE, ONLY: Nrain
      USE PRMS_BASIN, ONLY: Basin_area_inv, Hru_area, Active_hrus,
     +    DNEARZERO, Hru_route_order
      USE PRMS_CLIMATEVARS, ONLY: Solrad_tmax, Solrad_tmin, Basin_temp,
     +    Basin_tmax, Basin_tmin, Tmaxf, Tminf, Tminc, Tmaxc, Tavgf,
     +    Tavgc, Tmin_aspect_adjust, Tmax_aspect_adjust
      USE PRMS_OBS, ONLY: Tmax, Tmin
      USE PRMS_SET_TIME, ONLY: Nowmonth
      IMPLICIT NONE
! Functions
      INTRINSIC ABS, SNGL, DBLE
      EXTERNAL temp_set
      REAL, EXTERNAL :: c_to_f
! Arguments
!   Declared Parameters
      REAL, INTENT(IN) :: Max_lapse(MAXLAPSE), Min_lapse(MAXLAPSE)
!   Private Variables
      DOUBLE PRECISION, INTENT(IN) :: Meantmax, Meantmin
      DOUBLE PRECISION, INTENT(IN) :: Temp_meanx, Temp_meany, Temp_meanz
! Local Variables
      INTEGER :: i, j, ntmin, ntmax, ii
      REAL :: tmax_hru, tmin_hru, x1
      DOUBLE PRECISION :: ztmax, ztmin
      DOUBLE PRECISION :: intmax, intmin, xmax, xmin, intercept
      DOUBLE PRECISION :: sumtmin, sumtmax, xtmin, xtmax, ytmin
      DOUBLE PRECISION :: ytmax, stmax, stmin
      DOUBLE PRECISION :: zrain, xrain, yrain
      DOUBLE PRECISION :: maxlapse1, maxlapse2, maxlapse3
      DOUBLE PRECISION :: minlapse1, minlapse2, minlapse3
!***********************************************************************
      sumtmin = 0.0D0
      sumtmax = 0.0D0
      ntmin = 0
      ntmax = 0
      xtmin = 0.0D0
      xtmax = 0.0D0
      ytmin = 0.0D0
      ytmax = 0.0D0
      ztmin = 0.0D0
      ztmax = 0.0D0
      maxlapse1 = DBLE( Max_lapse(1) )
      maxlapse2 = DBLE( Max_lapse(2) )
      maxlapse3 = DBLE( Max_lapse(3) )
      minlapse1 = DBLE( Min_lapse(1) )
      minlapse2 = DBLE( Min_lapse(2) )
      minlapse3 = DBLE( Min_lapse(3) )

! Do not Transform the coordinates of the temperature stations
! until after summing

      DO j = 1, Temp_nsta
        i = Temp_nuse(j)

        IF ( Tmax(i)>-99.0 ) THEN
          ntmax = ntmax + 1
!         sumtmax = sumtmax + ((Tmax(i)+Tmax_add)/Tmax_div)
          sumtmax = sumtmax + DBLE( Tmax(i) )
          xtmax = xtmax + DBLE( Temp_STAx(i) )
          ytmax = ytmax + DBLE( Temp_STAy(i) )
          ztmax = ztmax + DBLE( Temp_STAelev(i) )
        ENDIF

        IF ( Tmin(i)>-99.0 ) THEN
          ntmin = ntmin + 1
!         sumtmin = sumtmin + ((Tmin(i)+DBLE(Tmin_add))/DBLE(Tmin_div))
          sumtmin = sumtmin + DBLE( Tmin(i) )
          xtmin = xtmin + DBLE( Temp_STAx(i) )
          ytmin = ytmin + DBLE( Temp_STAy(i) )
          ztmin = ztmin + DBLE( Temp_STAelev(i) )
        ENDIF
      ENDDO
!
! calculate means
!
      IF ( ntmin>0 ) THEN
        stmin = sumtmin/DBLE(ntmin)
        xtmin = xtmin/DBLE(ntmin)
        ytmin = ytmin/DBLE(ntmin)
        ztmin = ztmin/DBLE(ntmin)
        stmin = (stmin+DBLE(Tmin_add))/DBLE(Tmin_div)
        xtmin = (xtmin+DBLE(X_add))/DBLE(X_div)
        ytmin = (ytmin+DBLE(Y_add))/DBLE(Y_div)
        ztmin = (ztmin+DBLE(Z_add))/DBLE(Z_div)
      ELSE
! these are already transformed
        stmin = Meantmin
        xtmin = Temp_meanx
        ytmin = Temp_meany
        ztmin = Temp_meanz
      ENDIF

      IF ( ntmax>0 ) THEN
        stmax = sumtmax/DBLE(ntmax)
        xtmax = xtmax/DBLE(ntmax)
        ytmax = ytmax/DBLE(ntmax)
        ztmax = ztmax/DBLE(ntmax)
        stmax = (stmax+DBLE(Tmax_add))/DBLE(Tmax_div)
        xtmax = (xtmax+DBLE(X_add))/DBLE(X_div)
        ytmax = (ytmax+DBLE(Y_add))/DBLE(Y_div)
        ztmax = (ztmax+DBLE(Z_add))/DBLE(Z_div)
      ELSE
! these are already transformed
        stmax = Meantmax
        xtmax = Temp_meanx
        ytmax = Temp_meany
        ztmax = Temp_meanz
      ENDIF
!
! adjust the values if not using all the stations
!
      IF ( ABS(Temp_meanz-ztmin)>DNEARZERO ) THEN
        intercept = stmin - minlapse3*ztmin - minlapse1*xtmin -
     +              minlapse2*ytmin
        stmin = minlapse3*Temp_meanz + minlapse1*Temp_meanx +
     +          minlapse2*Temp_meany + intercept
      ENDIF
!
      IF ( ABS(Temp_meanz-ztmax)>DNEARZERO ) THEN
        intercept = stmax - maxlapse3*ztmax - maxlapse1*xtmax -
     +              maxlapse2*ytmax
        stmax = maxlapse3*Temp_meanz + maxlapse1*Temp_meanx +
     +          maxlapse2*Temp_meany + intercept
      ENDIF
!
! now redistribute based on lapse rates redistribute to HRUs
!
      xmax = stmax
      xmin = stmin
!     --------------
      intmax = xmax -
     +         (maxlapse3*Temp_meanz) -
     +         (maxlapse1*Temp_meanx) -
     +         (maxlapse2*Temp_meany)
!
      intmin = xmin -
     +         (minlapse3*Temp_meanz) -
     +         (minlapse1*Temp_meanx) -
     +         (minlapse2*Temp_meany)
!
      Basin_tmax = 0.0D0
      Basin_tmin = 0.0D0
      Basin_temp = 0.0D0
!
!  Compute maximum temperature at XY centroid of basin at the elevation
!  of the solrad station used to develop the DD solrad curves.
!
      Solrad_tmax = ((maxlapse1*Basin_centroid_x)+
     +              (maxlapse2*Basin_centroid_y)+
     +              (maxlapse3*Solradelev+intmax))*DBLE(Tmax_div)
     +              - DBLE(Tmax_add)
      Solrad_tmin = ((minlapse1*Basin_centroid_x)+
     +              (minlapse2*Basin_centroid_y)+
     +              (minlapse3*Solradelev+intmin))*DBLE(Tmin_div)
     +              - DBLE(Tmin_add)

!
!  Compute temperatures at distribution stations.
!
      DO i = 1, Nrain
        zrain = DBLE( Pstaelev(i) )
        xrain = DBLE( Pstax(i) )
        yrain = DBLE( Pstay(i) )
        Tmax_rain_sta(i) = SNGL( (maxlapse1*xrain) +
     +                     (maxlapse2*yrain) +
     +                     (maxlapse3*zrain) + intmax )

        Tmin_rain_sta(i) = SNGL( (minlapse1*xrain) +
     +                     (minlapse2*yrain) +
     +                     (minlapse3*zrain) + intmin )

        Tmax_rain_sta(i) = (Tmax_rain_sta(i)*Tmax_div) - Tmax_add
        Tmin_rain_sta(i) = (Tmin_rain_sta(i)*Tmin_div) - Tmin_add

        IF ( Tmax_rain_sta(i)<Tmin_rain_sta(i) ) THEN
          x1 = Tmax_rain_sta(i)
          Tmax_rain_sta(i) = Tmin_rain_sta(i)
          Tmin_rain_sta(i) = x1
        ENDIF

      ENDDO

      DO ii = 1, Active_hrus
        i = Hru_route_order(ii)
!
!  At this point, all temperatures are in the units
!  of the temperatures in the data file.
!
        tmax_hru = (Max_lapse(1)*MRUx(i)) +
     +             (Max_lapse(2)*MRUy(i)) +
     +             (Max_lapse(3)*MRUelev(i)) + SNGL( intmax )

        tmin_hru = (Min_lapse(1)*MRUx(i)) +
     +             (Min_lapse(2)*MRUy(i)) +
     +             (Min_lapse(3)*MRUelev(i)) + SNGL( intmin )
!
!  Transform back
!
        tmax_hru = (tmax_hru*Tmax_div) - Tmax_add
        tmin_hru = (tmin_hru*Tmin_div) - Tmin_add
!
!  Temperature adjustment by HRU
!
        tmax_hru = tmax_hru + Tmax_aspect_adjust(i, Nowmonth)
        tmin_hru = tmin_hru + Tmin_aspect_adjust(i, Nowmonth)

!
!  If max is less than min, switch
!
        IF ( tmax_hru<tmin_hru ) THEN
          x1 = tmax_hru
          tmax_hru = tmin_hru
          tmin_hru = x1
        ENDIF
!
!  Now sort out units.
!
        CALL temp_set(i, tmax_hru, tmin_hru, Tmaxf(i), Tminf(i),
     +       Tavgf(i), Tmaxc(i), Tminc(i), Tavgc(i), Hru_area(i))
      ENDDO

      Basin_tmax = Basin_tmax*Basin_area_inv
      Basin_tmin = Basin_tmin*Basin_area_inv
      Basin_temp = Basin_temp*Basin_area_inv

      END SUBROUTINE xyz_temp_run

      SUBROUTINE xyz_rain_run(Ppt_lapse, Rain_meanx, Rain_meany,
     +                        Rain_meanz, Meanppt, Rain_code)
      USE PRMS_XYZ_DIST, ONLY: MRUx, MRUy, Rain_STAx, Rain_STAy,
     +    Rain_nuse, Ppt_add, Ppt_div, Rain_nsta, Tmax_rain_sta,
     +    Tmin_rain_sta, Is_rain_day, Psta_freq_nuse, X_div, Y_div,
     +    Z_div, X_add, Y_add, Z_add, Precip_xyz, MAXLAPSE, MRUelev,
     +    Tmax_allsnow_dist, Tmax_allrain_dist, Adjust_snow, Adjust_rain
      USE PRMS_MODULE, ONLY: Nrain
      USE PRMS_BASIN, ONLY: Hru_area, Basin_area_inv, Active_hrus,
     +    Hru_route_order, NEARZERO, DNEARZERO, MM2INCH
      USE PRMS_CLIMATEVARS, ONLY: Tmaxf, Tminf, Newsnow, Pptmix,
     +    Hru_ppt, Hru_rain, Hru_snow, Basin_rain,
     +    Basin_ppt, Prmx, Basin_snow, Psta_elev, Basin_obs_ppt,
     +    Precip_units, Tmax_allsnow_f, Adjmix_rain, Tmax_allrain_f
      USE PRMS_SET_TIME, ONLY: Nowmonth
      USE PRMS_OBS, ONLY: Precip, Rain_day
      IMPLICIT NONE
! Functions
      INTRINSIC ABS, SNGL, DBLE
      EXTERNAL precip_form
! Arguments
!   Declared Parameters
      INTEGER, INTENT(IN) :: Rain_code
      REAL, INTENT(IN) :: Ppt_lapse(MAXLAPSE)
!   Undeclared Static Variables
      DOUBLE PRECISION, INTENT(IN) :: Rain_meanx, Rain_meany
      DOUBLE PRECISION, INTENT(IN) :: Rain_meanz, Meanppt
! Local Variables
      INTEGER :: i, j, err_chk, nppt, nsta_used, ii
!dbg  INTEGER :: iy, id, im
      REAL :: ppt
      DOUBLE PRECISION :: intppt, intercept, sppt
      DOUBLE PRECISION :: sum_obs, sumppt, xppt, yppt, zppt
      DOUBLE PRECISION :: pptlapse1, pptlapse2, pptlapse3
!***********************************************************************
      pptlapse1 = DBLE( Ppt_lapse(1) )
      pptlapse2 = DBLE( Ppt_lapse(2) )
      pptlapse3 = DBLE( Ppt_lapse(3) )
!
! Code to check the rain_code parameter to determine if it is
! raining in the basin.
!
      Is_rain_day = 0
      IF ( Rain_code==1 ) THEN
        DO j = 1, Rain_nsta
          i = Rain_nuse(j)
          IF ( Precip(i)>0.0 ) Is_rain_day = 1
        ENDDO

      ELSEIF ( Rain_code==2 ) THEN
        DO i = 1, Nrain
          IF ( Precip(i)>0.0 ) Is_rain_day = 1
        ENDDO

      ELSEIF ( Rain_code==3 ) THEN
        Is_rain_day = 1

      ELSEIF ( Rain_code==4 ) THEN
        IF ( Rain_day==1 ) Is_rain_day = 1

      ELSEIF ( Rain_code==5 ) THEN
        DO i = 1, Nrain
          IF ( Psta_freq_nuse(i)==1 ) THEN
            IF ( Precip(i)>0.0 ) Is_rain_day = 1
          ENDIF
        ENDDO
      ENDIF
!
! add adjust_snow and adjust_rain here
!
      IF ( Rain_code==1 ) THEN
        nsta_used = Rain_nsta
      ELSE
        nsta_used = Nrain
      ENDIF

      Precip_xyz = Precip
      DO j = 1, nsta_used
        IF ( Rain_code==1 ) THEN
          i = Rain_nuse(j)
        ELSE
          i = j
        ENDIF

        IF ( Precip_xyz(i)<0.0 ) THEN
        !make sure negative precipitation values are < -99.0
          Precip_xyz(i) = -100.0
        ELSEIF ( Precip_xyz(i)>NEARZERO ) THEN
      !     IF ( Tmax_rain_sta(i)<=Tmax_allsnow_sta(i,Nowmonth) ) THEN
      !       err_chk = 1
      !     ELSEIF ( Tmin_rain_sta(i)>Tmax_allsnow_sta(i,Nowmonth) .OR.
      !+             Tmax_rain_sta(i)>=Tmax_allrain_sta(i,Nowmonth) ) THEN
      !       err_chk = 0
      !     ELSE
      !       err_chk = 1
      !     ENDIF
          IF ( Tmax_rain_sta(i)<=Tmax_allsnow_dist ) THEN
            err_chk = 1
          ELSEIF ( Tmin_rain_sta(i)>Tmax_allsnow_dist .OR.
     +             Tmax_rain_sta(i)>=Tmax_allrain_dist(Nowmonth) ) THEN
            err_chk = 0
          ELSE
            err_chk = 1
          ENDIF

          IF ( err_chk==1 ) THEN
            Precip_xyz(i) = (Precip_xyz(i)*Adjust_snow(i,Nowmonth))
     +                      + Precip_xyz(i)
          ELSE
            Precip_xyz(i) = (Precip_xyz(i)*Adjust_rain(i,Nowmonth))
     +                      + Precip_xyz(i)
          ENDIF
        ELSE !ignore very small precipitation values
          Precip_xyz(i) = 0.0
        ENDIF

      ENDDO

      sumppt = 0.0D0
      nppt = 0
      xppt = 0.0D0
      yppt = 0.0D0
      zppt = 0.0D0

      DO j = 1, Rain_nsta
        i = Rain_nuse(j)
        IF ( Precip_xyz(i)>=0.0 ) THEN
          nppt = nppt + 1
!         sumppt = sumppt + DBLE( (Precip_xyz(i)+Ppt_add)/Ppt_div) )
          sumppt = sumppt + DBLE( Precip_xyz(i) )
          xppt = xppt + DBLE( Rain_STAx(i) )
          yppt = yppt + DBLE( Rain_STAy(i) )
          zppt = zppt + DBLE( Psta_elev(i) )
        ENDIF
      ENDDO
!
! calculate means
!
      IF ( nppt>0 ) THEN
        sppt = sumppt/DBLE(nppt)
        xppt = xppt/DBLE(nppt)
        yppt = yppt/DBLE(nppt)
        zppt = zppt/DBLE(nppt)
        sppt = (sppt+DBLE(Ppt_add))/DBLE(Ppt_div)
        xppt = (xppt+DBLE(X_add))/DBLE(X_div)
        yppt = (yppt+DBLE(Y_add))/DBLE(Y_div)
        zppt = (zppt+DBLE(Z_add))/DBLE(Z_div)
      ELSE
! these are already transformed
        sppt = Meanppt
        xppt = Rain_meanx
        yppt = Rain_meany
        zppt = Rain_meanz
      ENDIF
!
! adjust the values if not using all the stations
!
      IF ( ABS(Rain_meanz-zppt)>DNEARZERO ) THEN

        intercept = (sppt) -
     +              (pptlapse3*zppt) -
     +              (pptlapse1*xppt) -
     +              (pptlapse2*yppt)
        sppt = (pptlapse3*Rain_meanz) +
     +         (pptlapse1*Rain_meanx) +
     +         (pptlapse2*Rain_meany) + intercept

      ENDIF

      xppt = sppt

      intppt = xppt -
     +         (pptlapse3*Rain_meanz) -
     +         (pptlapse1*Rain_meanx) -
     +         (pptlapse2*Rain_meany)

      Basin_ppt = 0.0D0
      Basin_rain = 0.0D0
      Basin_snow = 0.0D0

      sum_obs = 0.0D0
      DO ii = 1, Active_hrus
        i = Hru_route_order(ii)
        Pptmix(i) = 0
        Newsnow(i) = 0
        Hru_ppt(i) = 0.0
        Hru_rain(i) = 0.0
        Hru_snow(i) = 0.0
        Prmx(i) = 0.0

        IF ( Is_rain_day/=0 ) THEN
          ppt = (Ppt_lapse(1)*MRUx(i)) + (Ppt_lapse(2)*MRUy(i)) +
     +          (Ppt_lapse(3)*MRUelev(i)) + SNGL( intppt )

          ppt = ppt*Ppt_div - Ppt_add

!******Ignore small amounts of precipitation on HRU
          IF ( ppt>NEARZERO ) THEN
            IF ( Precip_units==1 ) ppt = ppt*MM2INCH
            CALL precip_form(ppt, Hru_ppt(i), Hru_rain(i),
     +           Hru_snow(i), Tmaxf(i), Tminf(i), Pptmix(i),
     +           Newsnow(i), Prmx(i), Tmax_allrain_f(i,Nowmonth), 1.0,
     +           1.0, Adjmix_rain(i,Nowmonth), Hru_area(i), sum_obs,
     +           Tmax_allsnow_f(i,Nowmonth))

          ENDIF
        ENDIF
      ENDDO
      Basin_ppt = Basin_ppt*Basin_area_inv
      Basin_rain = Basin_rain*Basin_area_inv
      Basin_snow = Basin_snow*Basin_area_inv
      Basin_obs_ppt = sum_obs*Basin_area_inv

      END SUBROUTINE xyz_rain_run

!***********************************************************************
!     calculate the station mean by month
!***********************************************************************
      SUBROUTINE mean_by_month (PptMTH, TminMTH, TmaxMTH, Meanppt,
     +                          Meantmax, Meantmin, Rain_meanx,
     +                          Rain_meany, Rain_meanz, Temp_meanx,
     +                          Temp_meany, Temp_meanz)
      USE PRMS_XYZ_DIST, ONLY: Rain_nsta, Rain_nuse, Temp_nsta,
     +    Temp_nuse, Rain_STAx, Rain_STAy, Temp_STAx, Tmax_div,
     +    Temp_STAy, Ppt_add, Ppt_div, Tmin_add, Tmin_div, Tmax_add,
     +    X_div, Y_div, Z_div, X_add, Y_add, Z_add, Temp_STAelev
      USE PRMS_MODULE, ONLY: Ntemp, Nrain
      USE PRMS_CLIMATEVARS, ONLY: Psta_elev
      IMPLICIT NONE
      INTRINSIC DBLE
! Arguments
      REAL, INTENT(IN) :: PptMTH(Nrain), TminMTH(Ntemp), TmaxMTH(Ntemp)
      DOUBLE PRECISION, INTENT(OUT) :: Meanppt, Meantmax, Meantmin
      DOUBLE PRECISION, INTENT(OUT) :: Rain_meanx, Rain_meany,Rain_meanz
      DOUBLE PRECISION, INTENT(OUT) :: Temp_meanx, Temp_meany,Temp_meanz
! Local Variables
      INTEGER :: i, j
      DOUBLE PRECISION :: rain_n, temp_n
!***********************************************************************
      rain_n = 1.0D0/DBLE(Rain_nsta)
      temp_n = 1.0D0/DBLE(Temp_nsta)

      Meanppt = 0.0D0
      Meantmax = 0.0D0
      Meantmin = 0.0D0
      Rain_meanx = 0.0D0
      Rain_meany = 0.0D0
      Rain_meanz = 0.0D0
      Temp_meanx = 0.0D0
      Temp_meany = 0.0D0
      Temp_meanz = 0.0D0

      DO j = 1, Rain_nsta
        i = Rain_nuse(j)

!       PptMTH(i) = (PptMTH(i)+Ppt_add)/Ppt_div

        Meanppt = Meanppt + DBLE( PptMTH(i) )
        Rain_meanx = Rain_meanx + DBLE( Rain_STAx(i) )
        Rain_meany = Rain_meany + DBLE( Rain_STAy(i) )
        Rain_meanz = Rain_meanz + DBLE( Psta_elev(i) )
      ENDDO

      DO j = 1, Temp_nsta
        i = Temp_nuse(j)

!       TminMTH(i) = (TminMTH(i)+Tmin_add)/Tmin_div
!       TmaxMTH(i) = (TmaxMTH(i)+Tmax_add)/Tmax_div

        Meantmin = Meantmin + DBLE( TminMTH(i) )
        Meantmax = Meantmax + DBLE( TmaxMTH(i) )
        Temp_meanx = Temp_meanx + DBLE( Temp_STAx(i) )
        Temp_meany = Temp_meany + DBLE( Temp_STAy(i) )
        Temp_meanz = Temp_meanz + DBLE( Temp_STAelev(i) )
      ENDDO

      Meanppt = Meanppt*rain_n
      Meantmin = Meantmin*temp_n
      Meantmax = Meantmax*temp_n

      Rain_meanx = Rain_meanx*rain_n
      Rain_meany = Rain_meany*rain_n
      Rain_meanz = Rain_meanz*rain_n

      Temp_meanx = Temp_meanx*temp_n
      Temp_meany = Temp_meany*temp_n
      Temp_meanz = Temp_meanz*temp_n

      Meanppt = (Meanppt+DBLE(Ppt_add))/DBLE(Ppt_div)
      Meantmin = (Meantmin+DBLE(Tmin_add))/DBLE(Tmin_div)
      Meantmax = (Meantmax+DBLE(Tmax_add))/DBLE(Tmax_div)

      Rain_meanx = (Rain_meanx+DBLE(X_add))/DBLE(X_div)
      Rain_meany = (Rain_meany+DBLE(Y_add))/DBLE(Y_div)
      Rain_meanz = (Rain_meanz+DBLE(Z_add))/DBLE(Z_div)

      Temp_meanx = (Temp_meanx+X_add)/DBLE(X_div)
      Temp_meany = (Temp_meany+Y_add)/DBLE(Y_div)
      Temp_meanz = (Temp_meanz+Z_add)/DBLE(Z_div)

      END SUBROUTINE mean_by_month

!***********************************************************************
!     xyz_restart - write or read xyz restart file
!***********************************************************************
      SUBROUTINE xyz_restart(In_out)
      USE PRMS_MODULE, ONLY: Restart_outunit, Restart_inunit
      USE PRMS_XYZ_DIST
      IMPLICIT NONE
      ! Argument
      INTEGER, INTENT(IN) :: In_out
      EXTERNAL check_restart
      ! Local Variable
      CHARACTER(LEN=8) :: module_name
!***********************************************************************
      IF ( In_out==0 ) THEN
        WRITE ( Restart_outunit ) MODNAME
        WRITE ( Restart_outunit ) Is_rain_day
        WRITE ( Restart_outunit ) Tmax_rain_sta
        WRITE ( Restart_outunit ) Tmin_rain_sta
      ELSE
        READ ( Restart_inunit ) module_name
        CALL check_restart(MODNAME, module_name)
        READ ( Restart_inunit ) Is_rain_day
        READ ( Restart_inunit ) Tmax_rain_sta
        READ ( Restart_inunit ) Tmin_rain_sta
      ENDIF
      END SUBROUTINE xyz_restart

max_stars_repo_name: iyanuobidele/prms-on-pegasus

max_issues_repo_name: iyanuobidele/prms-on-pegasus

max_forks_repo_name: iyanuobidele/prms-on-pegasus


text:
! Sub-routine that computes the Hessian matrix of a function,
! through the finite difference method.
! Igor Lopes, February 2015
SUBROUTINE HESSIAN(H,X,NDIM)
    IMPLICIT NONE
    ! Parameters
    REAL(8) R1 /1.0D0/
    REAL(8) R2 /2.0D0/
    REAL(8) R4 /4.0D0/
    REAL(8) PERT /1.0D-4/
    ! Arguments
    INTEGER     :: NDIM
    REAL(8),DIMENSION(NDIM) :: X
    REAL(8),DIMENSION(NDIM,NDIM) :: H
    ! Locals
    INTEGER :: I,J
    REAL(8) :: FO,F1,F2,F3,F4,EVALFUNC
    REAL(8),DIMENSION(NDIM) :: X1,X2,X3,X4
    ! Algorithm
    FO=EVALFUNC(X,NDIM)
    DO I=1,NDIM
        DO J=1,NDIM
            IF(I.EQ.J)THEN
                X1=X
                X2=X
                X1(I)=X(I)+PERT
                X2(I)=X(I)-PERT
                F1=EVALFUNC(X1,NDIM)
                F2=EVALFUNC(X2,NDIM)
                H(I,I)=(F1-R2*FO+F2)/(PERT*PERT)
            ELSE
                X1=X
                X2=X
                X3=X
                X4=X
                X1(I)=X(I)+PERT
                X1(J)=X(J)+PERT
                X2(I)=X(I)+PERT
                X2(J)=X(J)-PERT
                X3(I)=X(I)-PERT
                X3(J)=X(J)+PERT
                X4(I)=X(I)-PERT
                X4(J)=X(J)-PERT
                F1=EVALFUNC(X1,NDIM)
                F2=EVALFUNC(X2,NDIM)
                F3=EVALFUNC(X3,NDIM)
                F4=EVALFUNC(X4,NDIM)
                H(I,J)=(F1-F2-F3+F4)/(R4*PERT*PERT)
            ENDIF
        ENDDO
    ENDDO
END SUBROUTINE

max_stars_repo_name: iarlopes/GPROPT

max_issues_repo_name: iarlopes/GPROPT

max_forks_repo_name: iarlopes/GPROPT


text:
	SUBROUTINE SFXCHR  ( iloc, parms, icolor, cdata, ntime, xval,
     +			     iret )
C************************************************************************
C* SFXCHR								*
C*									*
C* This subroutine plots character data for SFGRAM.			*
C*									*
C* SFXCHR  ( ILOC, PARMS, ICOLOR, CDATA, NTIME, XVAL, IRET )		*
C*									*
C* Input parameters:							*
C*	ILOC		INTEGER		Location on axis (1,2,3)	*
C*	PARMS		CHAR*		Parm*condition			*
C*	ICOLOR		INTEGER		Color				*
C*	CDATA (NTIME)	REAL		Character data			*
C*	NTIME		INTEGER		Number of times			*
C*	XVAL (NTIME)	REAL		Points on x axis		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 5/90						*
C* S. Jacobs/NCEP	 9/97	Changed call to GQTEXT and GSTEXT	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
C*
	CHARACTER*(*)	parms, cdata (*)
	REAL		xval (*)
C*
	CHARACTER	condtn*8, ppp*4
C*
C------------------------------------------------------------------------
	iret = 0
C
C*	Get the condition and parameter parts.
C
	condtn = parms (5: )
	ppp    = parms
C
C*	Set color.
C
	CALL GSCOLR  ( icolor, ier )
C
C*	Get value on y axis.
C
	CALL SFXLOC  ( iloc, 'C', yval, ier )
C
C*	Get conditions.
C
	CALL SFXCND  ( condtn, size, iwidth, ier )
C
C*	Set the character size and width.
C
	CALL GQTEXT  ( itxfn, itxhw, sztext, itxwid,
     +		       ibrdr, irrotn, ijust, ier )
	CALL GSTEXT  ( 0, 0, size, iwidth, 0, 0, 0, ier )
C
C*	Loop through points writing text.
C
	DO  i = 1, ntime
	    IF  ( cdata (i) .ne. ' ' )  THEN
		CALL GTEXT  ( 'M', xval (i), yval, cdata (i), 90.,
     +			       0, 0, ier )
	    END IF
	END DO
C
C*	Reset text attributes.
C
	CALL GSTEXT  ( 0, 0, sztext, itxwid, 0, 0, 0, ier )
C*
	RETURN
	END