! -*- Mode: Fortran; -*-
!
!  (C) 2003 by Argonne National Laboratory.
!      See COPYRIGHT in top-level directory.
!
! This tests spawn_mult by using the same executable and no command-line
! options.  The attribute MPI_APPNUM is used to determine which
! executable is running.
!
       program main
       use mpi_f08
       integer (kind=MPI_ADDRESS_KIND) aint

       integer errs, err
       integer rank, size, rsize, wsize, i
       integer np(2)
       type(MPI_Info) infos(2)
       integer errcodes(2)
       type(MPI_Comm) parentcomm, intercomm
       type(MPI_Status) status
       character*(30) cmds(2)
       integer appnum
       logical flag
       integer ierr
       integer can_spawn

       errs = 0

       call MTest_Init( ierr )

       call MTestSpawnPossible( can_spawn, errs )
        if ( can_spawn .eq. 0 ) then
            call MTest_Finalize( errs )
            goto 300
        endif

       call MPI_Comm_get_parent( parentcomm, ierr )

       if (parentcomm .eq. MPI_COMM_NULL) then
!       Create 2 more processes
           cmds(1) = "./spawnmult2f90"
           cmds(2) = "./spawnmult2f90"
           np(1)   = 1
           np(2)   = 1
           infos(1)= MPI_INFO_NULL
           infos(2)= MPI_INFO_NULL
           call MPI_Comm_spawn_multiple( 2, cmds, MPI_ARGVS_NULL,            &
      &             np, infos, 0,                                             &
      &             MPI_COMM_WORLD, intercomm, errcodes, ierr )
        else
           intercomm = parentcomm
        endif

!       We now have a valid intercomm

        call MPI_Comm_remote_size( intercomm, rsize, ierr )
        call MPI_Comm_size( intercomm, size, ierr )
        call MPI_Comm_rank( intercomm, rank, ierr )

        if (parentcomm .eq. MPI_COMM_NULL) then
!           Master
            if (rsize .ne. np(1) + np(2)) then
                errs = errs + 1
                print *, "Did not create ", np(1)+np(2),                    &
      &          " processes (got ", rsize, ")"
            endif
! Allow a multi-process parent
            if (rank .eq. 0) then
               do i=0, rsize-1
                  call MPI_Send( i, 1, MPI_INTEGER, i, 0, intercomm,        &
      &                 ierr )
               enddo
!       We could use intercomm reduce to get the errors from the
!       children, but we'll use a simpler loop to make sure that
!       we get valid data
              do i=0, rsize-1
                 call MPI_Recv( err, 1, MPI_INTEGER, i, 1, intercomm,       &
      &                       MPI_STATUS_IGNORE, ierr )
                 errs = errs + err
              enddo
            endif
         else
!       Child
!       FIXME: This assumes that stdout is handled for the children
!       (the error count will still be reported to the parent)
           if (size .ne. 2) then
              errs = errs + 1
              print *, "(Child) Did not create ", 2,                        &
      &             " processes (got ",size, ")"
              call mpi_comm_size( MPI_COMM_WORLD, wsize, ierr )
              if (wsize .eq. 2) then
                 errs = errs + 1
                 print *, "(Child) world size is 2 but ",                   &
      &          " local intercomm size is not 2"
              endif
           endif

         call MPI_Recv( i, 1, MPI_INTEGER, 0, 0, intercomm, status, ierr    &
      &     )
         if (i .ne. rank) then
            errs = errs + 1
            print *, "Unexpected rank on child ", rank, "(",i,")"
         endif
!
!       Check for correct APPNUM
         call mpi_comm_get_attr( MPI_COMM_WORLD, MPI_APPNUM, aint,         &
      &        flag, ierr )
!        My appnum should be my rank in comm world
         if (flag) then
            appnum = aint
            if (appnum .ne. rank) then
                errs = errs + 1
                print *, "appnum is ", appnum, " but should be ", rank
             endif
         else
             errs = errs + 1
             print *, "appnum was not set"
         endif

!       Send the errs back to the master process
        call MPI_Ssend( errs, 1, MPI_INTEGER, 0, 1, intercomm, ierr )
        endif

!       It isn't necessary to free the intercomm, but it should not hurt
        call MPI_Comm_free( intercomm, ierr )

!       Note that the MTest_Finalize get errs only over COMM_WORLD
        if (parentcomm .eq. MPI_COMM_NULL) then
            call MTest_Finalize( errs )
        endif

 300    continue
        call MPI_Finalize( ierr )
        end


text:
!-----------------------------------------------------------------------
!     Subroutine:    fieldlines_init_hint
!     Authors:       S. Lazerson (samuel.lazerson@ipp.mpg.de)
!     Date:          07/27/2021
!     Description:   This subroutine handles setting up the magnetic
!                    field based on HINT equilibira.
!-----------------------------------------------------------------------
      SUBROUTINE fieldlines_init_hint
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE fieldlines_runtime, extcur_fieldlines => extcur
      USE fieldlines_grid, ONLY: raxis_g => raxis, phiaxis, &
                                 zaxis_g => zaxis, nr, nphi, nz, &
                                 rmin, rmax, zmin, zmax, phimin, &
                                 phimax, vc_adapt_tol, B_R, B_Z, B_PHI,&
                                 BR_spl, BZ_spl
      USE read_hint_mod, ONLY: get_hint_grid, get_hint_B, get_hint_press, &
                               get_hint_maxp, read_hint_deallocate, &
                               get_hint_magaxis, get_hint_gridB
      USE mpi_params                                                    ! MPI
!      USE mpi
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER, PARAMETER :: BYTE_8 = SELECTED_INT_KIND (8)
#if defined(MPI_OPT)
      INTEGER(KIND=BYTE_8),ALLOCATABLE :: mnum(:), moffsets(:)
      INTEGER :: numprocs_local, mylocalid, mylocalmaster
      INTEGER :: MPI_COMM_LOCAL
#endif
      INTEGER :: ier, s, i, j, k, mystart, myend
      REAL(rprec) :: br, bphi, bz
      INTEGER :: nrh,nzh,nph
      REAL(rprec) :: rmin_hint, rmax_hint, zmin_hint, zmax_hint, &
                     pmax_hint, pres_max, betatot
!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      ! Handle what we do

      ! Divide up Work
#if defined(MPI_OPT)
      CALL MPI_COMM_DUP( MPI_COMM_SHARMEM, MPI_COMM_LOCAL, ierr_mpi)
      CALL MPI_COMM_RANK( MPI_COMM_LOCAL, mylocalid, ierr_mpi )              ! MPI
      CALL MPI_COMM_SIZE( MPI_COMM_LOCAL, numprocs_local, ierr_mpi )          ! MPI
#endif
      mylocalmaster = master

      ! Get the maximum pressure
      CALL get_hint_maxp(pres_max)

      ! Write info to screen
      IF (lverb) THEN
         betatot = 0
         CALL get_hint_grid(nrh,nzh,nph,rmin_hint,rmax_hint,zmin_hint,zmax_hint,pmax_hint)
         WRITE(6,'(A)')               '----- HINT Information -----'
         WRITE(6,'(A,F9.5,A,F9.5,A,I4)') '   R   = [',rmin_hint,',',rmax_hint,'];  NR:   ',nrh
         WRITE(6,'(A,F8.5,A,F8.5,A,I4)') '   PHI = [',0.0,',',pmax_hint,'];  NPHI: ',nph
         WRITE(6,'(A,F8.5,A,F8.5,A,I4)') '   Z   = [',zmin_hint,',',zmax_hint,'];  NZ:   ',nzh
         WRITE(6,'(A,F7.2,A)')           '   PRES_MAX = ',pres_max*1E-3,' [kPa]'
      END IF

      IF (lverb) THEN
         WRITE(6,'(5X,A,I3.3,A)',ADVANCE='no') 'Plasma Field Calculation [',0,']%'
         CALL FLUSH(6)
      END IF

      ! Break up the Work
      CALL MPI_CALC_MYRANGE(MPI_COMM_LOCAL,1, nr*nphi*nz, mystart, myend)

      IF (lafield_only) THEN
      ELSE
         DO s = mystart, myend
            i = MOD(s-1,nr)+1
            j = MOD(s-1,nr*nphi)
            j = FLOOR(REAL(j) / REAL(nr))+1
            k = CEILING(REAL(s) / REAL(nr*nphi))

            ! Bfield
            CALL get_hint_gridB(i,j,k,br,bphi,bz)
            B_R(i,j,k)   = B_R(i,j,k)   + br
            B_PHI(i,j,k) = B_PHI(i,j,k) + bphi
            B_Z(i,j,k)   = B_Z(i,j,k)   + bz

            IF (lverb .and. (MOD(s,nr) == 0)) THEN
               CALL backspace_out(6,6)
               WRITE(6,'(A,I3,A)',ADVANCE='no') '[',INT((100.*s)/(myend-mystart+1)),']%'
               CALL FLUSH(6)
            END IF
         END DO
      END IF

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi)
#endif

      ! Free variables
      CALL read_hint_deallocate

      IF (lverb) THEN
         CALL backspace_out(6,36)
         CALL FLUSH(6)
         WRITE(6,'(36X)',ADVANCE='no')
         CALL FLUSH(6)
         CALL backspace_out(6,36)
         CALL FLUSH(6)
      END IF

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi)
      CALL MPI_COMM_FREE(MPI_COMM_LOCAL,ierr_mpi)
      CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'beams3d_init_vmec',ierr_mpi)
#endif
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE fieldlines_init_hint

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
	SUBROUTINE TI_DTM4  ( dattim, dattm4, iret )
C************************************************************************
C* TI_DTM4								*
C*									*
C* This subroutine modifies a standard GEMPAK time to include a four-   *
C* digit year, instead of a 2-digit year.  Any 2-digit year less than or*
C* equal to 40 will be assumed to be in  the 21st century; years greater*
C* than 40 will be assumed to be in the 20th century.                   *
C*									*
C* TI_DTM4  ( DATTIM, DATTM4, IRET )					*
C*									*
C* Input parameters:							*
C*	DATTIM		CHAR*		GEMPAK time			*
C*									*
C* Output parameters:							*
C*	DATTM4    	CHAR*  		Time with 4-digit year          *
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -1 = invalid time              *
C*					 -7 = invalid year		*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 3/99               				*
C* D. Kidwell/NCEP	 4/99 	Added 4-digit year check; fixed prologue*
C* B. Hebbard/NCEP       3/18   Moved century break from 2020 to 2040   *
C************************************************************************
	CHARACTER*(*)	dattim, dattm4
C*
	CHARACTER	cent*2
C------------------------------------------------------------------------
	iret   = 0
	dattm4 = ' '
C
C*	Check to see if year is already four digits.
C
	CALL ST_LSTR ( dattim, length, ier )
	islash = INDEX ( dattim, '/' )
C
	IF ( ( length .ge. 13 ) .and. ( islash .eq. 9 ) ) THEN
	    dattm4 = dattim
	  ELSE
	    CALL TI_CCNT ( dattim, cent, iret )
	    IF ( iret .eq. 0 ) THEN
	        dattm4 = cent // dattim
	    END IF
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of Sandia Corporation nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C$Id: dbpini.f,v 1.1 1999/02/16 21:37:59 gdsjaar Exp $
C=======================================================================
      SUBROUTINE DBPINI (OPTION, NDB, TITLE, NDIM, NUMNP, NUMEL, NELBLK,
     &                   NUMNPS, LNPSNL, LNPSDF, NUMESS, LESSEL, LESSNL,
     &                   LESSDF, NVARGL, NVARNP, NVAREL, FILNAM)
C=======================================================================
C   --*** DBPINI *** (EXOLIB) Display database title and initial variables
C   --   Written by Amy Gilkey - revised 10/15/87
C   -- Revised 8/5/95 to work with EXODUSIIV2 database format
C   --
C   --DBPINI displays the database filename (optional) and title and
C   --the database initial variables.
C   --
C   --Parameters:
C   --   OPTION - IN - '*' to print all, else print options:
C   --      'N' to print database name
C   --      'T' to print title
C   --      'I' to print number of nodes, etc.
C   --      'S' to print node set and side set information
C   --      'V' to print number of variables
C   --   NDB    - IN - the database file (if OPTION)
C   --   TITLE  - IN - the database title (if OPTION)
C   --   NDIM   - IN - the number of coordinates per node (if OPTION)
C   --   NUMNP  - IN - the number of nodes (if OPTION)
C   --   NUMEL  - IN - the number of elements (if OPTION)
C   --   NELBLK - IN - the number of element blocks (if OPTION)
C   --   NUMNPS - IN - the number of node sets (if OPTION)
C   --   LNPSNL - IN - the length of the node sets node list (if OPTION)
C   --   LNPSDF - IN - the length of the node set distribution list
C   --   NUMESS - IN - the number of side sets (if OPTION)
C   --   LESSEL - IN - the length of the side sets element list
C   --                 (if OPTION)
C   --   LESSNL - IN - the length of the side sets node list (if OPTION)
C   --   NVARGL - IN - the number of global variables (if OPTION)
C   --   NVARNP - IN - the number of nodal variables (if OPTION)
C   --   NVAREL - IN - the number of element variables (if OPTION)

      include 'exodusII.inc'

      CHARACTER*(*) FILNAM
      CHARACTER*(*) OPTION
      CHARACTER*(MXLNLN) TITLE
      INTEGER NDB
      INTEGER NDIM, NUMNP, NUMEL, NELBLK, NUMNPS,
     &        LNPSNL, LNPSDF, NUMESS, LESSEL, LESSNL, LESSDF
      LOGICAL ALL
      ALL = (OPTION .EQ. '*')

      IF (ALL .OR. (INDEX (OPTION, 'I') .GT. 0)) THEN
         WRITE(*, 10000)
      ENDIF

      IF (ALL .OR. (INDEX (OPTION, 'N') .GT. 0)) THEN
         WRITE (*, 10010) FILNAM(:LENSTR(FILNAM))
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'T') .GT. 0)) THEN
         WRITE (*, 10020) TITLE(1:LENSTR(TITLE))
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'I') .GT. 0)) THEN
         WRITE (*, 10030, IOSTAT=IDUM)
     &      NDIM, NUMNP, NUMEL, NELBLK
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'S') .GT. 0)) THEN
         If (NUMNPS .LE. 0) THEN
            WRITE (*, 10040, IOSTAT=IDUM) NUMNPS
         ELSE
            WRITE (*, 10040, IOSTAT=IDUM) NUMNPS, LNPSNL, LNPSDF
         END IF
         IF (NUMESS .LE. 0) THEN
            WRITE (*, 10050, IOSTAT=IDUM) NUMESS
         ELSE
            WRITE (*, 10050, IOSTAT=IDUM)
     &             NUMESS, LESSEL, LESSNL, LESSDF
         END IF
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'V') .GT. 0)) THEN
         WRITE (*, 10060, IOSTAT=IDUM) NVARGL, NVARNP, NVAREL
      END IF

10000 FORMAT (/, 1X, 'DATABASE INITIAL VARIABLES')
10010 FORMAT (/, 1X, 'Database:  ', A)
10020 FORMAT (/, 1X, A)
10030 FORMAT (
     &      /, 1X, 'Number of coordinates per node         =', I10
     &      /, 1X, 'Number of nodes                        =', I10
     &      /, 1X, 'Number of elements                     =', I10
     &      /, 1X, 'Number of element blocks               =', I10
     &      )
10040 FORMAT (
     &      /, 1X, 'Number of node sets                    =', I10, :
     &      /, 1X, '   Length of node list                 =', I10
     &      /, 1X, '   Length of distribution factors list =', I10
     &      )
10050 FORMAT
     &      (  1X, 'Number of side sets                    =', I10, :
     &      /, 1X, '   Length of element list              =', I10
     &      /, 1X, '   Length of node list                 =', I10
     &      /, 1X, '   Length of distribution factors list =', I10
     &      )
10060 FORMAT (
     &      /, 1X, 'Number of global variables             =', I10
     &      /, 1X, 'Number of variables at each node       =', I10
     &      /, 1X, 'Number of variables at each element    =', I10
     &      )

      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
	subroutine partials(fn_srcpar,
     &	nsrc, src_cusp, src_lat, src_lon, src_dep,
     &	nsta, sta_lab, sta_lat, sta_lon,
     &	mod_nl, mod_ratio, mod_v, mod_top,
     &	tmp_ttp, tmp_tts,
     &	tmp_xp, tmp_yp, tmp_zp)

c Compute partial derivatives for straight ray paths (option for 1D layered
c models commented out). Velocity is taken from first entry in the model
c specification array in hypoDD.inp.

	implicit none

	include'hypoDD.inc'

c	Parameters:
	character	fn_srcpar*80	! Source-parameter file
	integer		nsrc		! No. of sources
	integer		src_cusp(MAXEVE)! [1..nsrc]
	doubleprecision	src_lat(MAXEVE)	! [1..nsrc]
	doubleprecision	src_lon(MAXEVE)	! [1..nsrc]
	real		src_dep(MAXEVE)	! [1..nsrc]
	integer		nsta		! No. of stations
	character	sta_lab(MAXSTA)*7! [1..nsta]
	real		sta_lat(MAXSTA)	! [1..nsta]
	real		sta_lon(MAXSTA)	! [1..nsta]
	integer		mod_nl		! No. of layers
	real		mod_ratio	! Vp/Vs
	real		mod_v(MAXLAY)	! [1..mod_nl]
	real		mod_top(MAXLAY)	! [1..mod_nl]
	real		tmp_ttp(MAXSTA,MAXEVE)! [1..nsta,1..nsrc]
	real		tmp_tts(MAXSTA,MAXEVE)! [1..nsta,1..nsrc]
	real		tmp_xp(MAXSTA,MAXEVE)! [1..nsta,1..nsrc]
	real		tmp_yp(MAXSTA,MAXEVE)! [1..nsta,1..nsrc]
	real		tmp_zp(MAXSTA,MAXEVE)! [1..nsta,1..nsrc]

c	Local variables:
	real		ain
	real		az
	real		del
	real		dist
	integer		i, j, k
	integer		iunit		! Output unit number
	real		pi
	integer		trimlen
	real		vs(20)

	parameter(pi=3.141593)

	character rcsid*150
	data rcsid /"$Header: /home1/crhet/julian/HYPODD/src/hypoDD/RCS/partials.f,v 1.7 2001/02/17 23:54:57 julian Exp julian $"/
	save rcsid

      iunit = 0
      if (trimlen(fn_srcpar).gt.1) then
c        Open source-parameter file
         call freeunit(iunit)
         open(iunit,file=fn_srcpar,status='unknown')
      endif

c     Make sure hypocenters don't fall on layer boundaries
      do i=1,nsrc
         do j=1,mod_nl
            if (abs(src_dep(i)-mod_top(j)).lt.0.0001)
     &         src_dep(i) = src_dep(i)-0.001
         enddo
      enddo

c     Get S velocity model
      do i=1,mod_nl
         vs(i) = mod_v(i)/mod_ratio
      enddo

c     Compute epicentral distances, azimuths, angles of incidence,
c     and P/S-travel times from sources to stations

c the following loop needs some more work. I replaced the ray tracing
c with a straight ray path approach, which you can probably get away in
c your case. However, you need to get the geometry between source and receiver
c straight; i.e. you need to include receiver depth (probably easiest to
c shift the source by the amount of receiver depth) to get proper incidence
c angle (ain) and distance (dist). I think this is the only place you need to
c change anything. of course, any ray trace can be plugged in here.
      do i=1,nsta
         do j=1,nsrc
            call delaz2(src_lat(j), src_lon(j), sta_lat(i), sta_lon(i),
     &                 del, dist, az)

cc           1D ray tracing
c            call ttime(dist, src_dep(j), mod_nl, mod_v, mod_top,
c     &                 tmp_ttp(i, j), ain)
c            call ttime(dist, src_dep(j), mod_nl, vs, mod_top,
c     &                 tmp_tts(i, j), ain)

c	    Straight ray path:
c here you need to adopt for receiver depth, i.e. src_dep-sta_dep
            tmp_ttp(i,j)= sqrt(dist**2 + src_dep(j)**2)/mod_v(1)
            tmp_tts(i,j)= sqrt(dist**2 + src_dep(j)**2)/vs(1)
            ain= 180-(atan(dist/src_dep(j))*57.2958)

cc           Determine wave speed (k) at the hypocenter
c            do k=1,mod_nl
c               if (src_dep(j).le.mod_top(k)) goto 10	! break
c            enddo
c10          continue

c	    Straight ray path approach, only first entry in velocity
c           array taken:
            k= 2

c           Depth derivative
            tmp_zp(i,j) = cos((ain * pi)/180.0)/mod_v(k-1)
c           Epicentral derivatives
	    tmp_xp(i,j) = (sin((ain * pi)/180.0) *
     &               cos(((az - 90) * pi)/180.0))/mod_v(k-1)
	    tmp_yp(i,j) = (sin((ain * pi)/180.0) *
     &               cos((az * pi)/180.0))/mod_v(k-1)

c           Write to source-parameter file
            if (iunit .ne. 0)
     &         write(iunit,'(i9,2x,f9.4,2x,f9.4,2x,a7,2x,f9.4,
     &         2x,f9.4,2x,f9.4)')
     &         src_cusp(j), src_lat(j), src_lon(j), sta_lab(i),
     &         dist, az, ain

         enddo
      enddo

      if (iunit .ne. 0) close(iunit)	! Source-parameter file

      end !of subroutine partials


text:
!!# KINDS MODULE: <<KND_MoCshort>>
MODULE KND_MoCshort

!!## PURPOSE
!! Define the kind <KIND_MCS> for the Method of Characteristics
!! with ``short'' (intra-cell) interpolation transport method
!! (with acronym MCS) and make accessible some kinds defined elsewhere.

!!## MODULE DEPENDENCIES
USE KND_IntrinsicTypes !!((01-A-KND_IntrinsicTypes.f90))
USE KND_AngularFluxes  !!((02-A-KND_AngularFluxes.f90))
USE KND_ScalarFluxes   !!((02-A-KND_ScalarFluxes.f90))
USE KND_Currents       !!((02-A-KND_Currents.f90))
USE KND_Source         !!((02-A-KND_Source.f90))

!!## OWNER(S)
!! [waw] William A. Wieselquist, william.wieselquist AT gmail.com

!!## HISTORY
!!
!! + (01/01/2007) created [waw]

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PRIVATE

!!## KIND DECLARATIONS
!! The general purpose kind <<KIND_MCS>> for this package.
INTEGER,PARAMETER :: KIND_MCS = KIND_Rdp

!!## PUBLIC ACCESS LIST
PUBLIC :: KIND_AngularFlux
PUBLIC :: KIND_Current
PUBLIC :: KIND_MCS
PUBLIC :: KIND_ExtSource
PUBLIC :: KIND_ScalarFlux

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
C    %W% %G%
C****************************************************************
C
C     File: ld_geown.f
C
C     Purpose: Routine to load GE Owner data from raw data file
C
c     Return code:  n = 0 : Success
c                   n = 1 : Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_ge.f
C
C****************************************************************
      integer function ld_geown (xbuf, file, options, error)
      integer error, file, options(*)
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/owner_cm.inc'

      character word(20)*40, capital*40, tempc*40
      integer ftn_atoi, add_ptio, status, read_ge_file

      ld_geown = 0
      error = 0

      last = lastch0 (xbuf)
      if (xbuf(last:last) .eq. '/') then
        status = read_ge_file (file, xbuf(last:))
        if (status .eq. 0) go to 340
        last = lastch0 (xbuf)
      endif
      call uscan (xbuf(1:last), word, nwrd, '~',  ' ')

      do i = nwrd+1, 9
        word(i) = ' '
      enddo
c
c     Allow ownership number "0"
c
      npti1 = ftn_atoi(word(1))
      if (npti1 .eq. 0 .and. word(1) .ne. '0') then
        if (ichar (xbuf(1:1)) .ge. ichar ('a') .and.
     &      ichar (xbuf(1:1)) .le. ichar ('z')) ld_geown = 1
      else
        tempc = capital(word(3))
        num = add_ptio (npti1, tempc(1:3))
c
c       If num < 0, the owner exists
c
        if (num .lt. 0) then
          last = lastch (word(2))
          write (errbuf(1), 10000) npti1, word(3)(1:3),
     &       word(2)(1:last)
10000     format (' Owner GE record (',
     &        i6, 1x, a, 1x, a, ') could not be hashed.')
          call prterx ('W',1)
          error = 1
        else
          owner_name(num) = word(2)
        endif
      endif
      go to 900

  340 continue
      write (errbuf(1), 10090)
10090 format (' Premature E-O-F encountered reading Owner Data segment')
      write (errbuf(2), 10100) xbuf(1:80)
10100 format(11x, '(', a, ')')
      call prterx ('W',2)
      ld_geown = 1
      error = 1

  900 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      subroutine gglu( nsp, i0, ipi, ipf, lipotx, g0, tmatrx, g0t, gg)

      implicit real (a-h,o-z)
      implicit integer (i-n)
c  output
c    gg:  (nsp*lx**2, nsp*lx**2, 0:nphasx) submatrix spanning the entire
c          angular momentum basis for each unique potential

      include '../HEADERS/dim.h'
      include 'xparam.h'
      integer  i0 (0:nphx),  lipotx(0:nphx)

      parameter (one = 1, zero = 0)
      integer   ipiv(istatx)

c**** array of state kets at current energy
      common /stkets/ lrstat(4, istatx), istate

      complex   tmatrx(nspx, istatx)
c     big work matrices
      complex   g0( istatx, istatx), g0t( istatx, istatx)
      complex   g0s( istatx, nspx*(lx+1)**2)
c     return matrix containing info about each unique potential
      complex   gg(nspx*(lx+1)**2, nspx*(lx+1)**2, 0:nphasx)

      character*3  cerr
      character*13 trans

 400  format(i4)

c -------------------- LU gg
c     multiply T and G0 matrices together, construct g0t = 1 - G0*T
c     notice that the signs below for g0t ARE correct since 1 is the
c     unit matrix
c     since t is tri-diagonal, this product can be computed in n^2 time
c     also fill up some work matrices for use in eigenvalue and
c     determinant calculations and elsewhere
      do icol = 1,istate
         do irow = 1,istate
c         T diagonal contribution
            g0t(irow, icol) = - g0(irow, icol) * tmatrx(1, icol)
c         T off-diagonal contribution
            l1   = lrstat(2,icol)
            m1   = lrstat(3,icol)
            isp1 = lrstat(4,icol)
            m2 = m1+isp1
            if (nsp.eq.2 .and. m2.gt.-l1+1 .and. m2.lt.l1+2) then
               ist2 = icol + (-1)**isp1
               g0t(irow, icol) = g0t(irow, icol)
     1              - g0(irow, ist2) * tmatrx(nsp, icol)
            endif
         enddo

         g0t(icol, icol) = g0t(icol, icol) + one
      enddo

c --- invert matrix by LU decomposition
c     call cgetrf from lapack.  this performs an LU decomposition on
c     the matrix g0t = 1 - g0*T
      call cgetrf( istate, istate, g0t, istatx, ipiv, info )
      if (info.lt.0) then
         call wlog('    *** Error in cgetrf when computing G')
         write(cerr,400)abs(info)
         call wlog('        Argument #'//cerr//
     $        ' had an illegal value.')
      elseif (info.gt.0) then
         call wlog('    *** Error in cgetrf when computing G')
         write(cerr,400)info
         call wlog('        g0t('//cerr// ','//cerr//
     $        ') is exactly 0 -- '//
     $        'this matrix cannot be decomposed.')
      endif

c     now we want g_c = (g0t)^-1 * g0.  Rather than calculating
c     the inverse of g0t from the LU decomposition, we can compute
c     g0t^-1 * g0 directly by backsubstituting the columns of G0.
c     See sect. 2.3 in Numerical Recipes or LAPACK Users' Guide
c     sect. 2.3

c     third arg in number of output columns, istate for full
c     matrix, ipart(ik) for just the parts of the matrix needed
c     to contruct fine structure + DOS functions

      do ip=ipi, ipf
         ipart = nsp*(lipotx(ip)+1)**2
         do is1 = 1, istate
            do is2 = 1, ipart
               g0s(is1,is2) = g0(is1, is2 + i0(ip))
            enddo
         enddo

         trans = 'NotTransposed'
         call cgetrs(trans, istate, ipart, g0t, istatx,
     $        ipiv, g0s, istatx, info)
         if (info.lt.0) then
            call wlog('    *** Error in cgetrf')
            write(cerr,400) abs(info)
            call wlog('        Argument #'//cerr//
     $              ' had an invalid value.')
         endif

c **** at this point g0s contains the full MS ****

c  pack FMS matrix into an nsp*(lx+1)^2 x nsp*(lx+1)^2 matrix for each
c  ipot

         do is2=1,ipart
            do is1=1,ipart
               gg( is1, is2, ip) = g0s( is1+i0(ip), is2)
            enddo
         enddo
      enddo

      return
      end


program read_cond

  implicit none

  character(len=80) filename

  integer i, j, npts, interval, nmpi, nintervals, n_intervals_skip

  double precision x, y, cond_temp, average

  double precision, allocatable :: cond(:)

  print*,'enter filename'
  read(*,*) filename

  print*,'how many total points'
  read(*,*) npts

  print*,'stats reset interval'
  read(*,*) interval

  print*,'how many sample to throw out?'
  read(*,*) n_intervals_skip

  allocate(cond(npts))

  open(unit=999,file=trim(filename))
  do i=1,npts
     read(unit=999, FMT=*) cond(i)
  end do

  nintervals = npts / (interval)
  average = 0.

  do j=n_intervals_skip+1,nintervals
     cond_temp = cond(j*interval)
     print*,cond_temp
     average = average + cond_temp
  end do

  average = average / (nintervals-n_intervals_skip)
  print*,'overall average',average

end program read_cond


text:
! -------------------------------------------------------------
      subroutine surfcur_diag
!................................................................
c     Purpose:
c................................................................

      USE Vmeshes
      USE NumParams
      use OutCtrl, ONLY: w_jsurf
      USE LoopCtrl
      USE Vvacuum3
      USE Vprecal1
      USE Vprecal2
      USE Vprecal3
      USE Vprecal6
      USE Vprecal7
      USE Vsurface1
      USE Vsurface2, ONLY: xu, yu, xv
      USE Vsurface3, ONLY: yv, zu, zv
      USE Vsurface4, ONLY: xuu, yuu, zuu
      USE Vsurface5, ONLY: xuv, yuv, zuv
      USE Vsurface6, ONLY: xvv, yvv, zvv
      USE Vsurface7, ONLY: snx, sny, snz
      USE Vdiagno2
      USE Vdiagno3
      use Vsurfcurdiag
      implicit none
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      integer :: i, m, n, ku, kv
      real(rprec), dimension(:), allocatable ::  wx, wy, wz,
     1     curvat, cude, potu, potv, potuu, potuv, potvv, vx, vy, vz
      real(rprec) :: cw,cwm,sw,swm,cm,cn,v2,w2,vw
c ----------------------------------------------------------------------
      i=0
      if (.not.allocated(wx))
     1 allocate (wx(nuv), wy(nuv), wz(nuv), curvat(nuv), cude(nuv),
     1   potuu(nuv), potuv(nuv), potvv(nuv), vx(nuv), vy(nuv),
     2   vz(nuv), potu(nuv), potv(nuv), pote(nuv), stat = i)
      if (i .ne. 0) return
c
      vx= zero
      vy= zero
      vz= zero
      curvat= zero
      cude= zero
      wx= zero
      wy= zero
      wz= zero

      pote = zero
      potu = cut
      potv = cup
      potuu = zero
      potuv = zero
      potvv = zero

      do 30 n= 0,nf
         cn      = pi2*n
      do 30 m= 0,mf
         cm      = pi2*m
         i       = 0
      do 30  kv = 1,nv
      do 30  ku = 1,nu
         i      = i+1
         cw     = comu(ku,m)*conv(kv,n)-simu(ku,m)*sinv(kv,n)
         cwm    = comu(ku,m)*conv(kv,n)+simu(ku,m)*sinv(kv,n)
         sw     = simu(ku,m)*conv(kv,n)+comu(ku,m)*sinv(kv,n)
         swm    = simu(ku,m)*conv(kv,n)-comu(ku,m)*sinv(kv,n)
         pote(i)  = pote(i)+ eps(m)*eps(n)*(pot(m,n)*sw + pot(m,-n)*swm)
         potu(i)  = potu(i) + cm   *eps(n)*(pot(m,n)*cw + pot(m,-n)*cwm)
         potv(i)  = potv(i) + cn   *eps(m)*(pot(m,n)*cw - pot(m,-n)*cwm)
         potuu(i) = potuu(i)- cm*cm*eps(n)*(pot(m,n)*sw + pot(m,-n)*swm)
         potuv(i) = potuv(i)- cm      *cn* (pot(m,n)*sw - pot(m,-n)*swm)
         potvv(i) = potvv(i)- cn*cn*eps(m)*(pot(m,n)*sw + pot(m,-n)*swm)
   30 continue

c     curvature  of the surface current lines
c
      do i=1,nuv
         vx(i)  = xu(i)*potv(i)-xv(i)*potu(i)
         vy(i)  = yu(i)*potv(i)-yv(i)*potu(i)
         vz(i)  = zu(i)*potv(i)-zv(i)*potu(i)
         wx(i)  = xu(i)  *(potuv(i)*potv(i)-potvv(i)*potu(i))
     $          + xv(i)  *(potuv(i)*potu(i)-potuu(i)*potv(i))
     $          - potu(i)*(  xuv(i)*potv(i)-  xvv(i)*potu(i))
     $          - potv(i)*(  xuv(i)*potu(i)-  xuu(i)*potv(i))
c
         wy(i)  = yu(i)  *(potuv(i)*potv(i)-potvv(i)*potu(i))
     $          + yv(i)  *(potuv(i)*potu(i)-potuu(i)*potv(i))
     $          - potu(i)*(  yuv(i)*potv(i)-  yvv(i)*potu(i))
     $          - potv(i)*(  yuv(i)*potu(i)-  yuu(i)*potv(i))
c
         wz(i)  = zu(i)  *(potuv(i)*potv(i)-potvv(i)*potu(i))
     $          + zv(i)  *(potuv(i)*potu(i)-potuu(i)*potv(i))
     $          - potu(i)*(  zuv(i)*potv(i)-  zvv(i)*potu(i))
     $          - potv(i)*(  zuv(i)*potu(i)-  zuu(i)*potv(i))
      end do

      do i=1,nuvh
         v2     = vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
         w2     = wx(i)*wx(i)+wy(i)*wy(i)+wz(i)*wz(i)
         vw     = vx(i)*wx(i)+vy(i)*wy(i)+vz(i)*wz(i)
         curvat(i) = sqrt(w2-vw*vw/v2)/v2
         cude(i) = sqrt(v2/(snx(i)*snx(i)+sny(i)*sny(i)+snz(i)*snz(i)))
      end do

      curvra = maxval(curvat(1:nuvh))
      if (curvra .ne. zero) curvra = 1./curvra
      curvri = minval(curvat(1:nuvh))
      if (curvri .ne. zero ) curvri = 1./curvri
      cudema = maxval(cude(1:nuvh))
      cudemi = minval(cude(1:nuvh))
      avcude = (sum(cude(1:nu)) + 2.*sum(cude(1+nu:nuvh-nu))
     1       +  sum(cude(nuvh-nu+1:nuvh)))/float(nuv)

c.....Always Write this surface current info to output
      write (inesc, 20) 'J Surf max, min, ave = ',cudema, cudemi, avcude
      write (inesc, 20) 'Curv_R of J min, max = ',curvra, curvri

c.....Write extra surface current info to output if asked to do so
      if(w_jsurf>0 .or. w_jsurf==-1) then
         write (inesc, '(a,i4,a)') '---- potential phi(i) for i = 1,',
     1        nuv,' ----'
         write (inesc,*) (pote(i), i = 1, nuv)
         write (inesc, '(a)') '---- end  potential phi(i) ----'
      endif

      if(w_jsurf>1 .or. w_jsurf==-2) then
         write (inesc, '(a,i4,a)')
     1      '---- current jsurf(i) for i = 1,',nuvh,' ----'
         write (inesc,*) (cude(i), i = 1, nuvh)
         write (inesc, '(a)') '---- end  current jsurf(i) ----'
      endif

      if (iloop .le. 0) deallocate (wx, wy, wz, curvat, cude,
     1   potuu, potuv, potvv, vx, vy, vz, potu, potv)

 20   format(a,1p3e16.8)

      end subroutine surfcur_diag


!
!     ------------------------------------------------------------------
!              Y K F
!     ------------------------------------------------------------------
!
!               k
!       Stores Y (i, j; r) in the array YK
!
!
      SUBROUTINE YKF(I, J, K, REL)
      use param_C
      use radial_C
      use nel_C
!
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!...Translated by Pacific-Sierra Research 77to90  4.3E  00:46:39  11/16/01
!...Switches:
!      PARAMETER(NOD=220)
!      COMMON/PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,ID
!     :   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
!      COMMON/RADIAL/R(NOD),RR(NOD),R2(NOD),YK(NOD),YR(NOD),X(NOD)
!      POINTER(IQP,P(NOD,1)),(IQN,N(1)),(IQL,L(1)),(IQAZ,AZ(1)),
!     :       (IQMAX,MAX_(1))
!      COMMON/NEL/IQP,IQN,IQL,IQAZ,IQMAX,IQ(7)
      DIMENSION :: F(NOD)
      LOGICAL :: REL
!
      CALL ZK (I, J, K)
      A = EH**(K + 1)
      C = 2*K + 1
      A2 = A*A
      H90 = C*H3/D30
      A3 = A2*A*H90
      AI = H90/A
      AN = 114.D0*A*H90
      A34 = 34.D0*H90
      DO M = ND - 1, 2, -1
         F(M) = AN*YK(M+1) + A34*(YK(M)+A2*YK(M+2)) - YK(M-1)*AI - YK(M+3)*A3
      END DO
      F(1) = C*H3*(YK(1)+D4*A*YK(2)+A2*YK(3))
      DO M = ND - 1, 1, -1
         YK(M) = YK(M+2)*A2 + F(M)
      END DO
      IF (.NOT.REL) RETURN
      MM = MAX0(MAX_(I),MAX_(J))
      C = C*FINE
      DO M = 1, MM
         YK(M) = YK(M) + C*P(M,I)*P(M,J)
      END DO
      RETURN
      END SUBROUTINE YKF


module init_n_module

  use bl_types
  use bl_constants_module
  use ml_layout_module
  use multifab_physbc_module
  use define_bc_module
  use bc_module
  use bl_rng_module
  use bl_random_module
  use BoxLibRNGs

  use probin_common_module, only: prob_lo, prob_hi, prob_type, initial_variance, &
                                  perturb_width, smoothing_width
  use probin_reactdiff_module, only: nspecies, n_init_in, model_file_init, &
                                     cross_section, integer_populations, use_bl_rng

  implicit none

  private

  public :: init_n, init_n_model

  ! prob_type codes for LowMach:
  ! 0=thermodynamic equilibrium, n=n_init_in(1,1:nspecies)
  ! 1=gaussian spreading (order of accuracy testing)
  ! 2=gradient along y, n=n_init_in(1,1:nspecies) on bottom (y=0) and n_init(2,1:nspecies) on top (y=Ly)
  ! 3=1+sin^2(pi*x)*sin^2(pi*y)*sin^2(pi*z) test problem
  ! 4=vertical stripe (thirds of domain)

contains

  subroutine init_n(mla,n_init,dx,the_bc_tower)

    ! initialize rho_i and umac in the valid region
    ! we first initialize c_i in the valid region
    ! then enforce that sum(c_i)=1 by overwriting the final concentration,
    ! and then use the EOS to compute rho_i

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(inout) :: n_init(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(bc_tower ), intent(in   ) :: the_bc_tower

    ! local variables
    integer                        :: lo(mla%dim), hi(mla%dim)
    integer                        :: i, dm, n, nlevs, ng_n
    real(kind=dp_t), pointer       :: np(:,:,:,:)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "init_n")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_n = n_init(1)%ng

    ! looping over boxes
    do n=1,nlevs
       do i=1,nfabs(n_init(n))
          np => dataptr(n_init(n),i)
          lo = lwb(get_box(n_init(n),i))
          hi = upb(get_box(n_init(n),i))
          select case (dm)
          case (2)
             call init_n_2d(np(:,:,1,:),ng_n,lo,hi,dx(n,:))
          case (3)
             call init_n_3d(np(:,:,:,:),ng_n,lo,hi,dx(n,:))
          end select
       end do
    end do

    call destroy(bpt)

   do n=1,nlevs
      call multifab_fill_boundary(n_init(n))
      call multifab_physbc(n_init(n),1,scal_bc_comp,nspecies, &
                           the_bc_tower%bc_tower_array(n),dx_in=dx(n,:))
   end do

  end subroutine init_n

  subroutine init_n_2d(n_init,ng_n,lo,hi,dx)

    integer          :: lo(:), hi(:), ng_n
    real(kind=dp_t)  :: n_init(lo(1)-ng_n:,lo(2)-ng_n:,:)
    real(kind=dp_t)  :: dx(:)

    ! local varables
    integer         :: i,j,species
    real(kind=dp_t) :: x,y,r,cen(2),L(2)
    real(kind=dp_t) :: one_fraction_domain1,one_fraction_domain2,x1,x2,stripe_ratio,rad

    L(1:2) = prob_hi(1:2)-prob_lo(1:2) ! Domain length

    select case (prob_type)

    case(0)
       !============================================================
       ! Thermodynamic equilibrium
       !============================================================

       do j=lo(2),hi(2)
          do i=lo(1),hi(1)

             n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies)

          end do
       end do

    case(1)
       !=============================================================
       ! Initializing from a Gaussian
       !=============================================================

       cen(1:2) = 0.6d0*prob_lo(1:2) + 0.4d0*prob_hi(1:2)

       do j=lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
          do i=lo(1),hi(1)
             x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

             r = sqrt((x-cen(1))**2 + (y-cen(2))**2)

             n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies)*exp(-100.d0*r**2)

          end do
       end do

    case(2)
       !=========================================================
       ! Initializing with constant gradient
       !=========================================================

       do j=lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
          do i=lo(1),hi(1)
             x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

             ! linear gradient in rho
             n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies) + &
                  (n_init_in(2,1:nspecies) - n_init_in(1,1:nspecies))*(y-prob_lo(2))/L(2)

          end do
       end do

    case(3)
       !=========================================================
       ! 1+sin^2(pi*x)*sin^2(pi*y)*sin^2(pi*z) test problem
       !=========================================================

       do j=lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
          do i=lo(1),hi(1)
             x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

             n_init(i,j,1:nspecies) = 1.d0 + sin(M_PI*x)**2 * sin(M_PI*y)**2

          end do
       end do

    case(4)
       !=================================================================
       ! vertical stripe having width = perturb_width*dx(1)
       ! n_init = n_init_in(1,:) inside, n_init = n_init_in (2,:) outside
       ! can be discontinous or smooth depending on smoothing_width
       !=================================================================

       stripe_ratio = perturb_width*dx(1)/L(1)

       ! prob_lo(1) < one_fraction_domain2 < center < one_fraction_domain1 < prob_hi(1)
       one_fraction_domain1=0.5d0*(1.0d0-stripe_ratio)*prob_lo(1)+0.5d0*(1.0d0+stripe_ratio)*prob_hi(1)
       one_fraction_domain2=0.5d0*(1.0d0+stripe_ratio)*prob_lo(1)+0.5d0*(1.0d0-stripe_ratio)*prob_hi(1)

       do i=lo(1),hi(1)
          x1 =(prob_lo(1) + dx(1)*(dble(i)+0.5d0) - one_fraction_domain1)
          x2 =(prob_lo(1) + dx(1)*(dble(i)+0.5d0) - one_fraction_domain2)

          if (smoothing_width > 0.d0) then
             ! smooth interface
             do j=lo(2),hi(2)
                n_init(i,j,1:nspecies) = n_init_in(2,1:nspecies) + &
                     0.5d0*(n_init_in(2,1:nspecies)-n_init_in(1,1:nspecies)) * &
                     (tanh(x1/(smoothing_width*dx(1))) - tanh(x2/(smoothing_width*dx(1))))
             end do
          else
             ! discontinuous interface
             if (x2 < 0.d0 .or. x1 > 0.d0) then   ! outside the stripe
                do j=lo(2),hi(2)
                   n_init(i,j,1:nspecies) = n_init_in(2,1:nspecies)
                end do
             else                                 ! inside the stripe
                do j=lo(2),hi(2)
                   n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies)
                end do
             end if
          end if
       end do

    case(5)
       !=================================================================
       ! bubble having radius = 0.5*perturb_width*dx(1)
       ! n_init = n_init_in(1,:) inside, n_init = n_init_in (2,:) outside
       ! can be discontinous or smooth depending on smoothing_width
       !=================================================================

       rad = 0.5d0*perturb_width*dx(1)

       do j=lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+0.5d0)*dx(2) - 0.5d0*(prob_lo(2)+prob_hi(2))
          do i=lo(1),hi(1)
             x = prob_lo(1) + (dble(i)+0.5d0)*dx(1) - 0.5d0*(prob_lo(1)+prob_hi(1))

             r = sqrt(x**2 + y**2)

             if (smoothing_width .eq. 0) then
                ! discontinuous interface
                if (r .lt. rad) then
                   n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies)
                else
                   n_init(i,j,1:nspecies) = n_init_in(2,1:nspecies)
                end if
             else
                ! smooth interface
                n_init(i,j,1:nspecies) = n_init_in(1,1:nspecies) + &
                     (n_init_in(2,1:nspecies) - n_init_in(1,1:nspecies))* &
                     0.5d0*(1.d0 + tanh((r-rad)/(smoothing_width*dx(1))))
             end if

          end do
       end do

    case default

       call bl_error("init_n_2d: prob_type not supported")

    end select

    if(integer_populations) then ! Ensure that the initial number of molecules are integers
       if(initial_variance<0.0d0) then ! Distribute the particles on the box using a multinomial sampler

          ! To do this really correctly one should do this for the whole system, not box per box
          if (parallel_IOProcessor()) write(*,*) "Using multinomial initial distribution PER BOX, not domain"

          do species=1, nspecies
            call sample_integers(n_init(lo(1):hi(1),lo(2):hi(2),species), &
                     ncells=size(n_init(lo(1):hi(1),lo(2):hi(2),species)), &
                     dv=dx(1)*dx(2)*cross_section)
          end do

       else ! Make the number of molecules in each cell Poisson distributed with desired mean

          do j=lo(2),hi(2)
             do i=lo(1),hi(1)
                call round_to_integers(n_init(i,j,1:nspecies), dv=dx(1)*dx(2)*cross_section)
             end do
          end do

       end if
    end if

  end subroutine init_n_2d

  subroutine init_n_3d(n_init,ng_n,lo,hi,dx)

    integer          :: lo(:), hi(:), ng_n
    real(kind=dp_t)  :: n_init(lo(1)-ng_n:,lo(2)-ng_n:,lo(3)-ng_n:,:)
    real(kind=dp_t)  :: dx(:)

    ! local varables
    integer         :: i,j,k,species
    real(kind=dp_t) :: x,y,z,r,cen(3),L(3)
    real(kind=dp_t) :: one_fraction_domain1,one_fraction_domain2,x1,x2,stripe_ratio,rad

    L(1:3) = prob_hi(1:3)-prob_lo(1:3) ! Domain length

    select case (prob_type)

    case(0)
       !================================================================================
       ! Thermodynamic equilibrium
       !================================================================================

       do k=lo(3),hi(3)
          do j=lo(2),hi(2)
             do i=lo(1),hi(1)

                n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies)

             end do
          end do
       end do

    case(1)
       !================================================================================
       ! Initializing from a Gaussian
       !================================================================================
       cen(1:3) = 0.6d0*prob_lo(1:3) + 0.4d0*prob_hi(1:3)

       do k=lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+0.5d0)*dx(3)
          do j=lo(2),hi(2)
             y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
             do i=lo(1),hi(1)
                x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

                r = sqrt((x-cen(1))**2 + (y-cen(2))**2 + (z-cen(3))**2)

                n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies)*exp(-100.d0*r**2)

             end do
          end do
       end do

    case(2)
       !========================================================
       ! Initializing with constant gradient
       !========================================================

       do k=lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+0.5d0)*dx(3)
          do j=lo(2),hi(2)
             y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
             do i=lo(1),hi(1)
                x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

                n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies) + &
                     (n_init_in(2,1:nspecies) - n_init_in(1,1:nspecies))*(y-prob_lo(2))/L(2)

             end do
          end do
       end do

    case(3)
       !========================================================
       ! 1+sin^2(pi*x)*sin^2(pi*y)*sin^2(pi*z) test problem
       !========================================================

       do k=lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+0.5d0)*dx(3)
          do j=lo(2),hi(2)
             y = prob_lo(2) + (dble(j)+0.5d0)*dx(2)
             do i=lo(1),hi(1)
                x = prob_lo(1) + (dble(i)+0.5d0)*dx(1)

                n_init(i,j,k,1:nspecies) = 1.d0 + sin(M_PI*x)**2 * sin(M_PI*y)**2 * sin(M_PI*z)**2

             end do
          end do
       end do

    case(4)
       !=================================================================
       ! vertical stripe having width = perturb_width*dx(1)
       ! n_init = n_init_in(1,:) inside, n_init = n_init_in (2,:) outside
       ! can be discontinous or smooth depending on smoothing_width
       !=================================================================

       stripe_ratio = perturb_width*dx(1)/L(1)

       ! prob_lo(1) < one_fraction_domain2 < center < one_fraction_domain1 < prob_hi(1)
       one_fraction_domain1=0.5d0*(1.0d0-stripe_ratio)*prob_lo(1)+0.5d0*(1.0d0+stripe_ratio)*prob_hi(1)
       one_fraction_domain2=0.5d0*(1.0d0+stripe_ratio)*prob_lo(1)+0.5d0*(1.0d0-stripe_ratio)*prob_hi(1)

       do i=lo(1),hi(1)
          x1 =(prob_lo(1) + dx(1)*(dble(i)+0.5d0) - one_fraction_domain1)
          x2 =(prob_lo(1) + dx(1)*(dble(i)+0.5d0) - one_fraction_domain2)

          if (smoothing_width > 0.d0) then
             ! smooth interface
             do k=lo(3),hi(3)
             do j=lo(2),hi(2)
                n_init(i,j,k,1:nspecies) = n_init_in(2,1:nspecies) + &
                     0.5d0*(n_init_in(2,1:nspecies)-n_init_in(1,1:nspecies)) * &
                     (tanh(x1/(smoothing_width*dx(1))) - tanh(x2/(smoothing_width*dx(1))))
             end do
             end do
          else
             ! discontinuous interface
             if (x2 < 0.d0 .or. x1 > 0.d0) then   ! outside the stripe
                do k=lo(3),hi(3)
                do j=lo(2),hi(2)
                   n_init(i,j,k,1:nspecies) = n_init_in(2,1:nspecies)
                end do
                end do
             else                                 ! inside the stripe
                do k=lo(3),hi(3)
                do j=lo(2),hi(2)
                   n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies)
                end do
                end do
             end if
          end if
       end do

    case(5)
       !=================================================================
       ! bubble having radius = 0.5*perturb_width*dx(1)
       ! n_init = n_init_in(1,:) inside, n_init = n_init_in (2,:) outside
       ! can be discontinous or smooth depending on smoothing_width
       !=================================================================

       rad = 0.5d0*perturb_width*dx(1)

       do k=lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+0.5d0)*dx(3) - 0.5d0*(prob_lo(3)+prob_hi(3))
          do j=lo(2),hi(2)
             y = prob_lo(2) + (dble(j)+0.5d0)*dx(2) - 0.5d0*(prob_lo(2)+prob_hi(2))
             do i=lo(1),hi(1)
                x = prob_lo(1) + (dble(i)+0.5d0)*dx(1) - 0.5d0*(prob_lo(1)+prob_hi(1))

                r = sqrt(x**2 + y**2 + z**2)

                if (smoothing_width .eq. 0) then
                   ! discontinuous interface
                   if (r .lt. rad) then
                      n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies)
                   else
                      n_init(i,j,k,1:nspecies) = n_init_in(2,1:nspecies)
                   end if
                else
                   ! smooth interface
                   n_init(i,j,k,1:nspecies) = n_init_in(1,1:nspecies) + &
                        (n_init_in(2,1:nspecies) - n_init_in(1,1:nspecies))* &
                        0.5d0*(1.d0 + tanh((r-rad)/(smoothing_width*dx(1))))
                end if
             end do
          end do
       end do

    case default

       call bl_error("init_n_3d: prob_type not supported")

    end select

    if(integer_populations) then ! Ensure that the initial number of molecules are integers
       if(initial_variance<0.0d0) then ! Distribute the particles on the box using a multinomial sampler

          ! To do this really correctly one should do this for the whole system, not box per box
          if (parallel_IOProcessor()) write(*,*) "Using multinomial initial distribution PER BOX, not domain"

          do species=1, nspecies
            call sample_integers(n_init(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),species), &
                     ncells=size(n_init(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),species)), &
                     dv=dx(1)*dx(2)*dx(3)*cross_section)
          end do

       else ! Make the number of molecules in each cell Poisson distributed with desired mean

          do k=lo(3),hi(3)
             do j=lo(2),hi(2)
                do i=lo(1),hi(1)
                   call round_to_integers(n_init(i,j,k,1:nspecies), dv=dx(1)*dx(2)*dx(3)*cross_section)
                end do
             end do
          end do

       end if
    end if

  end subroutine init_n_3d

  subroutine round_to_integers(n, dv)
    real(dp_t), intent(inout) :: n(nspecies)
    real(dp_t), intent(in) :: dv ! cell volume

    integer :: comp, nparticles
    real(dp_t) :: random(nspecies)

    if(initial_variance>0.0d0) then
       do comp=1, nspecies
          ! Generate the initial fluctuations using a Poisson random number generator
          ! This assumes that the distribution of initial conditions is a product Poisson measure
          if (use_bl_rng) then
             call PoissonRNG(number=nparticles, mean=n(comp)*dv, engine=rng_eng_init%p)
          else
             call PoissonRNG(number=nparticles, mean=n(comp)*dv)
          end if
          n(comp) = nparticles/dv
       end do
    else ! Minimize fluctuations but ensure the number density is an integer
       ! It is important here to use smart rounding so that the average concentration is preserved in expectation
       ! Using nearest integer n = nint(n*dv)/dv will not work except for lots of molecules
       if (use_bl_rng) then
          call UniformRNGs(random, nspecies, engine=rng_eng_init%p)
       else
          call UniformRNGs(random, nspecies)
       end if
        ! Round to nearest integer
       n = floor( n*dv + random ) / dv
    end if

  end subroutine round_to_integers

  subroutine sample_integers(n, ncells, dv)
    integer, intent(in) :: ncells
    real(dp_t), intent(inout) :: n(ncells)
    real(dp_t), intent(in) :: dv ! cell volume

    integer, dimension(ncells) :: n_molecules ! Temporary array of integer number of molecules
    real(dp_t), dimension(ncells) :: p

    p=n/sum(n) ! Probability proportional to mean density
    p=p/sum(p) ! Make sure it sums to 1 to roundoff again

    if (use_bl_rng) then
       call MultinomialRNG(samples=n_molecules, n_samples=ncells, N=nint(sum(n)*dv), p=p, engine=rng_eng_init%p)
    else
       call MultinomialRNG(samples=n_molecules, n_samples=ncells, N=nint(sum(n)*dv), p=p)
    end if

    n = n_molecules/dv ! Convert back to number density

  end subroutine sample_integers

  subroutine init_n_model(mla,n_init,dx,the_bc_tower,input_array,comp) ! Read from a file

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(inout) :: n_init(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(bc_tower ), intent(in   ) :: the_bc_tower
    real(kind=dp_t), intent(in   ) :: input_array(:,:,:)
    integer        , intent(in   ) :: comp

    ! local variables
    integer                        :: lo(mla%dim), hi(mla%dim)
    integer                        :: i, dm, n, nlevs, ng_n
    real(kind=dp_t), pointer       :: np(:,:,:,:)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "init_n")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_n = n_init(1)%ng

    ! looping over boxes
    do n=1,nlevs
       do i=1,nfabs(n_init(n))
          np => dataptr(n_init(n),i)
          lo = lwb(get_box(n_init(n),i))
          hi = upb(get_box(n_init(n),i))
          select case (dm)
          case (2)
             call init_n_model_2d(np(:,:,1,:),ng_n,lo,hi,input_array(:,:,1),comp)
          case (3)
             call init_n_model_3d(np,ng_n,lo,hi,input_array,comp)
          end select
       end do
    end do

    call destroy(bpt)

    if (comp .eq. nspecies) then
       do n=1,nlevs
          call multifab_fill_boundary(n_init(n))
          call multifab_physbc(n_init(n),1,scal_bc_comp,nspecies, &
                               the_bc_tower%bc_tower_array(n),dx_in=dx(n,:))
       end do
    end if

  end subroutine init_n_model

  subroutine init_n_model_2d(n_init,ng_n,lo,hi,input_array,comp)

    integer         :: lo(:), hi(:), ng_n, comp
    real(kind=dp_t) :: n_init(lo(1)-ng_n:,lo(2)-ng_n:,:)
    real(kind=dp_t) :: input_array(0:,0:) ! This argument is replicated so it is the whole box, not just our patch!

    ! local varables
    integer         :: i,j

    do j=lo(2),hi(2)
    do i=lo(1),hi(1)

       if(model_file_init>0) then
          n_init(i,j,comp) = input_array(i,j)
       else
          n_init(i,j,comp) = input_array(j,i)
       end if

    end do
    end do

  end subroutine init_n_model_2d

  subroutine init_n_model_3d(n_init,ng_n,lo,hi,input_array,comp)

    integer         :: lo(:), hi(:), ng_n, comp
    real(kind=dp_t) :: n_init(lo(1)-ng_n:,lo(2)-ng_n:,lo(3)-ng_n:,:)
    real(kind=dp_t) :: input_array(0:,0:,0:) ! This argument is replicated so it is the whole box, not just our patch!

    ! local varables
    integer         :: i,j,k

    do k=lo(3),hi(3)
    do j=lo(2),hi(2)
    do i=lo(1),hi(1)

       if(model_file_init>0) then
          n_init(i,j,k,comp) = input_array(i,j,k)
       else
          n_init(i,j,k,comp) = input_array(k,j,i)
       end if

    end do
    end do
    end do

  end subroutine init_n_model_3d

end module init_n_module


text:
      SUBROUTINE HSURFP( NOUT, CVECT,
     &                 STRMAX, STRNOP, STRNON, STRCEN, STRRAD,
     &                 SHORTO, SAMPLE, PI,
     &                 ATMAX, ATNO, ATBRK, ATRES,
     &                 ATCHN, ATRNO, ATXYZ, ATVDW, ATBND, CUTSIZE)
      IMPLICIT NONE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Birkbeck College. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Birkbeck College are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 1996 Oliver Smart & Birkbeck College, All rights reserved    *
C * (c) 1997 Oliver Smart                    *
C *                                                                  *
C ********************************************************************
C
C Modification history:
C
C Date	Author		Modification
C 11/96 O.S. Smart      First version
C 12/97 O.S.S.		cut off lists introduced to speed calculations
C			change here is vble cutsize
C
C 5/11/96 Routine to calculate property graphs i.e. to find as function of
C coord along cvect the % of surface closest to an oxygen atom.
C Adapted from s/r H2DMAP which does a very similar thing - work out
C 2D property maps.

C passed variables *************

C output stream no. (usually set to 6) RETURN UNCHANGED
      INTEGER			NOUT

C The channel vector - already united
      DOUBLE PRECISION		CVECT(3)

C store for sphere centres and radii
C store the initial point as no 0 and subseqently points in the +ve
C bit. Points below the initial point are stored in -ve bit.
C maximum number of entries
      INTEGER                   STRMAX
C the no. of entries, +Ve entries -Ve entries
      INTEGER                   STRNOP, STRNON
C the centres and radii
      DOUBLE PRECISION          STRCEN(3,-STRMAX:STRMAX),
     &                          STRRAD(-STRMAX:STRMAX)

C 11/95 introduce new variable to reduce text output -
C if oshort =
C 0: full text output
C 1: no in run progress reported
C 2: ditto plus no graph type output - only leaving
C    Minimum radius & conductance calculations.
C 3: All text output other than input card mirroring
C    and error messages turned off
      INTEGER                   SHORTO

C sampling distance used
      DOUBLE PRECISION		SAMPLE

C pi
      DOUBLE PRECISION		PI

C maximum no. of atoms
      INTEGER                   ATMAX

C number of atoms read in from each file:
      INTEGER                   ATNO

C atom names found in Brookhaven file:
      CHARACTER*4               ATBRK(ATMAX)

C residue name in brook
      CHARACTER*3               ATRES(ATMAX)

C chain identifier in brookhaven pdb file
      CHARACTER*1               ATCHN(ATMAX)

C integer residue no.
      INTEGER                   ATRNO(ATMAX)

C co-ordinates
      DOUBLE PRECISION          ATXYZ( 3, ATMAX)

C vdw and bond radii of each atoms
      DOUBLE PRECISION          ATVDW(ATMAX), ATBND(ATMAX)

C cut off list control vble, -ve: no cutoff lists used in calculation
C otherwise this is the additional radius above the sphere radius
C to be used for cutoff list (see s/r holeen) needed here for calls to
C holeen
      DOUBLE PRECISION		CUTSIZE

C internal vbles ***************

C unit vectors used to define rotation
      DOUBLE PRECISION		PERPV1(3), PERPV2(3)

C dot product between cvect and perpv1
      DOUBLE PRECISION		CDOTP

C max// number of divisions to make in the angle
C (73 for 5 degrees)
      INTEGER			AMAX, ADIV
      PARAMETER(		AMAX= 140)

C loop counts angles, positions
      INTEGER			ACOUNT, SCOUNT

C value for angle, direction of spike, distance
C from sphere centre to any atoms vdw surface at that point,
C store for angle step
      DOUBLE PRECISION		ANGLE, DIRECT(3), DTOUCH,
     & 				ANSTEP

C need some dummy variables for s/r HOLEEN call
      INTEGER                   IAT1, IAT2, IAT3
      DOUBLE PRECISION          DAT2, DAT3
C also a test position and energy (to do spike option)
      DOUBLE PRECISION          NEWCEN(3), NEWENG

C y value for a line
      DOUBLE PRECISION		YLINE

C vbles introduced in writing hsurfp
C proportion of surface contacted by
C oxygen, nitrogen, oxygen or nitrogen, other.
C acid residue, base type residue, polar type residue,
C acid base or polar type residue, other (non-polar) type
C residue.
      INTEGER			PPO,
     &				PPN,
     &				PPOorN,
     &				PPOTH,
     &				RACID,
     &				RBASE,
     &				RAorB,
     &				RPOLAR,
     &				RABorP,
     &				RNONPO

C end of decs ******************

C do nothing if SHORTO set above 2
      IF (SHORTO.GT.2) GOTO 55555

C default number of angle sample - 73 gives 5 degree step
      ADIV = 73

C first thing to do is to define a vector at right angles
C to cvect from which the angle can be measured
      IF ( (ABS(CVECT(2)).LT.1E-06) .AND.
     &     (ABS(CVECT(3)).LT.1E-06)       ) THEN
        PERPV1(1) = 0.
        PERPV1(2) = 1.
        PERPV1(3) = 0.
      ELSE
        PERPV1(1) = 1.
        PERPV1(2) = 0.
        PERPV1(3) = 0.
      ENDIF
C now PERPV1 must be at right angles to CVECT - so find the component
C in common
      CALL DDOT( CVECT, PERPV1, CDOTP)
C as CVECT is a unit vector then can easily correct
      PERPV1(1) = PERPV1(1) - CDOTP*CVECT(1)
      PERPV1(2) = PERPV1(2) - CDOTP*CVECT(2)
      PERPV1(3) = PERPV1(3) - CDOTP*CVECT(3)
C make sure that this is a unit vector
      CALL DUVECT2( PERPV1)
C now get a vector which is normal both to CVECT and PERPV1
      CALL DCROSS( PERPV1, CVECT, PERPV2)

C tell user vectors used to define rotation
      WRITE( NOUT, '(A)')
     &' S/r HSURFP will output information to draw ',
     &'  Graph of surface properties vs coordinate along cvect',
     &'  A residue is regarded as',
     &'      Acid if it has id of: GLU or ASP;',
     &'                     id of: ARG or LYS;',
     &'                     id of: HIS, HID, HIE, GLN,',
     &'                                 THR, ASN or SER.',
     &' (This list is hard-coded in s/r HSURFP)'

C column titles
      WRITE( NOUT, '(20A10)')
     &    'cvect_coord', 'radius',
     &    'prop_O',     'prop_N',
     &    'prop_OorN',  'pro_other',
     &    'prop_Acid',  'prop_Base',
     &    'prop_AorB',  'pro_Polar',
     &    'pro_ABorP',  'prop_nonP'

C flush output stream
      IF (NOUT.EQ.6) CALL FLUSH6

C go through stored points from lowest to highest
      DO 11 SCOUNT = -STRNON, STRNOP
C work y value for this line and place in yline
        CALL DDOT( CVECT, STRCEN(1,SCOUNT), YLINE)

C zero the counts
        PPO = 0
        PPN = 0
        PPOorN = 0
        PPOTH = 0
        RACID = 0
        RBASE = 0
        RAorB = 0
        RPOLAR = 0
        RABorP = 0
        RNONPO = 0

C step in angle used
        ANSTEP = 2.*PI/FLOAT(ADIV-1)

C go through angle from -180 to 180
        DO 21 ACOUNT = 1, ADIV
C angle starts -pi goes to +pi
          ANGLE = 2.*PI*FLOAT(ACOUNT-1)/FLOAT(ADIV-1) - PI

C the direction of the point to be considered from the
C sphere centre is given by
          DIRECT(1) = COS(ANGLE)*PERPV1(1) +
     &                SIN(ANGLE)*PERPV2(1)
          DIRECT(2) = COS(ANGLE)*PERPV1(2) +
     &                SIN(ANGLE)*PERPV2(2)
          DIRECT(3) = COS(ANGLE)*PERPV1(3) +
     &                SIN(ANGLE)*PERPV2(3)
C load up point on surface of sphere at this point
          DTOUCH = STRRAD(SCOUNT)
          NEWCEN(1) = DIRECT(1)*DTOUCH + STRCEN(1,SCOUNT)
          NEWCEN(2) = DIRECT(2)*DTOUCH + STRCEN(2,SCOUNT)
          NEWCEN(3) = DIRECT(3)*DTOUCH + STRCEN(3,SCOUNT)
C do single hole calculation
C find energy of point newcen
C - energy is minus radius the largest sphere
C which can be centered at this point without vdW overlap
C put in fixed large cutsize
          CALL HOLEEN( NEWCEN, NEWENG, IAT1, IAT2, IAT3, DAT2, DAT3,
     &                         ATMAX, ATNO, ATXYZ, ATVDW, 5D0*CUTSIZE)
C what we want is info in iat1/iat2 which gives the closest
C and second closest atoms
C if closest atom is a hydrogen then disregard
          IF ( (ATBRK(IAT1)(1:1).EQ.'h') .OR.
     &         (ATBRK(IAT1)(1:1).EQ.'H')      ) IAT1 = IAT2

C store value of 1 if we have O, -1 if N and zero if other
          IF     ( (ATBRK(IAT1)(1:1).EQ.'o') .OR.
     &             (ATBRK(IAT1)(1:1).EQ.'O')      ) THEN
            PPO = PPO + 1
            PPOorN = PPOorN + 1

          ELSEIF ( (ATBRK(IAT1)(1:1).EQ.'n') .OR.
     &             (ATBRK(IAT1)(1:1).EQ.'N')      ) THEN
            PPN = PPN + 1
            PPOorN = PPOorN + 1

          ELSE
C other
            PPOTH = PPOTH + 1

          ENDIF

C look at residue polarity
          IF (     (ATRES(IAT1).EQ.'GLU') .OR.
     &             (ATRES(IAT1).EQ.'ASP')      ) THEN
            RACID  = RACID + 1
            RAorB  = RAorB + 1
            RABorP = RABorP + 1

          ELSEIF ( (ATRES(IAT1).EQ.'ARG') .OR.
     &             (ATRES(IAT1).EQ.'LYS')     ) THEN
            RBASE =  RBASE + 1
            RAorB  = RAorB + 1
            RABorP = RABorP + 1

          ELSEIF ( (ATRES(IAT1).EQ.'TYR') .OR.
     &             (ATRES(IAT1).EQ.'HIS') .OR.
     &             (ATRES(IAT1).EQ.'HID') .OR.
     &             (ATRES(IAT1).EQ.'HIE') .OR.
     &             (ATRES(IAT1).EQ.'GLN') .OR.
     &             (ATRES(IAT1).EQ.'THR') .OR.
     &             (ATRES(IAT1).EQ.'ASN') .OR.
     &             (ATRES(IAT1).EQ.'SER')     ) THEN
C not acid or base but polar
            RPOLAR = RPOLAR + 1
            RABorP = RABorP + 1
          ELSE
C other - non polar
            RNONPO = RNONPO + 1

          ENDIF

C end of going thru angle (acount loop)
21      CONTINUE

C write out information
        WRITE( NOUT, '(20F10.4)')
     &    YLINE,                     STRRAD(SCOUNT),
     &    FLOAT(PPO)/FLOAT(ADIV),    FLOAT(PPN)/FLOAT(ADIV),
     &    FLOAT(PPOorN)/FLOAT(ADIV), FLOAT(PPOTH)/FLOAT(ADIV),
     &    FLOAT(RACID)/FLOAT(ADIV),  FLOAT(RBASE)/FLOAT(ADIV),
     &    FLOAT(RAorB)/FLOAT(ADIV),  FLOAT(RPOLAR)/FLOAT(ADIV),
     &    FLOAT(RABorP)/FLOAT(ADIV), FLOAT(RNONPO)/FLOAT(ADIV)

C end of going thru stored points (scount loop)
11    CONTINUE

C return here
55555 CONTINUE
C flush output stream
      IF (NOUT.EQ.6) CALL FLUSH6

      RETURN
      END


text:
module enkf_obsmod
!$$$  module documentation block
!
! module: obsmod           module to read obs, obs metadata and ob priors
!                          (written out by GSI forward operator)
!
! prgmmr: whitaker         org: esrl/psd               date: 2009-02-23
!
! abstract: holds observation and observation prior data, and associated
! metadata. Includes subroutines for reading the data from files written
! out by GSI forward operator code, and screening out data with gross errors.
!
! Public Subroutines:
!   readobs: reads obs, obs priors and associated metadata from files
!    diag* files output from GSI forward operator, distributes data to
!    all tasks.
!   obsmod_cleanup: deallocate allocated arrays.
!
! Private Subroutines:
!   screenobs: screen out obs with large observation errors or
!    that fail background check.
!   channelstats: compute number of obs per satellite sensor/channel
!    (numobspersat) and mean observation error (oberrvarmean) - excluding
!    screened obs. Called by public subroutine
!    readobs.
!
! Public Variables (all defined by subroutine readobs):
!   nobs_conv (integer scalar):  number of "conventional" obs (from prepbufr file).
!   nobs_oz (integer scalar): number of sbuv ozone obs.
!   nobs_sat (integer scalar): number of satellite radiance obs.
!   nobstot (integer scalar): total number of obs (=nobs_conv+nobs_oz+nobs_sat)
!   jpch_rad: (integer scalar) total number of satellite sensors/channels
!    (imported from module radinfo).
!   npred: (integer scalar) total number of adaptive bias correction terms
!    (imported from module radinfo).
!   indxsat(nobs_sat): maps satellite radiance observation number
!    (1,2,3...,nobs_sat) to satellite sensor/channel number (1,2,3...jpch_rad).
!    jpch_rad is total number of satellite sensor/channels, defined by module
!    radinfo.
!   ob(nobstot): real array containing observation values.
!   obloclon(nobstot): real array of observation longitudes (degrees)
!   obloclat(nobstot): real array of observation longitudes (degrees)
!   obpress(nobstot): real array of observation pressures (hPa). For satellite
!    radiances, these are assumed to be the maximum of the weighting function
!    for that sensor/channel.
!   oblnp(nobstot): as above, but log(pressure).
!   obtime(nobstot): real array of observation time offsets (in hours,
!    relative to analysis time).
!   obtype(nobstot): character(20) array of observation types ('ps','  u','  v',
!    'oz','t','amsua_n15', etc.).
!   oberrvar(nobstot): real array of observation error variances (which have
!    been modified by GSI forward operator).
!   oberrvar_orig(nobstot):  real array of unmodified observation error
!    variance (as defined by bufr file or external error table).
!   ensmean_ob(nobstot): real array containing ensemble mean ob prior,
!    including bias correction.
!   ensmean_obnobc(nobstot): as above, but not includeing bias correction.
!   obfit_prior(nobstot): obs(nob)-ensmean_ob(nob) for prior.
!   obfit_post(nobstot): obs(nob)-ensmean_ob(nob) for posterior (not
!    defined until after analysis).
!   obsprd_prior(nobstot): real array of ensemble variance of
!     observation prior perturbations.
!   obsprd_post(nobstot): real array of ensemble variance of
!     observation posterior perturbations (not defined until after analysis,
!     and then only on root task).
!   numobspersat(jpch_rad):  number of screened obs for each
!     satellite sensor/channel.
!   oberrvarmean(jpch_rad):  mean observation error variance for
!     each satellite sensor/channel for all screened obs
!   biaspreds(npred+1, nobs_sat):  real array of bias predictors for
!     each satellite radiance ob (includes non-adaptive scan angle
!     bias correction term in biaspreds(1,1:nobs_sat)).
!   deltapredx(npred,jpch_rad): real array of bias coefficient increments
!     (initialized to zero, updated by analysis).
!   obloc(3,nobstot): real array of spherical cartesian coordinates
!     (x,y,z) of screened observation locations.
!   stattype(nobstot):  integer array containing prepbufr report type
!     (e.g. 120 for radiosonde temp) for "conventional" obs (nob <= nobs_conv)
!     and satellite channel number for satellite radiance obs (nob >
!     nobs_conv+nobs_oz). For ozone obs (nobs_conv<nob<nobs_sat),
!     set to 700 + k, where k is the  level index of the ozone retrieval.
!
! Modules Used: mpisetup, params, kinds, constants, mpi_readobs
!
! program history log:
!   2009-02-23  Initial version.
!   2011-06-20  Added the option of observation box for LETKF.
!   2015=07-25  Removed observation boxes for LETKF (use kdtree instead)
!   2016-03-21  Buchard  Add aerosol variables
!
! attributes:
!   language: f95
!
!$$$

use mpisetup
use kinds, only : r_kind, r_double, i_kind, r_single
use constants, only: zero, one, deg2rad, rad2deg, rd, cp, pi
use params, only: &
      datestring,datapath,sprd_tol,nanals,saterrfact, &
      lnsigcutoffnh, lnsigcutoffsh, lnsigcutofftr, corrlengthnh,&
      corrlengthtr, corrlengthsh, obtimelnh, obtimeltr, obtimelsh,&
      lnsigcutoffsatnh, lnsigcutoffsatsh, lnsigcutoffsattr,&
      varqc, huber, zhuberleft, zhuberright,&
      lnsigcutoffpsnh, lnsigcutoffpssh, lnsigcutoffpstr

use mpi_readobs, only:  mpi_getobs

implicit none
private
public :: readobs, obsmod_cleanup

real(r_single), public, allocatable, dimension(:) :: obsprd_prior, ensmean_obnobc,&
 ensmean_ob, ob, oberrvar, obloclon, obloclat, &
 obpress, obtime, oberrvar_orig,&
 oblnp, obfit_prior, prpgerr, oberrvarmean, probgrosserr, &
 lnsigl,corrlengthsq,obtimel
integer(i_kind), public, allocatable, dimension(:) :: numobspersat
! posterior stats computed in enkf_update
real(r_single), public, allocatable, dimension(:) :: obfit_post, obsprd_post
real(r_single), public, allocatable, dimension(:,:) :: biaspreds
real(r_kind), public, allocatable, dimension(:,:) :: deltapredx
! arrays passed to kdtree2 routines must be single.
real(r_single), public, allocatable, dimension(:,:) :: obloc
integer(i_kind), public, allocatable, dimension(:) :: stattype, indxsat
real(r_single), public, allocatable, dimension(:) :: biasprednorm,biasprednorminv
character(len=20), public, allocatable, dimension(:) :: obtype
integer(i_kind), public, allocatable, dimension(:) :: iused
integer(i_kind), public ::  nobs_sat, nobs_oz, nobs_conv, nobs_aer, nobstot

! for serial enkf, anal_ob is only used here and in loadbal. It is deallocated in loadbal.
! for letkf, anal_ob used on all tasks in letkf_update (bcast from root in loadbal), deallocated
! in letkf_update.
real(r_single), public, allocatable, dimension(:,:) :: anal_ob

contains

subroutine readobs()
! reads obs, obs priors and associated metadata from
! diag* files output from GSI forward operator, distributes data to
! all tasks.  Ob prior perturbations for each ensemble member
! are written to a temp file, since the entire array can be
! very large.
use radinfo, only: npred,nusis,nuchan,jpch_rad,iuse_rad,radinfo_read,predx,pg_rad
use convinfo, only: convinfo_read, init_convinfo, cvar_pg, nconvtype, ictype,&
                    ioctype
use ozinfo, only: init_oz, ozinfo_read, pg_oz, jpch_oz, nusis_oz, nulev
use covlocal, only: latval
integer nob,n,j,ierr
real(r_double) t1
real(r_single) tdiff,tdiffmax,deglat,radlat,radlon
! read in conv data info
call init_convinfo()
call convinfo_read()
! read in oz data info
call init_oz()
call ozinfo_read()
! read radiance bias correction info (standard out redirected
! specified unit number)
call radinfo_read()
!==> read in obs, obs metadata and ob priors.
!  (reading occurs on root, data broadcast to other tasks).
! temp files hxprime_ens_YYYYMMDDHH* are created containing
! ob prior perturbations for each ensemble member.
! make bias predictors unitless and O(1)
! so bias coefficents have same units as radiance obs.
! (by computing RMS values over many analyses)
if (nproc == 0) print*, 'npred  = ', npred
allocate(biasprednorm(npred),biasprednorminv(npred))
!biasprednorm(1) = 0.01_r_single   ! constant term
!biasprednorm(2) = 2.6e-2_r_single ! scan angle path
!biasprednorm(3) = 1.6e-2_r_single ! total column water
!biasprednorm(5) = 1.9e-2_r_single ! integrated weighted (by weighting fns) lapse rate
!biasprednorm(4) = zero   ! IWLR**2, don't use this predictor (too co-linear)?
!biasprednorm(4) = 1.1e-3_r_single
!biasprednorm(4) = zero   ! don't use this predictor (too co-linear)?
! what the heck, just scale them all by 0.01!
!biasprednorm = 0.01_r_single
biasprednorm=one
biasprednorminv=zero
do n=1,npred
   if (nproc == 0) print *,n,'biasprednorm = ',biasprednorm(n)
   if (biasprednorm(n) > 1.e-7_r_single) biasprednorminv(n)=one/biasprednorm(n)
enddo
! scale bias coefficients.
do j=1,jpch_rad
   predx(:,j) = biasprednorm(:)*predx(:,j)
enddo
! allocate array for bias correction increments, initialize to zero.
allocate(deltapredx(npred,jpch_rad))
deltapredx = zero
t1 = mpi_wtime()
call mpi_getobs(datapath, datestring, nobs_conv, nobs_oz, nobs_sat, nobs_aer, nobstot, &
                obsprd_prior, ensmean_obnobc, ensmean_ob, ob, &
                oberrvar, obloclon, obloclat, obpress, &
                obtime, oberrvar_orig, stattype, obtype, biaspreds,&
                anal_ob,indxsat,nanals)
tdiff = mpi_wtime()-t1
call mpi_reduce(tdiff,tdiffmax,1,mpi_real4,mpi_max,0,mpi_comm_world,ierr)
if (nproc == 0) then
 print *,'max time in mpireadobs  = ',tdiffmax
 print *,'total number of obs ',nobstot
 allocate(iused(nobstot)); iused=0
endif
allocate(obfit_prior(nobstot))
! screen out some obs by setting ob error to a very large number
! set obfit_prior
call screenobs()

allocate(probgrosserr(nobstot),prpgerr(nobstot))
! initialize prob of gross error to 0.0 (will be reset by analysis if varqc is true)
probgrosserr = zero
if (varqc .and. .not. huber) then
   ! for flat-tail VarQC, read in a-prior prob of gross error.
   prpgerr = zero ! initialize to zero
   do nob=1,nobstot
      if (nob <= nobs_conv) then
         ! search for matching record in convinfo file.
         ! if match found, set prob. of gross error to nonzero value given in
         ! file.
         do j=1,nconvtype
            if (ictype(j) == stattype(nob) .and. &
                ioctype(j) == obtype(nob)(1:3)) then
                prpgerr(nob) = cvar_pg(j)
                exit
            end if
         enddo
      else if (nob <= nobs_conv+nobs_oz) then
         do j=1,jpch_oz
            if (stattype(nob)-700 == nulev(j) .and. obtype(nob) == nusis_oz(j)) then
                prpgerr(nob) = pg_oz(j)
                exit
            end if
         enddo
      else
         prpgerr(nob) = pg_rad(indxsat(nob-(nobs_conv+nobs_oz)))
      end if
      !if (nproc == 0) print *,nob,obtype(nob)(1:3),prpgerr(nob)
   enddo
endif
! compute number of usuable obs, average ob error for each satellite sensor/channel.
if (nobs_sat > 0) then
  do nob=1,nobs_sat
     do n=2,npred+1
       biaspreds(n,nob)=biaspreds(n,nob)* biasprednorminv(n-1)
     end do
  end do
  call channelstats()
end if

! calculate locations of obs that passed initial screening in cartesian coords.
allocate(obloc(3,nobstot))
allocate(oblnp(nobstot)) ! log(p) at ob locations.
allocate(corrlengthsq(nobstot),lnsigl(nobstot),obtimel(nobstot))
lnsigl=1.e10
do nob=1,nobstot
   oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units
   if (obloclon(nob) < zero) obloclon(nob) = obloclon(nob) + 360._r_single
   radlon=deg2rad*obloclon(nob)
   radlat=deg2rad*obloclat(nob)
! cartesian coordinates of 'good' obs.
   obloc(1,nob) = cos(radlat)*cos(radlon)
   obloc(2,nob) = cos(radlat)*sin(radlon)
   obloc(3,nob) = sin(radlat)
   deglat = obloclat(nob)
!  get limits on corrlength,lnsig,and obtime
   if (nob > nobs_conv+nobs_oz) then
      if (obtype(nob)(1:3) == 'aod') then
         lnsigl(nob)=-1.e5_r_single!latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh)
      else
         lnsigl(nob) = latval(deglat,lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh)
      endif
   else if (obtype(nob)(1:3) == ' ps') then
      lnsigl(nob) = latval(deglat,lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh)
   else
      lnsigl(nob)=latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh)
   end if
   corrlengthsq(nob)=latval(deglat,corrlengthnh,corrlengthtr,corrlengthsh)**2
   obtimel(nob)=latval(deglat,obtimelnh,obtimeltr,obtimelsh)
end do

! these allocated here, but not computed till after the state
! update in enkf_update.
allocate(obfit_post(nobstot))
allocate(obsprd_post(nobstot))
obsprd_post = zero
end subroutine readobs

subroutine screenobs()
! screen out obs with large observation errors or
! that fail background check.  For screened obs oberrvar is set to 1.e31_r_single
!use radbias, only: apply_biascorr
use radinfo, only: iuse_rad,nuchan,nusis,jpch_rad
real(r_single) fail,failm
integer nn,nob
fail=1.e31_r_single
failm=1.e30_r_single
! apply bias correction here just for debugging purposes.
!call apply_biascorr()
!==> pre-process obs, obs metadata.
do nob=1,nobstot
  if (nob > nobs_conv+nobs_oz) oberrvar(nob) = saterrfact*oberrvar(nob)
  ! empirical adjustment of obs errors for Huber norm from ECMWF RD tech memo
  if (varqc) oberrvar(nob) = oberrvar(nob)*(min(one,0.5_r_single+0.125_r_single*(zhuberleft+zhuberright)))**2

  obfit_prior(nob) = ob(nob)-ensmean_ob(nob)

  ! gross error check.
  if (obsprd_prior(nob) > 1.e9) then
    oberrvar(nob)= fail
    if (obtype(nob)(1:3) == '  u') oberrvar(nob)=fail
    if (obtype(nob)(1:3) == '  v' .and. oberrvar(nob-1) >= failm) oberrvar(nob-1)=fail
    if(nproc == 0)print *, ' ob rejected due to obsprd_prior > 1.e9 ',obsprd_prior(nob)
    cycle
  end if

  if (obtype(nob) == ' ps' .and. (stattype(nob) == 111 .or. stattype(nob) == 112)) then
      ! tcvitals obs get a free pass for background check.
      if (nproc == 0) then
      print *,'tcvitals ps ob found (lat,lon,time,ob,O-F,sprd,oberr_std):'
      write(6,9101) stattype(nob),obloclat(nob),obloclon(nob),obtime(nob),ob(nob),obfit_prior(nob), &
                    sqrt(obsprd_prior(nob)),sqrt(oberrvar(nob))
      endif
  else
     ! background checks.
     if (abs(obfit_prior(nob)) > sprd_tol*sqrt(obsprd_prior(nob)+oberrvar(nob))) then
       !print *,obtype(nob)(1:3),stattype(nob),ob(nob),ensmean_ob(nob),obsprd_prior(nob),oberrvar(nob)
       if(nproc == 0)print *, 'reject ob due to obsprd_prior > sprd_tol ',&
          obtype(nob)(1:3),obfit_prior(nob),sprd_tol*sqrt(obsprd_prior(nob)+oberrvar(nob))
       oberrvar(nob)= fail
       if (obtype(nob)(1:3) == '  u') oberrvar(nob)=fail
       if (obtype(nob)(1:3) == '  v' .and. oberrvar(nob-1) >= failm) oberrvar(nob-1)=fail
       cycle
     end if
  end if

enddo
if (nproc == 0) then
   nn = 0
   do nob=1,nobstot
     if(oberrvar(nob) >= failm)nn = nn + 1
   end do
   print *,nobstot-nn,'obs kept'
   print *,nn,'total obs rejected'
end if ! nproc=0
9101 format(i3,7(1x,f7.2))

end subroutine screenobs

subroutine channelstats
use radinfo, only: npred,nusis,nuchan,jpch_rad
implicit none
integer(i_kind) nob,nob2,i
! count number of obs per channel/sensor.
allocate(numobspersat(jpch_rad))
allocate(oberrvarmean(jpch_rad))
numobspersat = 0
oberrvarmean = zero
do nob=1,nobs_sat
   nob2=nob+nobs_conv+nobs_oz
   i=indxsat(nob)
   numobspersat(i) = numobspersat(i) + 1
   oberrvarmean(i) = oberrvarmean(i) + oberrvar(nob2)
enddo
!if (nproc == 0) then
!   print *,'numobspersat:',minval(numobspersat),maxval(numobspersat)
!   print *,'oberrvarmean:',minval(oberrvarmean),maxval(oberrvarmean)
!end if
! average ob error for each channel.
do i=1,jpch_rad
   if (numobspersat(i) > 0) then
      oberrvarmean(i) = oberrvarmean(i)/real(numobspersat(i),r_single)
   else
      oberrvarmean(i) = 9.9e31_r_single
   end if
enddo

end subroutine channelstats

subroutine obsmod_cleanup()
! deallocate module-level allocatable arrays
if (allocated(obsprd_prior)) deallocate(obsprd_prior)
if (allocated(obfit_prior)) deallocate(obfit_prior)
if (allocated(obsprd_post)) deallocate(obsprd_post)
if (allocated(obfit_post)) deallocate(obfit_post)
if (allocated(ensmean_ob)) deallocate(ensmean_ob)
if (allocated(ensmean_obnobc)) deallocate(ensmean_obnobc)
if (allocated(ob)) deallocate(ob)
if (allocated(oberrvar)) deallocate(oberrvar)
if (allocated(oberrvar_orig)) deallocate(oberrvar_orig)
if (allocated(oberrvarmean)) deallocate(oberrvarmean)
if (allocated(obloclon)) deallocate(obloclon)
if (allocated(obloclat)) deallocate(obloclat)
if (allocated(obpress)) deallocate(obpress)
if (allocated(obtime)) deallocate(obtime)
if (allocated(oblnp)) deallocate(oblnp)
if (allocated(numobspersat)) deallocate(numobspersat)
if (allocated(obloc)) deallocate(obloc)
if (allocated(biaspreds)) deallocate(biaspreds)
if (allocated(deltapredx)) deallocate(deltapredx)
if (allocated(indxsat)) deallocate(indxsat)
if (allocated(obtype)) deallocate(obtype)
if (allocated(probgrosserr)) deallocate(probgrosserr)
if (allocated(prpgerr)) deallocate(prpgerr)
if (allocated(biasprednorm)) deallocate(biasprednorm)
if (allocated(biasprednorminv)) deallocate(biasprednorminv)
if (allocated(anal_ob)) deallocate(anal_ob)
if (nproc==0) then
   if (allocated(iused)) deallocate(iused)
endif
end subroutine obsmod_cleanup

end module enkf_obsmod

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
MODULE MODULE_LANDSCAPES

USE MODULE_SORT

CONTAINS

!SUBROUTINE BOTTOM_UP_1D(x, y, T_start, T_ind, T_pn, n_nodes, n_edges)
!
!    IMPLICIT NONE
!
!    TYPE array_pointer
!        INTEGER,DIMENSION(:),POINTER::p1
!    END TYPE array_pointer
!
!
!    INTEGER,INTENT(IN)::n_nodes,n_edges
!    INTEGER,DIMENSION(n_edges),INTENT(IN)::T_ind
!    INTEGER,DIMENSION(n_nodes+1),INTENT(IN)::T_start
!    DOUBLE PRECISION, DIMENSION(n_nodes),INTENT(IN)::T_pn
!    DOUBLE PRECISION, DIMENSION(n_nodes),INTENT(OUT)::x,y
!
!    INTEGER, DIMENSION(n_nodes):: order, belongs_to
!    LOGICAL, DIMENSION(n_nodes):: checked
!
!    TYPE(array_pointer),DIMENSION(:),POINTER::basin
!    INTEGER, DIMENSION(:), ALLOCATABLE:: basin_lmin, basin_lmax, basin_coin
!    INTEGER::n_basins
!
!    INTEGER, DIMENSION
!
!    CALL ARGSORT_REVERSE(T_pn, order)
!
!    x(:)=0.0d0
!    y(:)=T_pn(:)
!
!    checked(:)=.FALSE.
!    belongs_to=0
!    n_basins=0
!
!    ALLOCATE(basin_lmin(n_basins), basin_lmax(n_basins), basin_coin(n_basins))
!    ALLOCATE()
!
!    DO ii=1, n_nodes
!
!        ind_node = order(ii)
!
!        new_basin = .True.
!
!        DO jj=T_start(ind_node)+1, T_start(ind_node+1)
!
!            IF (checked(jj).eqv..True.) THEN
!
!                new_basin = .False.
!
!            END IF
!
!        END DO
!
!        IF (new_basin.eqv..True.) THEN
!
!            n_basins = n_basins+1
!            belongs_to(ind_node)=n_basins
!
!
!
!
!END SUBROUTINE BOTTOM_UP_1D

END MODULE MODULE_LANDSCAPES

!! f2py --f90flags='-fast' -c -m pyn_fort_net pyn_fort_net.f90
!!f2py --f90flags='-fast' -c -m pyn_fort_net pyn_fort_net.f90 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lmkl_def -lpthread

max_stars_repo_name: uibcdf/OpenKinNet

max_issues_repo_name: uibcdf/OpenKinNet

max_forks_repo_name: uibcdf/OpenKinNet


text:
c=======================================================================
c
c     subroutine SKTHGENOPT
c
c     Computing the function and gradient of
c
c     Theta(dual) = 0.5*[wopt'*Q*wopt] + wopt'*plin + cste
c
c-----------------------------------------------------------------------
      SUBROUTINE skthgenopt ( n, Q, kapplin, rho, omega, omega0,
     &                        ndual, ydual, delta, mu, dist,
     &                        dwork, funct, grad )
c-----------------------------------------------------------------------
c
c     INPUT
c       n       : number of values                               integer
c       Q       : semi definite matrix (n*n)                      double
c       kapplin : vector(n) linear part of the objective funtcion double
c       rho     : mean returns (n)                                double
c       omega   : vector(n) wights of the portfolio               double
c       omega0  : vector(n) second order portfolio                double
c       ndual   : number of variables                            integer
c       ydual   : dual solution lambda (ndual)                    double
c       delta   : maximal distance (with respect to the L2 norm)
c                 between omega0 and the optimal portfolio        double
c       mu      : value of the target return                      double
c       dist    : distance (with respect to the L2 norm)
c                 between omega0 and omega                        double
c
c     WORKSPACE
c       dwork   : n                                               double
c
c     OUTPUT
c       funct   : dual function value                             double
c       grad    : (2)-array gradient value                        double
c
c     CALL
c       OVTMCV, XV, SVVX
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, ndual
      DOUBLE PRECISION funct, cste, delta, mu, dist, Q(*), kapplin(*),
     &                 omega(*), rho(*), ydual(*), omega0(*), grad(*)
c
c     workspaces
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER pdw
      DOUBLE PRECISION var, scal, lin, perf, NEGONE
      PARAMETER (NEGONE = -1.E0)
c
c     external functions
      EXTERNAL OVTMCV, XV, SVVX
c
c-----------------------------------------------------------------------
c
c     pointers for double precision work space : dwork
c     ------------------------------------------------
c     Double workspace
      pdw = 1
c
c     Total size of iwork array = n
c
c----------------------------------------------------------------
c
      CALL OVTMCV(n, Q, omega, var)  ! w'*Q*w
      CALL XV(n, omega, kapplin, lin)
      funct = 0.5*var + lin + ydual(1)*mu + ydual(2)*dist
      funct = -funct
      CALL XV(n, rho, omega, perf)
      grad(1) = -(mu - perf)
      CALL SVVX (n, omega, omega0, NEGONE, dwork(pdw))
      CALL XV(n, dwork(pdw), dwork(pdw), scal)
      grad(2)= -(scal - delta*delta)
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


! { dg-do run }
!
! PR fortran/18918
!
! this_image(coarray) run test,
! expecially for num_images > 1
!
! Tested are values up to num_images == 8,
! higher values are OK, but not tested for
!
implicit none
integer :: a(1)[2:2, 3:4, 7:*]
integer :: b(:)[:, :,:]
allocatable :: b
integer :: i

if (this_image(A, dim=1) /= 2) STOP 1
i = 1
if (this_image(A, dim=i) /= 2) STOP 2

select case (this_image())
  case (1)
    if (this_image(A, dim=2) /= 3) STOP 3
    if (this_image(A, dim=3) /= 7) STOP 4
    i = 2
    if (this_image(A, dim=i) /= 3) STOP 5
    i = 3
    if (this_image(A, dim=i) /= 7) STOP 6
    if (any (this_image(A) /= [2,3,7])) STOP 7

  case (2)
    if (this_image(A, dim=2) /= 4) STOP 8
    if (this_image(A, dim=3) /= 7) STOP 9
    i = 2
    if (this_image(A, dim=i) /= 4) STOP 10
    i = 3
    if (this_image(A, dim=i) /= 7) STOP 11
    if (any (this_image(A) /= [2,4,7])) STOP 12

  case (3)
    if (this_image(A, dim=2) /= 3) STOP 13
    if (this_image(A, dim=3) /= 8) STOP 14
    i = 2
    if (this_image(A, dim=i) /= 3) STOP 15
    i = 3
    if (this_image(A, dim=i) /= 8) STOP 16
    if (any (this_image(A) /= [2,3,8])) STOP 17

  case (4)
    if (this_image(A, dim=2) /= 4) STOP 18
    if (this_image(A, dim=3) /= 8) STOP 19
    i = 2
    if (this_image(A, dim=i) /= 4) STOP 20
    i = 3
    if (this_image(A, dim=i) /= 8) STOP 21
    if (any (this_image(A) /= [2,4,8])) STOP 22

  case (5)
    if (this_image(A, dim=2) /= 3) STOP 23
    if (this_image(A, dim=3) /= 9) STOP 24
    i = 2
    if (this_image(A, dim=i) /= 3) STOP 25
    i = 3
    if (this_image(A, dim=i) /= 9) STOP 26
    if (any (this_image(A) /= [2,3,9])) STOP 27

  case (6)
    if (this_image(A, dim=2) /= 4) STOP 28
    if (this_image(A, dim=3) /= 9) STOP 29
    i = 2
    if (this_image(A, dim=i) /= 4) STOP 30
    i = 3
    if (this_image(A, dim=i) /= 9) STOP 31
    if (any (this_image(A) /= [2,4,9])) STOP 32

  case (7)
    if (this_image(A, dim=2) /= 3) STOP 33
    if (this_image(A, dim=3) /= 10) STOP 34
    i = 2
    if (this_image(A, dim=i) /= 3) STOP 35
    i = 3
    if (this_image(A, dim=i) /= 10) STOP 36
    if (any (this_image(A) /= [2,3,10])) STOP 37

  case (8)
    if (this_image(A, dim=2) /= 4) STOP 38
    if (this_image(A, dim=3) /= 10) STOP 39
    i = 2
    if (this_image(A, dim=i) /= 4) STOP 40
    i = 3
    if (this_image(A, dim=i) /= 10) STOP 41
    if (any (this_image(A) /= [2,4,10])) STOP 42
end select

allocate (b(3)[-1:0,2:4,*])

select case (this_image())
  case (1)
    if (this_image(B, dim=1) /= -1) STOP 43
    if (this_image(B, dim=2) /= 2) STOP 44
    if (this_image(B, dim=3) /= 1) STOP 45
    i = 1
    if (this_image(B, dim=i) /= -1) STOP 46
    i = 2
    if (this_image(B, dim=i) /= 2) STOP 47
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 48
    if (any (this_image(B) /= [-1,2,1])) STOP 49

  case (2)
    if (this_image(B, dim=1) /= 0) STOP 50
    if (this_image(B, dim=2) /= 2) STOP 51
    if (this_image(B, dim=3) /= 1) STOP 52
    i = 1
    if (this_image(B, dim=i) /= 0) STOP 53
    i = 2
    if (this_image(B, dim=i) /= 2) STOP 54
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 55
    if (any (this_image(B) /= [0,2,1])) STOP 56

  case (3)
    if (this_image(B, dim=1) /= -1) STOP 57
    if (this_image(B, dim=2) /= 3) STOP 58
    if (this_image(B, dim=3) /= 1) STOP 59
    i = 1
    if (this_image(B, dim=i) /= -1) STOP 60
    i = 2
    if (this_image(B, dim=i) /= 3) STOP 61
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 62
    if (any (this_image(B) /= [-1,3,1])) STOP 63

  case (4)
    if (this_image(B, dim=1) /= 0) STOP 64
    if (this_image(B, dim=2) /= 3) STOP 65
    if (this_image(B, dim=3) /= 1) STOP 66
    i = 1
    if (this_image(B, dim=i) /= 0) STOP 67
    i = 2
    if (this_image(B, dim=i) /= 3) STOP 68
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 69
    if (any (this_image(B) /= [0,3,1])) STOP 70

  case (5)
    if (this_image(B, dim=1) /= -1) STOP 71
    if (this_image(B, dim=2) /= 4) STOP 72
    if (this_image(B, dim=3) /= 1) STOP 73
    i = 1
    if (this_image(B, dim=i) /= -1) STOP 74
    i = 2
    if (this_image(B, dim=i) /= 4) STOP 75
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 76
    if (any (this_image(B) /= [-1,4,1])) STOP 77

  case (6)
    if (this_image(B, dim=1) /= 0) STOP 78
    if (this_image(B, dim=2) /= 4) STOP 79
    if (this_image(B, dim=3) /= 1) STOP 80
    i = 1
    if (this_image(B, dim=i) /= 0) STOP 81
    i = 2
    if (this_image(B, dim=i) /= 4) STOP 82
    i = 3
    if (this_image(B, dim=i) /= 1) STOP 83
    if (any (this_image(B) /= [0,4,1])) STOP 84

  case (7)
    if (this_image(B, dim=1) /= -1) STOP 85
    if (this_image(B, dim=2) /= 2) STOP 86
    if (this_image(B, dim=3) /= 2) STOP 87
    i = 1
    if (this_image(B, dim=i) /= -1) STOP 88
    i = 2
    if (this_image(B, dim=i) /= 2) STOP 89
    i = 3
    if (this_image(B, dim=i) /= 2) STOP 90
    if (any (this_image(B) /= [-1,2,2])) STOP 91

  case (8)
    if (this_image(B, dim=1) /= 0) STOP 92
    if (this_image(B, dim=2) /= 2) STOP 93
    if (this_image(B, dim=3) /= 2) STOP 94
    i = 1
    if (this_image(B, dim=i) /= 0) STOP 95
    i = 2
    if (this_image(B, dim=i) /= 2) STOP 96
    i = 3
    if (this_image(B, dim=i) /= 2) STOP 97
    if (any (this_image(B) /= [0,2,2])) STOP 98
end select

end