text:
MODULE io2_module

  use parallel_module
  use rgrid_module, only: Igrid, dV, Hgrid
  use array_bound_module
  use wf_module
  use watch_module

  implicit none

  PRIVATE
  PUBLIC :: read_io2, read_data_io2

  integer :: nprocs_new, nmax
  integer,allocatable :: nmap(:),iomap(:,:,:)
  integer :: MB_0_IO, MB_1_IO

CONTAINS

  SUBROUTINE read_io2( rank )
    implicit none
    integer,intent(IN) :: rank
    integer,parameter  :: unit0=10
    integer :: i,j,ierr,idummy
    integer :: ML_tmp,ML1_tmp,ML2_tmp,ML3_tmp
    integer :: MB_tmp,MB1_tmp,MB2_tmp
    character(10) :: cbuf
!---
    if ( rank == 0 ) then
       open(unit0,file='input_wfiodir',status='old')
       read(unit0,*) nprocs_new
       allocate( nmap(nprocs_new) ) ; nmap=0
       do i=1,nprocs_new
          read(unit0,*) nmap(i)
          do j=1,nmap(i)
             read(unit0,*)
          end do
       end do
       nmax=maxval( nmap )
       allocate( iomap(0:6,nmax,nprocs_new) ) ; iomap=0
       rewind unit0
       read(unit0,*) nprocs_new
       do i=1,nprocs_new
          read(unit0,*) nmap(i)
          do j=1,nmap(i)
             read(unit0,*) iomap(0,j,i)
          end do
       end do
       do i=1,nprocs_new
          read(unit0,*)
          do j=1,nmap(i)
             read(unit0,*) cbuf,idummy,iomap(1:6,j,i)
          end do
       end do
       close(unit0)
    end if
!---
    call mpi_bcast(nmax,1,mpi_integer,0,mpi_comm_world,ierr)
    call mpi_bcast(nprocs_new,1,mpi_integer,0,mpi_comm_world,ierr)
    if ( rank /= 0 ) then
       allocate( nmap(nprocs_new) ) ; nmap=0
       allocate( iomap(0:6,nmax,nprocs_new) ) ; iomap=0
    end if
    call mpi_bcast(nmap,nprocs_new,mpi_integer,0,mpi_comm_world,ierr)
    call mpi_bcast(iomap,size(iomap),mpi_integer,0,mpi_comm_world,ierr)
!---
    if ( myrank == 0 ) then
       open(3,file="wf.dat1",status="old",form="unformatted")
       read(3) ML_tmp,ML1_tmp,ML2_tmp,ML3_tmp
       read(3) MB_tmp,MB1_tmp,MB2_tmp
       close(3)
       write(*,*) "MB_0_IO,MB_1_IO=",1,MB_tmp
    end if
    call mpi_bcast(MB_tmp,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    MB_0_IO=1
    MB_1_IO=MB_tmp
!---
  END SUBROUTINE read_io2

  SUBROUTINE read_data_io2( disp_switch )
    implicit none
    logical,intent(IN) :: disp_switch
    integer,parameter :: u1=10
    integer :: i,n,k,s,i1,i2,i3,i0,a1a,b1a,a2a,b2a,a3a,b3a
    integer :: a1c,b1c,a2c,b2c,a3c,b3c,a1b,b1b,a2b,b2b,a3b,b3b
    integer,allocatable :: LLL(:,:,:)
    character(5) :: crank_old
    character(32) :: filename
    real(8) :: ct0,ct1,et0,et1
#ifdef _DRSDFT_
    real(8),allocatable :: w(:,:,:)
    real(8),parameter :: z0=0.0d0
#else
    complex(8),allocatable :: w(:,:,:)
    complex(8),parameter :: z0=(0.0d0,0.0d0)
#endif
    if ( disp_switch_parallel ) then
       write(*,'(a40," read_data_io2")') repeat("-",60)
    end if
    call watch(ct0,et0)
    call read_io2(myrank)
    a1b=Igrid(1,1)
    b1b=Igrid(2,1)
    a2b=Igrid(1,2)
    b2b=Igrid(2,2)
    a3b=Igrid(1,3)
    b3b=Igrid(2,3)
    allocate( LLL(a1b:b1b,a2b:b2b,a3b:b3b) ) ; LLL=0
    i=ML_0_WF-1
    do i3=a3b,b3b
    do i2=a2b,b2b
    do i1=a1b,b1b
       i=i+1
       LLL(i1,i2,i3)=i
    end do
    end do
    end do
    do i=1,nmap(myrank+1)
       write(crank_old,'(i5.5)') iomap(0,i,myrank+1)
       filename="wf.dat1."//crank_old
       if ( disp_switch_parallel ) write(*,*) filename
       open(u1,file=filename,form='unformatted')
       a1a=iomap(1,i,myrank+1)
       b1a=iomap(2,i,myrank+1)
       a2a=iomap(3,i,myrank+1)
       b2a=iomap(4,i,myrank+1)
       a3a=iomap(5,i,myrank+1)
       b3a=iomap(6,i,myrank+1)
       allocate( w(a1a:b1a,a2a:b2a,a3a:b3a) ) ; w=z0
       a1c=max(a1b,a1a)
       b1c=min(b1b,b1a)
       a2c=max(a2b,a2a)
       b2c=min(b2b,b2a)
       a3c=max(a3b,a3a)
       b3c=min(b3b,b3a)
       do s=MS_0_WF,MS_1_WF
       do k=MK_0_WF,MK_1_WF
       do n=MB_0_IO,MB_1_IO
          w=z0
          read(u1) w
          do i3=a3c,b3c
          do i2=a2c,b2c
          do i1=a1c,b1c
             i0=LLL(i1,i2,i3)
             unk(i0,n,k,s)=w(i1,i2,i3)
          end do
          end do
          end do
       end do ! n
       end do ! k
       end do ! s
       deallocate( w )
       close(u1)
    end do ! i
    deallocate( LLL )
    call watch(ct1,et1)
    if ( disp_switch_parallel ) then
       write(*,*) "time(read_data_io2)=",ct1-ct0,et1-et0
       write(*,'(a40," read_data_io2(end)")') repeat("-",60)
    end if
  END SUBROUTINE read_data_io2

END MODULE io2_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
	SUBROUTINE UPDSXY
C************************************************************************
C* UPDSXY								*
C* 									*
C* This subroutine updates the coordinate systems for the S level.	*
C* All higher levels are also updated.					*
C* 									*
C* UPDSXY								*
C**									*
C* Log:									*
C* C. Lin/EAI	 	6/97						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'XYDEF.CMN'
C-------------------------------------------------------------------------
C	Check to see if device is set before doing computations.
C
	IF ( ddev .ne. ' ' ) THEN
	   xbndls = FLOAT ( ibndls )
	   ybndbs = FLOAT ( ibndbs )
	   xbndrs = FLOAT ( ibndrs )
	   ybndts = FLOAT ( ibndts )
C
	END IF
C
C*	Call the next level update routine.
C
	CALL UPDDXY
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE IN_BDTA  ( iret )
C************************************************************************
C* IN_BDTA                                                              *
C*                                                                      *
C* This subroutine serves as a BLOCKDATA statement, initializing        *
C* variables in GEMPAK common blocks.  This subroutine is called        *
C* by IP_INIT.  If a GEMPAK program does not call IP_INIT, it must      *
C* call IN_BDTA directly.                                               *
C*                                                                      *
C* IN_BDTA  ( IRET )                                                    *
C*                                                                      *
C* Output parameters:                                                   *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C**                                                                     *
C* Log:                                                                 *
C* M. desJardins/GSFC    2/87                                           *
C* M. desJardins/GSFC    6/88   Documentation                           *
C* M. desJardins/GSFC    7/90   Initialize LUNGEM for lun flags         *
C* M. desJardins/NMC     1/92   Initiaize mbchan                        *
C* K. Brill/NMC         01/92   Remove CALL IN_BOLD                     *
C* M. desJardins/NMC     8/94   Initialize IESTAT                       *
C* C. Lin/EAI            8/95   Added IMFTYP                            *
C* S. Jacobs/NCEP        7/96   Added IFINVD                            *
C* K. Tyle/GSC          12/96   Replace IESTAT with call to ER_STAT;    *
C*                              added ER_INIT call                      *
C* R. Tian/SAIC          1/04   Added GD_INIT call                      *
C* T. Piper/SAIC        12/07   Removed GD_INIT                         *
C************************************************************************
      INCLUDE         'GEMPRM.PRM'
      INCLUDE         'GMBDTA.CMN'
C       INCLUDE         'ADBUFF.CMN'
      INCLUDE         'IMGDEF.CMN'
C------------------------------------------------------------------------
      iret= 0
C
C*      Initialize IM variables.
C
      imftyp = IFINVD
C
C*      Initialize DP variables.
C
      DO  i = 1, MMFLDP
         iflgdp (i) = 0
      END DO
C
C*      Initialize IP variables.
C
      lunp = 0
      lunf = 0
C
C*      Initialize PC variables.
C
      tablrd = .false.
      dsflg  = .false.
      stnflg = .false.
      pptabl = .false.
C
C*      Initialize SF, SN, GD DM file numbers and LUNs.
C
      DO  i = 1, MMFILE
         isfcfn (i) = -1
         isndfn (i) = -1
         igrdfn (i) = -1
         lundm (i)  = -1
      END DO
C
C*      Initialize DM variables
C
      kclast = 0
      DO i = 1, MCACHE
         kwrite (i) = .false.
         kcflno (i) = 0
         kcrecn (i) = 0
      END DO
C
C*      Initialize logical unit number flags.
C
      DO  i = 1, 10
         lungem (i) = 0
      END DO
C
C*      Initialize channel number for GPLT mailbox, so that an error
C*      can be generated for missing call to GINITP.
C
C        mbchan = -999
C
C*      Initialize ERCMN.CMN and ERCMN.H files.
C
C       CALL ER_STAT ( 0, 0, .false., iret )
C       CALL ER_INIT ( iret )
C*
      RETURN
      END

max_stars_repo_name: sgdecker/pygempak

max_issues_repo_name: sgdecker/pygempak

max_forks_repo_name: sgdecker/pygempak


text:
      subroutine rdaeb(x0,y0,a,b,alpha,xf,yf)
C    rdaeb     01/04/76      01/01/80
      common/gfgel/gf(28)
      if(gf(9).ge.0.) call strmod(0)
      i=abs(gf(9))
      k=gf(i)
      if(k+20) 5,4,8
  8   if(k+2) 5,4,5
  4   th01=gf(i+6)*0.01745329
      thf1=gf(i+7)*0.01745329
      alpha=gf(i+5)
      alpha1=gf(i+5)*0.01745329
      a=gf(i+3)
      b=gf(i+4)
      a1=abs(a)
      b1=abs(b)
      d0=a1*b1/sqrt((a1*sin(th01))**2+(b1*cos(th01))**2)
      df=a1*b1/sqrt((a1*sin(thf1))**2+(b1*cos(thf1))**2)
      x0=gf(i+1)+d0*cos(th01+alpha1)
      y0=gf(i+2)+d0*sin(th01+alpha1)
      xf=gf(i+1)+df*cos(thf1+alpha1)
      yf=gf(i+2)+df*sin(thf1+alpha1)
      go to 6
  5   call gffals(6)
  6   if(gf(9).gt.0.) call strmod(0)
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
2009 is being planned.

Recap2007 Festival

The second annual Burning Bike Festival was held Friday May 25th, 2007 at The Domes after the May Critical Mass ride. The festival was hosted by the ministers of the Bike Church.

2007 Events
 Tour de Domes II, drunken bicycle relay race
 Bicycle UnJousting competition
 Tall Bike Jousting bike jousting
 Bicycle Dancing competition
 Ceremonial bike burning (a Roadmaster)
 A Short Rant with a Twist (a skit by Chris Congleton)

Recap2006 Festival

The first ever Burning Bike Festival was held on the evening of May 26th, 2006 at The Domes after a Critical Mass ride. The festival was conceived and hosted by the ministers of the Bike Church.

2006 Events
 Tour de Domes drunken bicycle relay race
 Bicycle UnJousting competition
 Bicycle Dancing competition
 Ceremonial bike burning (a Magna)
 Greed versus the Ordinary Bicyclist (a skit by Chris Congleton)

See Bicycling Activities and Clubs for other bikerelated events.

20060915 12:57:38 nbsp awesome! Users/MarieHuynh

20080805 12:59:49 nbsp is this not happening this summer??? ( It was so fun last year. Users/EliseKane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*
*     ------------------------------------------------------------------
*       B L O C K   D A T A    G L C O N S
*     ------------------------------------------------------------------
*
      BLOCK DATA GLCONS
C
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
C
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
C
C     SET GLOBAL REAL CONSTANTS
C
      DATA ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS/
     :     0.0D 00,0.1D 00,0.5D 00,
     :     1.0D 00,2.0D 00,3.0D 00,
     :     4.0D 00,
     :     7.0D 00,1.1D 01,1.0D-08/
C
      END


text:
module brent
	integer, parameter :: tolerance = 1e-3

contains

	function root(x) result(r)
		real, dimension(10) :: x
		real :: r
		integer :: i

		r = 0

		do i = 1, 10
			r = r + x(i)
		end do

	end function

end module

max_stars_repo_name: nicolasessisbreton/fython

max_issues_repo_name: nicolasessisbreton/fython

max_forks_repo_name: nicolasessisbreton/fython


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

program nested_types

  use nested_types_module
  implicit none

  type(type2) :: obj
  integer :: int1 = 1
  real :: ret(1)

  obj%member%func1=>used_function
  ret = obj%member%func1(obj%member,int1)
  if (ret(1) .eq. 2.0) then
    print *, "PASS"
  else
    print *, "FAIL"
  endif

end program nested_types

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
      SUBROUTINE READNS(LUNIT,SUBSET,JDATE,IRET)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    READNS
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS SUBROUTINE SHOULD ONLY BE CALLED WHEN LOGICAL UNIT
C   LUNIT HAS BEEN OPENED FOR INPUT OPERATIONS.  IT READS THE NEXT
C   SUBSET FROM LOGICAL UNIT NUMBER LUNIT INTO INTERNAL SUBSET ARRAYS.
C   BUFR MESSAGES IN LUNIT MAY BE EITHER COMPRESSED OR UNCOMPRESSED.
C   THIS SUBROUTINE IS ACTUALLY A COMBINATION OF BUFR ARCHIVE LIBRARY
C   SUBROUTINES READMG AND READSB.
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 1998-07-08  J. WOOLLEN -- REPLACED CALL TO CRAY LIBRARY ROUTINE
C                           "ABORT" WITH CALL TO NEW INTERNAL BUFRLIB
C                           ROUTINE "BORT"
C 1999-11-18  J. WOOLLEN -- THE NUMBER OF BUFR FILES WHICH CAN BE
C                           OPENED AT ONE TIME INCREASED FROM 10 TO 32
C                           (NECESSARY IN ORDER TO PROCESS MULTIPLE
C                           BUFR FILES UNDER THE MPI)
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- MAXJL (MAXIMUM NUMBER OF JUMP/LINK ENTRIES)
C                           INCREASED FROM 15000 TO 16000 (WAS IN
C                           VERIFICATION VERSION); UNIFIED/PORTABLE FOR
C                           WRF; ADDED DOCUMENTATION (INCLUDING
C                           HISTORY); OUTPUTS MORE COMPLETE DIAGNOSTIC
C                           INFO WHEN ROUTINE TERMINATES ABNORMALLY
C
C USAGE:    CALL READNS (LUNIT, SUBSET, JDATE, IRET)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C
C   OUTPUT ARGUMENT LIST:
C     SUBSET   - CHARACTER*8: TABLE A MNEMONIC FOR BUFR MESSAGE
C                CONTAINING SUBSET BEING READ
C     JDATE    - INTEGER: DATE-TIME STORED WITHIN SECTION 1 OF BUFR
C                MESSAGE CONTAINING SUBSET BEING READ, IN FORMAT OF
C                EITHER YYMMDDHH OR YYYYMMDDHH, DEPENDING ON DATELEN()
C                VALUE
C     IREADNS  - INTEGER: RETURN CODE:
C                       0 = normal return
C                      -1 = there are no more subsets in the BUFR file
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT     READMG   READSB   STATUS
C    THIS ROUTINE IS CALLED BY: IREADNS
C                               Also called by application programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /MSGCWD/ NMSG(NFILES),NSUB(NFILES),MSUB(NFILES),
     .                INODE(NFILES),IDATE(NFILES)
      COMMON /BTABLES/ MAXTAB,NTAB,TAG(MAXJL),TYP(MAXJL),KNT(MAXJL),
     .                JUMP(MAXJL),LINK(MAXJL),JMPB(MAXJL),
     .                IBT(MAXJL),IRF(MAXJL),ISC(MAXJL),
     .                ITP(MAXJL),VALI(MAXJL),KNTI(MAXJL),
     .                ISEQ(MAXJL,2),JSEQ(MAXJL)

      CHARACTER*10 TAG
      CHARACTER*8  SUBSET
      CHARACTER*3  TYP

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

C  REFRESH THE SUBSET AND JDATE PARAMETERS
C  ---------------------------------------

      CALL STATUS(LUNIT,LUN,IL,IM)
      IF(IL.EQ.0) GOTO 900
      IF(IL.GT.0) GOTO 901
      SUBSET = TAG(INODE(LUN))
      JDATE  = IDATE(LUN)

C  READ THE NEXT SUBSET IN THE BUFR FILE
C  -------------------------------------

1     CALL READSB(LUNIT,IRET)
      IF(IRET.NE.0) THEN
         CALL READMG(LUNIT,SUBSET,JDATE,IRET)
         IF(IRET.EQ.0) GOTO 1
      ENDIF

C  EXITS
C  -----

      RETURN
900   CALL BORT('BUFRLIB: READNS - INPUT BUFR FILE IS CLOSED, IT MUST'//
     . ' BE OPEN FOR INPUT')
901   CALL BORT('BUFRLIB: READNS - INPUT BUFR FILE IS OPEN FOR OUTPUT'//
     . ', IT MUST BE OPEN FOR INPUT')
      END

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
      subroutine adprpl(actw,aw0plo,aw0prn,aw1plo,aw1prn,dlaplo,
     $ dlaprn,dleplo,dleprn,dlhplo,dlhprn,dloplo,dloprn,dltpll,
     $ dltplo,dltprl,dltprn,dlxpll,dlxplo,dlxprl,dlxprn,dlxdmp,
     $ eh,eh0plo,eh0prn,eh1plo,eh1prn,fo2lg,fxprpl,o20plo,o20prn,
     $ o21plo,o21prn,ph,ph0plo,ph0prn,ph1plo,ph1prn,qplaw0,qplaw1,
     $ qpleh0,qpleh1,qplolt,qplolx,qplott,qplotx,qplo20,qplo21,
     $ qplph0,qplph1,qpraw0,qpraw1,qpreh0,qpreh1,qprnlx,qprnlt,
     $ qprntx,qprntt,qpro20,qpro21,qprph0,qprph1,qredox,time1,
     $ tiplol,tiplot,tiprnl,tiprnt,xidump,xiplol,xiplot,xiprnl,
     $ xiprnt,xi1)
c
c     This subroutine adjusts print, plot, and PRS transfer points
c     as needed for the next step along the reaction path.
c
c     This subroutine is called by:
c
c       EQ6/path.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      logical qplaw0,qplaw1,qpleh0,qpleh1,qplolt,qplolx,qplott,qplotx,
     $ qplo20,qplo21,qplph0,qplph1,qpraw0,qpraw1,qpreh0,qpreh1,qprnlx,
     $ qprnlt,qprntx,qprntt,qpro20,qpro21,qprph0,qprph1,qredox
c
      real*8 actw,aw0plo,aw0prn,aw1plo,aw1prn,dlaplo,dlaprn,dleplo,
     $ dleprn,dlhplo,dlhprn,dloplo,dloprn,dltpll,dltplo,dltprl,dltprn,
     $ dlxpll,dlxplo,dlxprl,dlxprn,dlxdmp,eh,eh0plo,eh0prn,eh1plo,
     $ eh1prn,fo2lg,fxprpl,o20plo,o20prn,o21plo,o21prn,ph,ph0plo,
     $ ph0prn,ph1plo,ph1prn,time1,tiplol,tiplot,tiprnl,tiprnt,xidump,
     $ xiplol,xiplot,xiprnl,xiprnt,xi1
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      real*8 tolx,tx,xx
c
      real*8 texp,tlg
c
c-----------------------------------------------------------------------
c
c     Reset print points as necessary.
c
      if (qprntx) then
        tolx = fxprpl*dlxprn
  100   xiprnt = xiprnt + dlxprn
        if ((xi1 + tolx) .ge. xiprnt) go to 100
      endif
c
      if (qprntt) then
        tolx = fxprpl*dlxprn
  110   tiprnt = tiprnt + dltprn
        if ((time1 + tolx) .ge. tiprnt) go to 110
      endif
c
      if (qprnlx) then
        tolx = fxprpl*dlxprl
  120   xiprnl = texp(tlg(xiprnl) + dlxprl)
        xx = texp(tlg(xi1) + tolx)
        if (xx .ge. xiprnl) go to 120
      endif
c
      if (qprnlt) then
        tolx = fxprpl*dltprl
  130   tiprnl = texp(tlg(tiprnl) + dltprl)
        tx = texp(tlg(time1) + tolx)
        if (tx .ge. tiprnl) go to 130
      endif
c
      if (qprph0) then
        tolx = fxprpl*dlhprn
  140   ph0prn = ph0prn - dlhprn
        if ((ph - tolx) .le. ph0prn) go to 140
        ph1prn = ph0prn
  150   ph1prn = ph1prn + dlhprn
        if ((ph + tolx) .ge. ph1prn) go to 150
      endif
c
      if (qprph1) then
        tolx = fxprpl*dlhprn
  160   ph1prn = ph1prn + dlhprn
        if ((ph + tolx) .ge. ph1prn) go to 160
        ph0prn = ph1prn
  170   ph0prn = ph0prn - dlhprn
        if ((ph - tolx) .le. ph0prn) go to 170
      endif
c
      if (qredox) then
        if (qpreh0) then
          tolx = fxprpl*dleprn
  180     eh0prn = eh0prn - dleprn
          if ((eh - tolx) .le. eh0prn) go to 180
          eh1prn = eh0prn
  190     eh1prn = eh1prn + dleprn
          if ((eh + tolx) .ge. eh1prn) go to 190
        endif
c
        if (qpreh1) then
          tolx = fxprpl*dleprn
  200     eh1prn = eh1prn + dleprn
          if ((eh + tolx) .ge. eh1prn) go to 200
          eh0prn = eh1prn
  210     eh0prn = eh0prn - dleprn
          if ((eh - tolx) .le. eh0prn) go to 210
        endif
c
        if (qpro20) then
          tolx = fxprpl*dloprn
  220     o20prn = o20prn - dloprn
          if ((fo2lg - tolx) .le. o20prn) go to 220
          o21prn = o20prn
  230     o21prn = o21prn + dloprn
          if ((fo2lg + tolx) .ge. o21prn) go to 230
        endif
c
        if (qpro21) then
          tolx = fxprpl*dloprn
  240     o21prn = o21prn + dloprn
          if ((fo2lg + tolx) .ge. o21prn) go to 240
          o20prn = o21prn
  250     o20prn = o20prn - dloprn
          if ((fo2lg - tolx) .le. o20prn) go to 250
        endif
c
      endif
c
      if (qpraw0) then
        tolx = fxprpl*dlaprn
  260   aw0prn = aw0prn - dlaprn
        if ((actw - tolx) .le. aw0prn) go to 260
        aw1prn = aw0prn
  270   aw1prn = aw1prn + dlaprn
        if ((actw + tolx) .ge. aw1prn) go to 270
      endif
c
      if (qpraw1) then
        tolx = fxprpl*dlaprn
  280   aw1prn = aw1prn + dlaprn
        if ((actw + tolx) .ge. aw1prn) go to 280
        aw0prn = aw1prn
  290   aw0prn = aw0prn - dlaprn
        if ((actw - tolx) .le. aw0prn) go to 290
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Reset plot points as necessary.
c
      if (qplotx) then
        tolx = fxprpl*dlxplo
  400   xiplot = xiplot + dlxplo
        if ((xi1 + tolx)  .ge. xiplot) go to 400
      endif
c
      if (qplott) then
        tolx = fxprpl*dltplo
  410   tiplot = tiplot + dltplo
        if ((time1 + tolx)  .ge. tiplot) go to 410
      endif
c
      if (qplolx) then
        tolx = fxprpl*dlxpll
  420   xiplol = texp(tlg(xiplol) + dlxpll)
        xx = texp(tlg(xi1) + tolx)
        if (xx  .ge. xiplol) go to 420
      endif
c
      if (qplolt) then
        tolx = fxprpl*dltpll
  430   tiplol = texp(tlg(tiplol) + dltpll)
        tx = texp(tlg(time1) + tolx)
        if (tx .ge. tiplol) go to 430
      endif
c
      if (qplph0) then
        tolx = fxprpl*dlhplo
  440   ph0plo = ph0plo - dlhplo
        if ((ph - tolx) .le. ph0plo) go to 440
        ph1plo = ph0plo
  450   ph1plo = ph1plo + dlhplo
        if ((ph + tolx) .ge. ph1plo) go to 450
      endif
c
      if (qplph1) then
        tolx = fxprpl*dlhplo
  460   ph1plo = ph1plo + dlhplo
        if ((ph + tolx) .ge. ph1plo) go to 460
        ph0plo = ph1plo
  470   ph0plo = ph0plo - dlhplo
        if ((ph - tolx) .le. ph0plo) go to 470
      endif
c
      if (qredox) then
        if (qpleh0) then
          tolx = fxprpl*dleplo
  480     eh0plo = eh0plo - dleplo
          if ((eh - tolx) .le. eh0plo) go to 480
          eh1plo = eh0plo
  490     eh1plo = eh1plo + dleplo
          if ((eh + tolx) .ge. eh1plo) go to 490
        endif
c
        if (qpleh1) then
          tolx = fxprpl*dleplo
  500     eh1plo = eh1plo + dleplo
          if ((eh + tolx) .ge. eh1plo) go to 500
          eh0plo = eh1plo
  510     eh0plo = eh0plo - dleplo
          if ((eh - tolx) .le. eh0plo) go to 510
        endif
c
        if (qplo20) then
          tolx = fxprpl*dloplo
  520     o20plo = o20plo - dloplo
          if ((fo2lg - tolx) .le. o20plo) go to 520
          o21plo = o20plo
  530     o21plo = o21plo + dloplo
          if ((fo2lg + tolx) .ge. o21plo) go to 530
        endif
c
        if (qplo21) then
          tolx = fxprpl*dloplo
  540     o21plo = o21plo + dloplo
          if ((fo2lg + tolx) .ge. o21plo) go to 540
          o20plo = o21plo
  550     o20plo = o20plo - dloplo
          if ((fo2lg - tolx) .le. o20plo) go to 550
        endif
c
      endif
c
      if (qplaw0) then
        tolx = fxprpl*dlaplo
  560   aw0plo = aw0plo - dlaplo
        if ((actw - tolx) .le. aw0plo) go to 560
        aw1plo = aw0plo
  570   aw1plo = aw1plo + dlaplo
        if ((actw + tolx) .ge. aw1plo) go to 570
      endif
c
      if (qplaw1) then
        tolx = fxprpl*dlaplo
  580   aw1plo = aw1plo + dlaplo
        if ((actw + tolx) .ge. aw1plo) go to 580
        aw0plo = aw1plo
  590   aw0plo = aw0plo - dlaplo
        if ((actw - tolx) .le. aw0plo) go to 590
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Reset the PRS transfer point as necessary.
c
      if (xidump. le. xi1) then
        tolx = fxprpl*dlxdmp
  700   xidump = xidump + dlxdmp
        if ((xi1 + tolx)  .ge. xidump) go to 700
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


text:
      DOUBLE PRECISION FUNCTION ZLA_GBRPVGRW( N, KL, KU, NCOLS, AB,
     $                                        LDAB, AFB, LDAFB )
*
*     -- LAPACK routine (version 3.2.1)                                 --
*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
*     -- Jason Riedy of Univ. of California Berkeley.                 --
*     -- April 2009                                                   --
*
*     -- LAPACK is a software package provided by Univ. of Tennessee, --
*     -- Univ. of California Berkeley and NAG Ltd.                    --
*
      IMPLICIT NONE
*     ..
*     .. Scalar Arguments ..
      INTEGER            N, KL, KU, NCOLS, LDAB, LDAFB
*     ..
*     .. Array Arguments ..
      COMPLEX*16         AB( LDAB, * ), AFB( LDAFB, * )
*     ..
*
*  Purpose
*  =======
*
*  ZLA_GBRPVGRW computes the reciprocal pivot growth factor
*  norm(A)/norm(U). The "max absolute element" norm is used. If this is
*  much less than 1, the stability of the LU factorization of the
*  (equilibrated) matrix A could be poor. This also means that the
*  solution X, estimated condition numbers, and error bounds could be
*  unreliable.
*
*  Arguments
*  =========
*
*     N       (input) INTEGER
*     The number of linear equations, i.e., the order of the
*     matrix A.  N >= 0.
*
*     KL      (input) INTEGER
*     The number of subdiagonals within the band of A.  KL >= 0.
*
*     KU      (input) INTEGER
*     The number of superdiagonals within the band of A.  KU >= 0.
*
*     NCOLS   (input) INTEGER
*     The number of columns of the matrix A.  NCOLS >= 0.
*
*     AB      (input) COMPLEX*16 array, dimension (LDAB,N)
*     On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
*     The j-th column of A is stored in the j-th column of the
*     array AB as follows:
*     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
*
*     LDAB    (input) INTEGER
*     The leading dimension of the array AB.  LDAB >= KL+KU+1.
*
*     AFB     (input) COMPLEX*16 array, dimension (LDAFB,N)
*     Details of the LU factorization of the band matrix A, as
*     computed by ZGBTRF.  U is stored as an upper triangular
*     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
*     and the multipliers used during the factorization are stored
*     in rows KL+KU+2 to 2*KL+KU+1.
*
*     LDAFB   (input) INTEGER
*     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, J, KD
      DOUBLE PRECISION   AMAX, UMAX, RPVGRW
      COMPLEX*16         ZDUM
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN, REAL, DIMAG
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Statement Function Definitions ..
      CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
*     ..
*     .. Executable Statements ..
*
      RPVGRW = 1.0D+0

      KD = KU + 1
      DO J = 1, NCOLS
         AMAX = 0.0D+0
         UMAX = 0.0D+0
         DO I = MAX( J-KU, 1 ), MIN( J+KL, N )
            AMAX = MAX( CABS1( AB( KD+I-J, J ) ), AMAX )
         END DO
         DO I = MAX( J-KU, 1 ), J
            UMAX = MAX( CABS1( AFB( KD+I-J, J ) ), UMAX )
         END DO
         IF ( UMAX /= 0.0D+0 ) THEN
            RPVGRW = MIN( AMAX / UMAX, RPVGRW )
         END IF
      END DO
      ZLA_GBRPVGRW = RPVGRW
      END FUNCTION

max_issues_repo_name: cpatrick/calatk

max_forks_repo_name: cpatrick/calatk


text:
      SUBROUTINE ICHKZ (CERIND, FUNCID)

C  ICHKZ does post-processing after an error-indicator function has
C  been called.

C  Input parameters:
C    CERIND: actual error indicator returned by function
C    FUNCID: numeric identifier of function that returned error indicator

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

      INTEGER     IX, FUNCID, ITRIM, CERIND
      CHARACTER   FUNNAM*80, TMSG*500, MSG*160, ACODE*4
      INTEGER     IARFND

      IF (EXPSIZ.LT.1) THEN
         CALL UNMSG ('ICHKZ is called with invalid EXPSIZ.')
      ENDIF

C  set up test message
      CALL NFUNM (FUNCID, FUNNAM)
      TMSG = 'If ' // FUNNAM
      TMSG(ITRIM(TMSG)+2:) = 'is called and ' // CURCON
      TMSG(ITRIM(TMSG)+1:) = ', it should return error indicator = '

C  may be several valid error code possibilities
      DO 50 IX = 1, EXPSIZ-1
         WRITE (ACODE, '(I4.3)') EXPERR(IX)
         TMSG(ITRIM(TMSG)+2:) = ACODE // ','
50    CONTINUE
      IF (EXPSIZ.GE.2) TMSG(ITRIM(TMSG):) = ' or '
      WRITE(ACODE, '(I4.3)') EXPERR(EXPSIZ)
      TMSG(ITRIM(TMSG)+2:) = ACODE // '.'
      CALL SETMSG (ERRSRS, TMSG)

      IF (IARFND(CERIND, EXPSIZ, EXPERR) .GT. 0) THEN
         CALL PASS
      ELSE
         CALL FAIL
         WRITE (MSG, '(A,I6.3,A)') 'Reported error indicator ',
     1          CERIND, ' is not one of the expected values.'
         CALL INMSG (MSG)
      ENDIF

      END


text:
!AC -- read KORAL files 2/2/17
!   -- modified from fluid_model_harm.f90
   module fluid_model_koral3d

      use class_four_vector
      use interpolate, only: interp
      use kerr, only: kerr_metric, lnrf_frame, uks2ubl, bl2ks
      use phys_constants, only: pi
      use math, only: zbrent
      implicit none

      !AC koral will read modified harm .in files for now
      !AC should change to special koral namelist
      namelist /harm/  dfile, hfile, nt, indf

      logical :: doingkoralnth,shortfile
      character(len=100) :: dfile, hfile
      integer :: nx1, nx2, nx3,n, ndumps, nt, indf, dlen, nhead
      integer :: nrelbin=0
      integer :: ndim=3
      integer :: minpolecell=2
      integer, dimension(:), allocatable :: dumps
      real :: tstep=20., h, r0, aa, bb, pp, asim, dx1, dx2, gam, startx1, startx2, toffset=0.
      real :: relgammamin=1., relgammamax=1.
      real(8) :: scalefac
      real, dimension(:), allocatable :: t
      real, dimension(:), allocatable :: x1_arr, x2_arr, x3_arr,r_arr, th_arr,ph_arr
      real, dimension(:), allocatable :: rho_arr, p_arr, Be_arr, u0_arr, vrl_arr, &
           vpl_arr, vtl_arr
      real, dimension(:,:), allocatable :: nth_arr
      real, dimension(:), allocatable :: b0_arr, br_arr, bth_arr, bph_arr

!      shortfile=.TRUE.

      interface init_koral3d_data
        module procedure init_koral3d_data
      end interface

      interface del_koral3d_data
        module procedure del_koral3d_data
      end interface

      interface initialize_koral3d_model
        module procedure initialize_koral3d_model
      end interface

      interface transformbl2mksh3
         module procedure transformbl2mksh3
      end interface

      interface transformmksh32bl
         module procedure transformmksh32bl
      end interface

      interface koral3d_vals
        module procedure koral3d_vals
      end interface

      contains

        ! transform Boyer-Lindquist theta coordinates to modified Kerr-Schild x2 at a corresponding array of r
        subroutine transformbl2mksh3(th,r,x2,h,a,b,p)

        real, intent(in), dimension(:) :: th,r
        real, intent(out), dimension(size(th)) :: x2
        real, intent(in) :: h,a,b,p

        x2=0.5*(1 + ((r**p)/(h*pi)) * (atan(tan(0.5*h*pi)*(1-2*th/pi))/((b-a)*(2**p)+(a-0.5)*(r**p))))
        end subroutine transformbl2mksh3

        ! transform modified Kerr-Schild x2 to Boyer-Lindquist theta at a corresponding array of r=r0+exp(x1)
        subroutine transformmksh32bl(x2,r,th,h,a,b,p)

        real, intent(in), dimension(:) :: x2,r
        real, intent(out), dimension(size(x2)) :: th
        real, intent(in) :: h,a,b,p

        th = 0.5*pi*(1 + tan(h*pi*(-0.5+x2+(1-2*x2)*(a+(2**p)*(b-a)/(r**p))))/tan(0.5*h*pi))
        !write(6,*) 'x2len: ', size(x2)
        end subroutine transformmksh32bl

        ! Interpolates KORAL data to input coordinates
        subroutine koral3d_vals(x0,a,rho,p,b,u,bmag,Be,nth,type)

        type (four_Vector), intent(in), dimension(:) :: x0
        real, intent(in) :: a
        integer, intent(in) :: type
        real(kind=8), dimension(size(x0)) :: done,nfac,pfac
        real, dimension(size(x0)) :: x3, x2,x1,zm,theta,fone, &
         vpl0,vrl0,vtl0,rd,td,pd,rttd,zr,dzero, &
         vr0,vth0,vph0,bth,dummy,tt,ttd,zt,zphi,zpp
        real, dimension(nx1) :: uniqx1,uniqr
        real, dimension(nx2) :: uniqx2
        real, dimension(nx3) :: uniqx3,uniqph
        real, dimension(size(x0),2**(ndim+1)) :: ppi,Bei,rhoi,vrli,vtli, &
             vpli,bri,bthi,bphi,b0i,u0i, nth_bini

        real, dimension(size(x0)) :: dth,thu,thl,thone,minph
        integer, dimension(size(x0)*2*(2**ndim)) :: indx
        integer, dimension(size(x0)) :: lx1,lx2,lx3,ux1,ux2,ux3,x1l,x1u,x2l,x2u, x3l,x3u, &
         one,tindx,lx2l,lx2u,ux2l,ux2u,lx3l,lx3u,ux3l,ux3u,low,high
        real :: x2mintrust, x2maxtrust,thmintrust,thmaxtrust
        integer :: npts,i,j,maxtindx,ie
        real, dimension(size(x0)) :: b1tmp,b2tmp,b3tmp,b4tmp,u1tmp,ntmp
        real, dimension(size(x0)), intent(out) :: rho,p,bmag,Be
        real, dimension(size(x0),nrelbin), intent(out) :: nth
        type (four_Vector), intent(out), dimension(size(x0)) :: u,b

        one=1
        dzero=0d0; done=1d0; fone=1

        !Get unique code coordinates x1,x2,x3
        uniqx2=x2_arr(1:nx2)
        uniqx1=x1_arr(1:nx2*(nx1-1)+1:nx2)
        uniqx3=x3_arr(1:nx2*(nx1-1)*nx3+1:nx2*nx1)

        !Get unique r,phi Boyer-Lindquist coordinates
        uniqr=r0+exp(uniqx1)
        uniqph=uniqx3
        !there is no uniqth because in mks3, theta grid is r dependent!

        !Transform to Kerr-Schild x1,x2,x3 coordinates
        npts=size(x0)
        theta=x0%data(3)
        zr=x0%data(2)
        zt=x0%data(1)
        tt=bl2ks(dble(zr),dble(zt),dble(a),1) - bl2ks(dble(zr(1)),0d0,dble(a),1)
        tt=-tt+1d-8

        zpp=x0%data(4)
        !!zphi=zpp
        !!zphi=mod(zphi,(2.*pi))

        !!  AC this is necessary for koral dumps that keep phi output in ks!!
        zphi=bl2ks(dble(zr), dble(zpp), dble(a))
        zphi=mod(zphi,(2.*pi))
        if(any(isnan(zphi))) then
           write(6,*) 'nan in zphi!'
        end if

        where(zphi.lt.0.)
            zphi=zphi+2.*pi
        endwhere
        where(zphi.gt.pi)
            zphi=zphi-2.*pi
        endwhere

        x1 = log(dble(zr)-r0)
        x3 = zphi
        call transformbl2mksh3(theta,zr,x2,h,aa,bb,pp)

        !Distance from current points to grid points
        lx1=floor((x1-uniqx1(1))/(uniqx1(nx1)-uniqx1(1))*(nx1-1))+1
        ux1=lx1+1
        lx1=merge(lx1,one,lx1.ge.1)
        ux1=merge(ux1,nx1,ux1.le.nx1)

        lx2=floor((x2-uniqx2(1))/(uniqx2(nx2)-uniqx2(1))*(nx2-1))+1
        ux2=lx2+1
        lx2=merge(lx2,one,lx2.ge.1)
        ux2=merge(ux2,nx2,ux2.le.nx2)

        lx3=floor((x3-uniqx3(1))/(uniqx3(nx3)-uniqx3(1))*(nx3-1))+1
        ux3=lx3+1
        !lx3=merge(lx3,one,lx3.ge.1)
        !ux3=merge(ux3,nx3,ux3.le.nx3)

        ! phi-distance
        ! periodic in phi
        minph=uniqph(lx3)
        ux3=merge(ux3,one,ux3.le.nx3)
        !where(ux3.gt.nx3)
        !    ux3=1
        !endwhere
        where(lx3.lt.1)
            lx3=nx3
            minph=uniqph(lx3)-2.*pi
        endwhere
        pd=(zphi-minph)/(uniqph(2)-uniqph(1))

        ! r-distance
        rd=(zr-uniqr(lx1))/(uniqr(ux1)-uniqr(lx1))

        ! When the geodesic is inside the innermost zone center located outside the horizon,
        ! use nearest neighbor rather than interpolate
        where((uniqr(lx1).le.(1.+sqrt(1.-a**2.))).or.(lx1.eq.1))
           rd=1.
           pfac=1e-3
           nfac=1e-3
        elsewhere
           pfac=1.
           nfac=1.
        endwhere

        ! theta-distance
        ! Deal with poles
        !AC for interpolation, choose rl for theta distance?
        call transformmksh32bl(uniqx2(lx2),uniqr(lx1),thl,h,aa,bb,pp)
        call transformmksh32bl(uniqx2(ux2),uniqr(lx1),thu,h,aa,bb,pp)
        call transformmksh32bl(uniqx2(one),uniqr(lx1),thone,h,aa,bb,pp)

        where(ux2.ne.lx2)
           td=abs((theta-thl)/(thu-thl))
        elsewhere
           td=abs((theta-thl)/thone)
        endwhere

        ! Indices for nearest neighbors in 3D
        ! th is fastest changing index
        x2l=lx2-1;
        x2u=ux2-1
        x1l=(lx1-1)*nx2 ;
        x1u=(ux1-1)*nx2
        x3l=(lx3-1)*nx2*nx1;
        x3u=(ux3-1)*nx2*nx1

        ! Now find timestep indices:
        maxtindx=nt-2
        where(floor(tt/tstep).le.maxtindx)
           tindx=floor(tt/tstep)
           ttd=(tt-tindx*tstep)/tstep
        elsewhere
           tindx=maxtindx+1
           ttd=0.
        endwhere

        ! Create index array across times and 3D nearest neighbors (first half at first time, second half at second time):
        ! Need to make sure there aren't integer problems when these indices get very large
        indx=(/(/x1l+x2l+x3l+1+n*tindx/),(/x1l+x2u+x3l+1+n*tindx/),(/x1u+x2l+x3l+1+n*tindx/), &
            (/x1u+x2u+x3l+1+n*tindx/),(/x1l+x2l+x3u+1+n*tindx/), &
            (/x1l+x2u+x3u+1+n*tindx/),(/x1u+x2l+x3u+1+n*tindx/),(/x1u+x2u+x3u+1+n*tindx/), &
            (/x1l+x2l+x3l+1+n*(tindx+1)/),(/x1l+x2u+x3l+1+n*(tindx+1)/),(/x1u+x2l+x3l+1+n*(tindx+1)/), &
            (/x1u+x2u+x3l+1+n*(tindx+1)/),(/x1l+x2l+x3u+1+n*(tindx+1)/),(/x1l+x2u+x3u+1+n*(tindx+1)/), &
            (/x1u+x2l+x3u+1+n*(tindx+1)/),(/x1u+x2u+x3u+1+n*(tindx+1)/)/)

        ! keep indx from going out of bounds when loading one time slice
        where(indx.gt.size(rho_arr))
           indx=indx-n
        endwhere

        rhoi=reshape(rho_arr(indx),(/npts,2**(ndim+1)/))
        vrli=reshape(vrl_arr(indx),(/npts,2**(ndim+1)/))
        vtli=reshape(vtl_arr(indx),(/npts,2**(ndim+1)/))
        vpli=reshape(vpl_arr(indx),(/npts,2**(ndim+1)/))
        ppi=reshape(p_arr(indx),(/npts,2**(ndim+1)/))
        Bei=reshape(Be_arr(indx),(/npts,2**(ndim+1)/))
        b0i=reshape(b0_arr(indx),(/npts,2**(ndim+1)/))
        bri=reshape(br_arr(indx),(/npts,2**(ndim+1)/))
        bthi=reshape(bth_arr(indx),(/npts,2**(ndim+1)/))
        bphi=reshape(bph_arr(indx),(/npts,2**(ndim+1)/))
        u0i=reshape(u0_arr(indx),(/npts,2**(ndim+1)/))

        !interpolate data
        rttd=0.
        rho=interp(rhoi,rttd,pd,rd,td)
        p=interp(ppi,rttd,pd,rd,td)
        Be=interp(Bei,rttd,pd,rd,td)
        vrl0=interp(vrli,rttd,pd,rd,td)
        vtl0=interp(vtli,rttd,pd,rd,td)
        b1tmp=interp(b0i,rttd,pd,rd,td)
        b2tmp=interp(bri,rttd,pd,rd,td)
        b3tmp=interp(bthi,rttd,pd,rd,td)
        b4tmp=interp(bphi,rttd,pd,rd,td)
        u1tmp=interp(u0i,rttd,pd,rd,td)
        vpl0=interp(vpli,rttd,pd,rd,td)

        !mask data outside trusted region
        rho=merge(rho,dzero,x1.gt.uniqx1(1))*nfac
        p=merge(p,fone,x1.gt.uniqx1(1))*pfac
        Be=merge(Be,dzero,x1.gt.uniqx1(1)) !AC -- does using zero make sense?
        vrl0=merge(vrl0,dzero,x1.gt.uniqx1(1))
        vtl0=merge(vtl0,fone,x1.gt.uniqx1(1))
        b1tmp=merge(b1tmp,fone,x1.gt.uniqx1(1))
        b2tmp=merge(b2tmp,fone,x1.gt.uniqx1(1))
        b3tmp=merge(b3tmp,fone,x1.gt.uniqx1(1))
        b4tmp=merge(b4tmp,fone,x1.gt.uniqx1(1))
        u1tmp=merge(dble(u1tmp),done,x1.gt.uniqx1(1))
        vpl0=merge(vpl0,dzero,x1.gt.uniqx1(1))

        !AC TODO !!!! minimum polar cell - change!
        !AC this is the minimum polar coordinate x2 we trust.
        x2mintrust = uniqx2(minpolecell)
        x2maxtrust = uniqx2(size(uniqx2)-(minpolecell-1))
        thmintrust = 0.
        thmaxtrust = pi

        !cut out 15 degrees away from poles for disk (also cut out high Be in fluid.f90)
        if (type.eq.1) then
           thmintrust = 0.25
           thmaxtrust = pi - 0.25

        !cut out background entirely -- all theta > pi/2 !AC -- assumes symmetric grid!!
        else if(type.eq.2) then
           thmaxtrust = 0.5*pi!0.25*pi.   !uniqx2(int(size(uniqx2)/2)-2)

        !cut out foreground entirely -- all theta > pi/2 !AC -- assumes symmetric grid!!
        else if(type.eq.3) then
           thmintrust = 0.5*pi!0.75*pi. !uniqx2(int(size(uniqx2)/2)+2)
        end if

        !theta based cuts for jet
        if((type.eq.1).or.(type.eq.2).or.(type.eq.3)) then
            rho=merge(rho,dzero,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            p=merge(p,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            Be=merge(Be,dzero,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))!AC -- does using zero make sense?
            vrl0=merge(vrl0,dzero,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            vtl0=merge(vtl0,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            b1tmp=merge(b1tmp,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            b2tmp=merge(b2tmp,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            b3tmp=merge(b3tmp,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            b4tmp=merge(b4tmp,fone,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            u1tmp=merge(dble(u1tmp),done,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
            vpl0=merge(vpl0,dzero,(theta.gt.thmintrust).and.(theta.lt.thmaxtrust))
        end if

        !x2 based cuts for polar axis
        rho=merge(rho,dzero,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        p=merge(p,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        Be=merge(Be,dzero,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust)) !AC -- does using zero make sense?
        vrl0=merge(vrl0,dzero,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        vtl0=merge(vtl0,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        b%data(1)=merge(b1tmp,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        b%data(2)=merge(b2tmp,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        b%data(3)=merge(b3tmp,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        b%data(4)=merge(b4tmp,fone,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        u%data(1)=merge(dble(u1tmp),done,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))
        vpl0=merge(vpl0,dzero,(x2.gt.x2mintrust).and.(x2.lt.x2maxtrust))

        !AC load nonthermal
        !AC do we have to do this in a loop?
        if(doingkoralnth) then
           do ie=1, nrelbin
              nth_bini=reshape(nth_arr(indx,ie),(/npts,2**(ndim+1)/))
              ntmp=interp(nth_bini,rttd,rd,td)
              ntmp=merge(ntmp,dzero,x1.gt.uniqx1(1))
              nth(:,ie)=merge(ntmp,dzero,x2.gt.x2mintrust)
           end do
        end if

        ! Compute magnitude of interpolated b-field and force b^2 > 0 (need to look into numerical issues here):
        call assign_metric(b,transpose(kerr_metric(zr, real(x0%data(3)),a)))
        bmag=b*b; bmag=merge(bmag,dzero,bmag.ge.0d0)
        bmag=sqrt(bmag)

        ! Protect azimuthal velocities at poles
        call lnrf_frame(vrl0,vtl0,vpl0,zr,a,real(x0%data(3)),vr0,vth0,vph0,1)
        u%data(2)=u%data(1)*dble(vr0)
        u%data(3)=u%data(1)*dble(vth0)
        u%data(4)=u%data(1)*dble(vph0)
        call assign_metric(u,transpose(kerr_metric(zr,real(x0%data(3)) ,a)))

        end subroutine koral3d_vals

        ! Read inputs from file
        !AC koral will share the harm input file structure for now
        subroutine read_koral3d_inputs(ifile)
        character(len=20), intent(in) :: ifile
        open(unit=8,file=ifile,form='formatted',status='old')
        read(8,nml=harm)
        write(6,*) 'read: ',nt
        close(unit=8)
        end subroutine read_koral3d_inputs

        ! Read header values
        subroutine read_koral3d_data_header(nhead)
        integer, intent(in) :: nhead
        real(8), dimension(nhead) :: header
        real(8), dimension(3) :: header2
        real(8) :: tcur

        write(6,*) 'read header: ', nhead
        open(unit=8,file=hfile,form='formatted',status='old')
        read(8,*) header

        tcur=header(1)
        nx1=header(2)
        nx2=header(3)
        nx3=header(4)
        asim=header(5)
        !mbh=header(6)
        r0=header(7)
        h=header(8)
        aa=header(9)
        bb=header(10)
        pp=header(11)

        !AC nonthermal
        if(doingkoralnth) then
           read(8,*) header2
           nrelbin=header2(1)
           relgammamin=header2(2)
           relgammamax=header2(3)
           write(6,*) 'read nth header vals ', nrelbin, relgammamin, relgammamax

        end if

        write(6,*) 'read header vals: ', asim, r0, h, aa, bb, pp
        close(unit=8)
        end subroutine read_koral3d_data_header

        !Read text files with  data values
        !AC TODO can we read binary data files instead???
        subroutine read_koral3d_data_file(data_file,tcur,rho,p,u,b,Be,nnth)
        character(len=100), intent(in) :: data_file
        integer :: nhead, nhead2, dlen
        integer :: rhopos,ppos,Bepos,vpos,bpos,gdetpos,gridpos,nthpos
        real(8), intent(out) :: tcur
        real(8), dimension(:), allocatable, intent(out) :: p,rho,Be
        real(8), dimension(:,:), allocatable, intent(out) :: nnth
        real(8), dimension(:), allocatable :: gdet, header, header2
        type (four_vector), dimension(:), allocatable, intent(out) :: u,b
        real(8), dimension(:,:), allocatable :: grid, data, tmetric
        integer :: i, nelem

        shortfile = .FALSE.

          !AC TODO still need to change these positions for 3d sim files!!
          if(nrelbin>0) then
             dlen=41+nrelbin;
          else
             dlen = 38
          endif

          if(shortfile) dlen = 23

          nhead=11; nhead2=3
          allocate(header(nhead))
          allocate(header2(nhead2))

          ! M87 library vs. Andrew format difference
          ! FOR SHORTFILE BE IS THE ION TEMPERATURE TO RECONSTRUCT P FOR POST-PROCESSING
          if(shortfile) then
             rhopos=10; vpos=12; bpos=16; ppos=21; gridpos=4; Bepos=22
          else
             rhopos=10; Bepos=32; ppos=33; vpos=12; bpos=25; gdetpos=16; gridpos=4; nthpos=43 !AC: is nthpos correct?? change for 3D?
          endif

          write(6,*) 'dfile: ',dfile
          write(6,*) 'data file: ',data_file
          open(unit=8,file=data_file,form='formatted',status='old',action='read')

          allocate(data(dlen,nx2)); nelem=nx2 !x2 is the fastest changing index
          allocate(grid(nx1*nx2*nx3,6));      !AC (x1,x2,x3,r,theta,phi)
          allocate(p(nx1*nx2*nx3)); allocate(rho(nx1*nx2*nx3)); allocate(Be(nx1*nx2*nx3))
          allocate(u(nx1*nx2*nx3)); allocate(b(nx1*nx2*nx3)); allocate(gdet(nx1*nx2*nx3))
          allocate(nnth(nx1*nx2*nx3, nrelbin));

          write(6,*) 'read koral sizes', dlen, nx1, nx2, nx3, size(data)
          read(8,*) header
          tcur=header(1)

          if(doingkoralnth) then
             read(8,*) header2
          end if

          deallocate(header)
          deallocate(header2)

          write(6,*) 'read koral past header: ', nx1, nx2, nx3, nelem, dlen

          !AC read the data
          do i=0, (nx1*nx2*nx3)/nelem-1
             read(8,*) data
             rho(1+i*nelem:(i+1)*nelem)=data(rhopos,:)
!             if(.not.shortfile) Be(1+i*nelem:(i+1)*nelem)=data(Bepos,:)
             Be(1+i*nelem:(i+1)*nelem)=data(Bepos,:)
             p(1+i*nelem:(i+1)*nelem)=data(ppos,:) !AC p is actually Te

             b(1+i*nelem:(i+1)*nelem)%data(1)=(data(bpos,:))
             b(1+i*nelem:(i+1)*nelem)%data(2)=(data(bpos+1,:))
             b(1+i*nelem:(i+1)*nelem)%data(3)=(data(bpos+2,:))
             b(1+i*nelem:(i+1)*nelem)%data(4)=(data(bpos+3,:))

             u(1+i*nelem:(i+1)*nelem)%data(1)=(data(vpos,:))
             u(1+i*nelem:(i+1)*nelem)%data(2)=(data(vpos+1,:))
             u(1+i*nelem:(i+1)*nelem)%data(3)=(data(vpos+2,:))
             u(1+i*nelem:(i+1)*nelem)%data(4)=(data(vpos+3,:))

             !AC grid is x1,x2,r,theta
             grid(1+i*nelem:(i+1)*nelem,:)=transpose(data(gridpos:gridpos+5,:))
             gdet(1+i*nelem:(i+1)*nelem)=(data(gdetpos,:))

             if(doingkoralnth) then
                nnth(1+i*nelem:(i+1)*nelem, : ) = transpose(data(nthpos:nthpos+nrelbin-1, : ))
             endif

          end do
          write(6,*) 'after read loop'
          close(unit=8)
          deallocate(data)
          x1_arr=grid(:,1)
          x2_arr=grid(:,2)
          x3_arr=grid(:,3)
          r_arr=grid(:,4)
          th_arr=grid(:,5)
          ph_arr=grid(:,6)
          write(6,*) 'read koral assign grid'
          deallocate(grid)

          ! Transform velocities, magnetic fields from MKS to KS and then BL:
          if(doingkoralnth) then
             write(6,*) 'read koral nnth 1 ', minval(nnth(:,1)), maxval(nnth(:,1))
          endif
          write(6,*) 'read koral grid sizes', size(x1_arr), size(x2_arr), size(x3_arr), size(r_arr), size(th_arr), size(ph_arr)
          write(6,*) 'read koral Te ', minval(p), maxval(p)
          write(6,*) 'read koral Be ', minval(Be), maxval(Be)
          write(6,*) 'read koral gdet ', minval(gdet), maxval(gdet)
          write(6,*) 'read koral rho ', minval(rho), maxval(rho)
          write(6,*) 'read koral b0', minval(b%data(1)), maxval(b%data(1))
          write(6,*) 'read koral u0 ', minval(u%data(1)), maxval(u%data(1))

          write(6,*) 'read koral transform coords r ', minval(r_arr), maxval(r_arr), asim
          write(6,*) 'read koral transform coords theta  ', minval(th_arr), maxval(th_arr), asim
          write(6,*) 'read koral transform coords phi',minval(x3_arr), maxval(x3_arr)

          write(6,*) 'WARNING: ignoring minpolecell=',minpolecell, ' cells on polar axis!!'

          deallocate(gdet)
        end subroutine read_koral3d_data_file

        subroutine initialize_koral3d_model(a,ifile,doingnth,df,hf,ntt,indft,sfac,nrelbin_bk,relgammamin_bk,relgammamax_bk)
        real(kind=8), intent(in) :: a
        integer,intent(out),optional :: nrelbin_bk
        real,intent(out),optional :: relgammamin_bk, relgammamax_bk
        integer :: nx, status, nhead
        logical, intent(in) :: doingnth
        real(8),intent(in),optional :: sfac
        character(len=20), intent(in), optional :: ifile
        character(len=20) :: default_ifile='koral.in' !AC
        character(len=100), intent(in), optional :: df,hf
        integer, intent(in), optional :: ntt,indft
        nhead=11

        doingkoralnth = doingnth

        ! if all inputs are given, assign variables rather than reading from file
        if (present(df)) then
           dfile = df
           hfile = hf
           nt = ntt
           indf = indft
        else
           if (present(ifile)) then
              call read_koral3d_inputs(ifile)
           else
              call read_koral3d_inputs(default_ifile)
           endif
        endif

        if (present(sfac)) then
            scalefac=sfac
        else
            scalefac=1.
        endif

        call read_koral3d_data_header(nhead)

        nrelbin_bk = nrelbin
        relgammamin_bk = relgammamin
        relgammamax_bk = relgammamax

        if (abs(asim-a).gt.1e-4) then
           write(6,*) 'ERROR -- Different simulation and grtrans spin values!'
           return
        endif
        n=nx1*nx2*nx3
        call init_koral3d_data(n,n*nt)
        call load_koral3d_data(nt)
        write(6,*) 'nrelbin',nrelbin
        end subroutine initialize_koral3d_model

        !Read data from file and then transfer to global arrays
        subroutine load_koral3d_data(nt)
        real(8), dimension(:), allocatable :: rho,p,Be
        real(8), dimension(:,:), allocatable :: nnth
        real, dimension(:), allocatable :: vrl, vtl, vpl
        type (four_vector), dimension(:), allocatable :: u, b
        integer, intent(in) :: nt
        character(len=250) :: append
        character(len=100) :: data_file
        integer :: k
        real(8) :: tcur, tstep_test

        allocate(nnth(n,nrelbin))
        allocate(rho(n)); allocate(p(n)); allocate(Be(n))
        allocate(vrl(n)); allocate(vtl(n)); allocate(vpl(n))
        allocate(u(n)); allocate(b(n))
        do k=1,nt
           !read data from  file
           write(append, fmt='(I4.4)') indf-(k-1)
           data_file = trim(dfile) // append
           call read_koral3d_data_file(data_file,tcur,rho,p,u,b,Be,nnth)
           t(k)=tcur

           ! transform velocities to LNRF frame
           call lnrf_frame(real(u%data(2)/u%data(1)),real(u%data(3)/u%data(1)), &
                           real(u%data(4)/u%data(1)), &
                           r_arr,asim,th_arr,vrl,vtl,vpl)
           write(6,*) 'lnrf transform', size(b0_arr), size(b%data(1)), size(b0_arr((k-1)*n+1:k*n))
           write(6,*) 'lnrf transform', size(vrl), size(vtl), size(vpl), size(u0_arr), size(nth_arr(:,1))

           ! now assign to data arrays
           write(6,*) scalefac, sqrt(scalefac)
           rho_arr((k-1)*n+1:k*n)=rho*scalefac
           p_arr((k-1)*n+1:k*n)=p !AC p_arr is actually the ELECTRON TEMPERATURE straight from koral -- no scale fac!
           Be_arr((k-1)*n+1:k*n)=Be
           b0_arr((k-1)*n+1:k*n)=b%data(1)*sqrt(scalefac)
           br_arr((k-1)*n+1:k*n)=b%data(2)*sqrt(scalefac)
           bth_arr((k-1)*n+1:k*n)=b%data(3)*sqrt(scalefac)
           bph_arr((k-1)*n+1:k*n)=b%data(4)*sqrt(scalefac)
           u0_arr((k-1)*n+1:k*n)=u%data(1)
           vrl_arr((k-1)*n+1:k*n)=real(vrl)
           vpl_arr((k-1)*n+1:k*n)=real(vpl)
           vtl_arr((k-1)*n+1:k*n)=real(vtl)

           if(doingkoralnth) then
              nth_arr((k-1)*n+1:k*n, :) = nnth*scalefac;
           endif
        end do

        ! Need to be able to do light curves even when nload (nt) = 1... Maybe read headers for first two files.
        if (nt.gt.1) then
           ! use default sim time step unless told otherwise
           tstep=t(1)-t(2)
           tstep_test=t(nt-1)-t(nt)
           if (abs((tstep-tstep_test)/tstep).gt.0.1) then
              write(6,*) 'WARNING -- Time step changes over dumps!'
           endif
        endif

        deallocate(rho)
        deallocate(p)
        deallocate(Be)
        deallocate(vrl); deallocate(vtl); deallocate(vpl)
        deallocate(u); deallocate(b)
        deallocate(nnth)
        end subroutine  load_koral3d_data

        ! Advance to next timestep file
        subroutine advance_koral3d_timestep(dtobs)
        real(8), intent(in) :: dtobs
        integer :: nupdate
        nupdate = floor(dtobs / tstep)
        toffset = toffset+dtobs-tstep*nupdate
        nupdate=nupdate+floor(toffset/tstep)
        toffset=toffset-floor(toffset/tstep)*tstep
        indf=indf+nupdate
        write(6,*) 'advance koral timestep: ',dtobs,tstep,toffset,indf,nupdate !AC
        call update_koral3d_data(nupdate)
        end subroutine advance_koral3d_timestep

        subroutine update_koral3d_data(nupdate)
        integer, intent(in) :: nupdate
        integer :: nshift
        if(nupdate.gt.nt) then ! this is case where we just load all new data
           call load_koral3d_data(nt)
        else ! this is case where we shift existing data back and then load more
           nshift=n*nupdate+1
           rho_arr(nshift:n*nt)=rho_arr(1:n*nt-nshift)
           p_arr(nshift:n*nt)=p_arr(1:n*nt-nshift)
           Be_arr(nshift:n*nt)=Be_arr(1:n*nt-nshift)
           br_arr(nshift:n*nt)=br_arr(1:n*nt-nshift)
           bth_arr(nshift:n*nt)=bth_arr(1:n*nt-nshift)
           bph_arr(nshift:n*nt)=bph_arr(1:n*nt-nshift)
           b0_arr(nshift:n*nt)=b0_arr(1:n*nt-nshift)
           vrl_arr(nshift:n*nt)=vrl_arr(1:n*nt-nshift)
           vtl_arr(nshift:n*nt)=vtl_arr(1:n*nt-nshift)
           vpl_arr(nshift:n*nt)=vpl_arr(1:n*nt-nshift)
           if(doingkoralnth) then
              nth_arr(nshift:n*nt,:)=nth_arr(1:n*nt-nshift,:)
           endif
           call load_koral3d_data(nupdate)
        end if
        end subroutine update_koral3d_data

        ! allocate arrays
        subroutine init_koral3d_data(nx,n)
          integer, intent(in) :: n,nx
          integer :: i
          real :: logbinspace, binspace
        allocate(rho_arr(n)); allocate(p_arr(n)); allocate(Be_arr(n)); allocate(u0_arr(n))
        allocate(vrl_arr(n)); allocate(vtl_arr(n)); allocate(vpl_arr(n))
        allocate(b0_arr(n)); allocate(br_arr(n)); allocate(bth_arr(n))
        allocate(bph_arr(n))
        allocate(x1_arr(nx)); allocate(x2_arr(nx));
        allocate(r_arr(nx)); allocate(th_arr(nx));
        allocate(x3_arr(nx)); allocate(ph_arr(nx));
        allocate(t(nt))
        if(doingkoralnth) then
           allocate(nth_arr(n,nrelbin))
        end if

        end subroutine init_koral3d_data

        ! deallocate arrays
        subroutine del_koral3d_data()
        deallocate(rho_arr); deallocate(p_arr); deallocate(Be_arr); deallocate(u0_arr)
        deallocate(vrl_arr); deallocate(vtl_arr); deallocate(vpl_arr)
        deallocate(b0_arr); deallocate(br_arr); deallocate(bth_arr)
        deallocate(bph_arr); deallocate(x1_arr); deallocate(x2_arr)
        deallocate(x3_arr); deallocate(ph_arr);
        deallocate(r_arr); deallocate(th_arr); deallocate(t)
        if(doingkoralnth) then
           deallocate(nth_arr)
        end if
        end subroutine del_koral3d_data

      end module fluid_model_koral3d

max_stars_repo_name: RAnantua/grtrans

max_issues_repo_name: RAnantua/grtrans

max_forks_repo_name: RAnantua/grtrans


text:
A great Japanesestyle restaurant which Departed Businesses used to be at Anderson Plaza.  Japanesestyle restaurant Sushi Unlimited now resides in its place. Known to some as Mongs for its excellent wiki:WikiPedia:Mongolian barbecue

20110810 16:24:41 nbsp I believe this was the first real Japanese place in town (old timers, correct me if Im wrong on this one).  They used to give Botan rice candy and finger puzzles to kids. Users/OldDavis73

20110810 16:25:57 nbsp The owner of Tokyo was the mother of the owners of Ding How. She was quite a character.  Users/DonShor


text:
!==============================================================================!
  subroutine Control_Mod_Monitoring_Points(verbose)
!------------------------------------------------------------------------------!
!   Reads monintoring points from control file.
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Monitor_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  logical, optional :: verbose
!-----------------------------------[Locals]-----------------------------------!
  integer           :: n
  real              :: xyz(3), def(3)
  character(len=80) :: point_name
!==============================================================================!

  call Control_Mod_Read_Int_Item('NUMBER_OF_MONITORING_POINTS', 0, &
                                  monitor % n_points, verbose)

  if(monitor % n_points .eq. 0) return

  ! Allocate memory for points
  allocate(monitor % x(monitor % n_points))
  allocate(monitor % y(monitor % n_points))
  allocate(monitor % z(monitor % n_points))

  do n = 1, monitor % n_points
    write(point_name, '(a,i3.3)') 'MONITORING_POINT_', n

    def = 0.  ! don't have a better idea what to set
    call Control_Mod_Read_Real_Array(point_name, 3, def, xyz, verbose)

    monitor % x(n) = xyz(1)
    monitor % y(n) = xyz(2)
    monitor % z(n) = xyz(3)
  end do

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


module swd_write_shape_3_def

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

implicit none
private

! Please adjust wp to the working precision applied in the wave generator
integer, parameter :: wp = kind(1.0d0) ! (double precision)

! This module provides a class for writing SWD files for shape 3.
!
! Written by Jens Bloch Helmers, August, 1. 2019
!
!------------------------------------------------------------------------------

!##############################################################################
!
!              B E G I N    P U B L I C    Q U A N T I T I E S
!
!------------------------------------------------------------------------------
!
public :: swd_write_shape_3
!
!------------------------------------------------------------------------------
!
!                E N D    P U B L I C    Q U A N T I T I E S
!
!##############################################################################

type :: swd_write_shape_3
    integer :: amp    ! Type of spectral amplitudes to be provided
    integer :: lu     ! Unit number associated with swd file
    integer :: nsteps ! Expected number of time steps
    integer :: it     ! Current time steps
    integer :: n      ! Number of spectral components in x-dir
    integer :: nh     ! Number of auxilary shapes (hat(n) in theory manual)
    integer :: nsf    ! Number of sea floor points
contains
    procedure :: init    ! Open the swd file and write the header section
    procedure :: update  ! Add spectral data at current time step
    procedure :: close   ! close the file
end type swd_write_shape_3

real(wp), parameter :: pi = 3.14159265358979323846264338327950288419716939937510582097494_wp

contains

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

subroutine init(self, file, prog, cid, grav, Lscale, amp, &
                n, nh, order, dk, dt, nsteps, isf, nsf, xsf, zsf)
class(swd_write_shape_3),  intent(out) :: self  ! Object to create
character(len=*),  intent(in)  :: file  ! Requested name for new SWD file
character(len=*),  intent(in)  :: prog  ! Name of actual wave generator including version number
character(len=*),  intent(in)  :: cid   ! Text describing actual input parameters to the
                                        ! wave generator. Hint: write(cid, nml=wave_input_namelist)
real(wp),          intent(in)  :: grav  ! Acceleration of gravity applied in wave generator
real(wp),          intent(in)  :: lscale! Number of length units per meter applied in wave generator
                                        ! E.g. if US feet was applied: Lscale = 1/0.3048 = 3.2808
integer,           intent(in)  :: amp   ! Flag to indicate type of provided spectral amplitudes
integer,           intent(in)  :: n     ! n as defined in the theory documentation
integer,           intent(in)  :: nh    ! hat(n) as defined in the theory documentation
integer,           intent(in)  :: order ! Expansion order applied in wave generator. (<0 for fully non-linear)
real(wp),          intent(in)  :: dk    ! Spacing of k-numbers in x-direction
real(wp),          intent(in)  :: dt    ! Spacing of time steps in the swd file
integer,           intent(in)  :: nsteps ! Total number of time steps to be stored in the swd file.
integer,           intent(in)  :: isf    ! Flag to indicate geometric description of the sea floor
                                         ! 0: Piecewise linear sea floor
integer,           intent(in)  :: nsf    ! Number of offset points describing the sea floor
real(wp),          intent(in)  :: xsf(:) ! X-positions of sea floor points
real(wp),          intent(in)  :: zsf(:) ! Z-positions of sea floor points
!
integer :: nid, i
integer, parameter :: fmt = 100
integer, parameter :: shp = 3
integer, parameter :: nstrip = 0
! We need some C compatible characters
character(kind=c_char, len=:), allocatable :: ccid
character(kind=c_char, len=30) :: cprog
character(kind=c_char, len=20) :: cdate

call input_sanity_checks

open( newunit=self % lu, file=file, access='stream', form='unformatted', &
      status='replace', action='write', convert='little_endian')
! The CONVERT option is an Intel, gfortran, HP and IBM extension.
! Modify this parameter for other compilers if your compiler complains.

!
self % amp = amp
self % nsteps = nsteps
self % it = 0
self % n = n
self % nh = nh
self % nsf = nsf
!
ccid = trim(cid) // c_null_char
nid = len(ccid)
cdate = timestamp() // c_null_char
cprog = trim(prog) // c_null_char
!
write(self % lu) 37.0221_c_float
write(self % lu) int(fmt, c_int)
write(self % lu) int(shp, c_int)
write(self % lu) int(amp, c_int)
write(self % lu) cprog
write(self % lu) cdate
write(self % lu) int(nid, c_int)
write(self % lu) ccid
write(self % lu) real(grav, c_float)
write(self % lu) real(lscale, c_float)
write(self % lu) int(nstrip, c_int)
write(self % lu) int(nsteps, c_int)
write(self % lu) real(dt, c_float)
write(self % lu) int(order, c_int)
write(self % lu) int(n, c_int)
write(self % lu) int(nh, c_int)
write(self % lu) real(dk, c_float)
write(self % lu) int(isf, c_int)
write(self % lu) int(nsf, c_int)
do i = 1, nsf
    write(self % lu) real(xsf(i), c_float)
end do
do i = 1, nsf
    write(self % lu) real(zsf(i), c_float)
end do

contains

  function timestamp() result(timestr)
    ! Return a string with the current time in the form YYYY-MM-DD hh-mm-ss
    character (len = 19) :: timestr

    integer :: y, m, d, h, n, s
    integer :: values(8)

    call date_and_time (values = values)

    y = values(1)
    m = values(2)
    d = values(3)
    h = values(5)
    n = values(6)
    s = values(7)

    write (timestr, '(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)' ) &
         y,'-',m,'-',d,' ',h,':',n,':',s

  end function timestamp
!
  subroutine input_sanity_checks
    if (nsf > 0) then
        ! Sea floor is completely submerged...
        if (any(zsf(1:nsf) >= 0.0_wp)) then
            print*, "ERROR: from swd_write_shape_3 :: init"
            print*, "All zsf(:) should be negative"
            print*, "zsf(1:nsf) = ", zsf(1:nsf)
            stop
        end if
        ! Sea floor x-range is in domain [0, 2pi/dk]
        if (abs(xsf(1)) > epsilon(xsf)) then
            print*, 'ERROR from swd_write_shape_3 :: init'
            print*, 'xsf(1) should be 0.0'
            print*, 'xsf(1) = ', xsf(1)
            stop
        end if
        if (abs(xsf(nsf) - 2*pi/dk) > 1.0e-4_wp * 2*pi/dk) then
            print*, 'ERROR from swd_write_shape_3 :: init'
            print*, 'xsf(nsf) should be 2*pi/dk'
            print*, 'xsf(nsf) = ', xsf(nsf)
            stop
        end if
        ! Sea floor x-values should be monotonic increasing
        do i = 2, nsf
            if (xsf(i-1) >= xsf(i)) then
                print*, 'ERROR from swd_write_shape_3 :: init'
                print*, 'xsf(1:nsf) should be monotonic increasing'
                print*, 'xsf(1:nsf) = ', xsf(1:nsf)
                stop
            end if
        end do
    end if
    if (nsf==0 .and. nh/=0) then
        print*, 'ERROR from swd_write_shape_3 :: init'
        print*, 'nh should be 0 if nsf is 0'
        print*, 'nsf, nh = ', nsf, nh
        stop
    end if
    if (nsf>0 .and. nh<=0) then
        print*, 'ERROR from swd_write_shape_3 :: init'
        print*, 'nh should be positive if nsf is positive'
        print*, 'nsf, nh = ', nsf, nh
        stop
    end if

  end subroutine input_sanity_checks

end subroutine init

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

subroutine update(self, h, ht, c, ct, ch, cht)
! Output of temporal functions as defined in the shape 3 class
class(swd_write_shape_3), intent(inout) :: self   ! Object to update
complex(wp),           intent(in) :: h(0:)   ! h(0:n) temporal amp
complex(wp),           intent(in) :: ht(0:)  ! dh/dt(0:n) temporal amp
complex(wp), optional, intent(in) :: c(0:)   ! c(0:n) temporal amp
complex(wp), optional, intent(in) :: ct(0:)  ! dc/dt(0:n) temporal amp
complex(wp), optional, intent(in) :: ch(0:)  ! hat(c)(0:nh) temporal amp
complex(wp), optional, intent(in) :: cht(0:) ! d hat(c)/dt(0:nh) temporal amp
!
self % it = self % it + 1
call dump(h, self % n)
call dump(ht, self % n)
if (self % amp < 3) then
    call dump(c, self % n)
    call dump(ct, self % n)
    if (self % nsf > 1) then
        call dump(ch, self % nh)
        call dump(cht, self % nh)
    end if
end if

contains

    subroutine dump(array, n)
        complex(wp), intent(in) :: array(:)
        integer,     intent(in) :: n
        integer :: j
        do j = 0, n
            write(self % lu) cmplx(array(j), kind=c_float)
        end do
    end subroutine dump
!
end subroutine update

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

subroutine close(self)
class(swd_write_shape_3), intent(inout) :: self   ! Object to update
!
if (self % it /= self % nsteps) then
    print*, "WARNING: from swd_write_shape_3 :: close"
    print*, "Specified number of time steps = ", self % nsteps
    print*, "Number of provided time steps = ", self % it
end if
!
close(self % lu)
!
end subroutine close

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

end module swd_write_shape_3_def