text:
C
C
      SUBROUTINE ANEOSV (NUMQQ,TQQ,RHOQQ,MATQQ,
     &  PQQ,EQQ,SQQ,CVQQ,DPDTQQ,DPDRQQ,FKROQQ,CSQQ,KPAQQ,R2PLQQ,R2PHQQ,
     &  ZBARQQ)
C
C***********************************************************************
C
C     ANEOS PACKAGE - VECTOR RUNNING ENTRY POINT
C
C     WARNING - depending on the computer system
C     All variables are single or double precision.
C
C     INPUTS (ALL UNITS CGS-EV)
C
C     NUMQQ =  dimensions of input and output arrays
C              (normally limited to 64)
C     TQQ   =  temperature array
C     RHOQQ =  density array
C     MATQQ =  material number array
C              skip any materials with matqq <= 0
C
C     Output arrays
C
C     PQQ   =  pressure
C     EQQ   =  specific energy
C     SQQ   =  specific entropy
C     CVQQ  =  specific heat capacity
C     DPDTQQ=  dP / dT
C     DPDRQQ=  dP / dRHO
C     FKROQQ=  Rosseland mean opacity (including conduction)
C     CSQQ  =  sound speed
C     KPAQQ =  phase indicator (see below)
C     R2PLQQ=  lower phase density for mixed-phase state
C              not used for single phase
C     R2PHQQ=  upper phase density for mixed-phase state
C              not used for single phase
C     ZBARQQ=  ionization number
C
C                                TABLE          ANEOS
C     KPAQQ=STATE INDICATOR      =1, 1p    =1, 1p    (eos without melt)
C                                =2, 2p lv =2, 2p liquid/solid plus vapor
C                                          =4, 1p solid  (eos with melt)
C                                          =5, 2p melt   (eos with melt)
C                                          =6, 1p liquid (eos with melt)
C                                =-1 bad value of temperature
C                                =-2 bad value of density
C                                =-3 bad value of material number
C
C                         Note: some versions stop for bad inputs
C                         instead of returning negative values of kpa.
C
C************************************** 6/89 version slt ***************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TQQ(NUMQQ),RHOQQ(NUMQQ),MATQQ(NUMQQ)
      DIMENSION PQQ(NUMQQ),EQQ(NUMQQ),SQQ(NUMQQ),CVQQ(NUMQQ)
      DIMENSION DPDTQQ(NUMQQ),DPDRQQ(NUMQQ),FKROQQ(NUMQQ),CSQQ(NUMQQ)
      DIMENSION KPAQQ(NUMQQ),R2PLQQ(NUMQQ),R2PHQQ(NUMQQ),ZBARQQ(NUMQQ)
      COMMON /FILEOS/ KLST, KINP
C     maximum number of aneos eos - must be 10 for csq
      PARAMETER (MAXMAT=10)
      COMMON /MIXCON/  RHOMPH(MAXMAT),RHOMPL(MAXMAT),IPEOS
      COMMON /ANES/  ACK(99*MAXMAT),ZZS(30*MAXMAT),COT(30*MAXMAT)
     1 ,FNI(30*MAXMAT),RCT(MAXMAT+1),TCT(MAXMAT+1),RSOL(100*MAXMAT)
     2 ,RVAP(100*MAXMAT),TTWO(100*MAXMAT),SAVER(92),BOLTS,EIP(4370)
     3 ,LOCSV(MAXMAT+1),LOCKP(MAXMAT+1),LOCKPL(MAXMAT+1),NOTRAD
      COMMON /MELTTRANS/ TMC(1000*MAXMAT),RSMC(1000*MAXMAT)
     1 ,RLMC(1000*MAXMAT),NUMTMC(MAXMAT),MELT_CURVE_TABLE
      COMMON /HPPTRANS/ THPP(1000*MAXMAT)
     2 ,RLPP(1000*MAXMAT),RHPP(1000*MAXMAT),NUMHPP(MAXMAT),NUMPHASE
      COMMON /ANE2/ CRHS2P,CRLS2P,ZBAR,FKROS
      COMMON /BNES/ PM,EM,SM,CVM,DPDTM,DPDRM
      COMMON /FNES/ CMLT7,ACK46
C     vector buffer size
      PARAMETER (MATBUF=64)
C     the following contains sqrt(t(i)) for vector aneos entry
C     ipsqts points to current value
      COMMON /ANESQT/ SQTS(MATBUF),IPSQTS
C     counters for aneos package
      PARAMETER (KDIMTR=15)
      COMMON /ANCNTR/ LNKNTR(MAXMAT+1,KDIMTR),LTKNTR,MATKTR
C     next three commons are here for multitask only
      COMMON /ANZB/ ZB(92), ZRAT
C                    inputs
      COMMON /ANEEL/ TEVX,RHOX,ABARX,ZBARM,T32X,FNX
C        outputs
     &  ,PE,EE,SE,CVE,DPTE,DPRE
C        inputs
     &  ,NMATSX,IIZX
      COMMON /ANZLTF/ ZLTOMF
      COMMON /TPMEOS/ NPORFE(MAXMAT)
      SAVE   /FILEOS/,/MIXCON/,/ANES/,/ANE2/,/BNES/,/FNES/,/ANCNTR/,
     & /ANZB/,/ANEEL/,/TPMEOS/
      DIMENSION KOC(MATBUF), PHOFG(MATBUF), NMATSV(MATBUF)
      PARAMETER (ZERO=0.D0)
      PARAMETER (HALF=0.5D0)
      PARAMETER (ONE=1.D0)
      PARAMETER (TWO=2.D0)
      PARAMETER (THREE=3.D0)
      PARAMETER (FOUR=4.D0)
      PARAMETER (FIVE=5.D0)
      PARAMETER (ATHIRD=1.D0/3.D0)
      PARAMETER (QCB1=1.D-20)
      PARAMETER (QCB2=1.D-10)
      PARAMETER (QCB4=-100.00001D0)
      PARAMETER (QCB5=8.D0)
      PARAMETER (QCB6=15.D0)
      PARAMETER (QCB7=24.D0)
      PARAMETER (QCB8=0.2D0)
      PARAMETER (CSFLAG=-77.66D11)
      PARAMETER (TQQMN=0.001D0)
      PARAMETER (RQQMN=1.D-22)
      CHARACTER*50 MESS
      DATA L8BAD /0/

C
      DO LMAT=1,NUMQQ
         KOC(LMAT)=-1
         PHOFG(LMAT)=ZERO
         NMATSV(LMAT)=-77
         IF(MATQQ(LMAT).GT.0) THEN
            KOC(LMAT)=LOCSV(MATQQ(LMAT))
            CSQQ(LMAT)=CSFLAG
            NMATSV(LMAT)=NINT(ACK(KOC(LMAT)+30))
            PHOFG(LMAT)=ACK(KOC(LMAT)+22)
            KPAQQ(LMAT)=1
         END IF
      END DO
      KUSER=0
      KIDG=0
      KFSQT=0
      KKSOL=0
C
      DO 410 LMAT=1,NUMQQ
      IF(MATQQ(LMAT).LE.0) GO TO 410
      IF((TQQ(LMAT).LE.ZERO).OR.(RHOQQ(LMAT).LE.ZERO)) GO TO 480
      MATKTR=MATQQ(LMAT)
      LNKNTR(MATKTR,2)=LNKNTR(MATKTR,2)+1
C
C     check for special forms that are not processed in this loop
C
      IF(NMATSV(LMAT).LE.-1) THEN
        IF(NMATSV(LMAT).EQ.-1) THEN
          KKSOL=1
        ELSEIF(NMATSV(LMAT).EQ.-4) THEN
          KUSER=1
          LNKNTR(MATKTR,15)=LNKNTR(MATKTR,15)+1
        ELSE
          KIDG=1
        END IF
        GO TO 410
      END IF
      IF(KFSQT.EQ.0) THEN
        DO 337 KMAT=1,NUMQQ
        SQTS(KMAT)=MAX(QCB2,TQQ(KMAT))
        SQTS(KMAT)=SQRT(SQTS(KMAT))
  337   CONTINUE
        KFSQT=1
      END IF
      IPSQTS=LMAT
      LOC=LOCSV(MATKTR)
      CMLT7=ACK(LOC+71)
      ACK46=ACK(LOC+46)
      FKROS=1.D10
      CRHS2P=RHOQQ(LMAT)
      CRLS2P=RHOQQ(LMAT)
    5 IF(NMATSV(LMAT)-2) 30,80,6
C
C     TEST FOR POROUS MATERIAL (OLD MODEL)
C
    6 IF(NPORFE(MATKTR)) 7,10,7
    7 IF((RHOQQ(LMAT).GE.ACK(LOC+19)).OR.
     &     (TQQ(LMAT).GE.TCT(MATKTR))) GO TO 30
      GO TO 20
C
C     CHECK FOR LIQUID-VAPOR OR SOLID-VAPOR STATE
C
   10 IF((RHOQQ(LMAT).GE.ACK(LOC+47)).OR.          ! If density is greater than rho_lm
     &    (TQQ(LMAT).GE.TCT(MATKTR))) GO TO 30     ! or temperature is above critical point then not two-phase
      IF (TQQ(LMAT).GT.ACK(LOC+18)) GO TO 20       ! If temperature is greater than melt temp., check for mixed phase
      IF (RHOQQ(LMAT).GE.ACK(LOC+23)) GO TO 30     ! If density is greater than minimum density skip two-phase test
   20 CALL ANE2PH (TQQ(LMAT),RHOQQ(LMAT),MATKTR,PQQ(LMAT),EQQ(LMAT),
     &  SQQ(LMAT),CVQQ(LMAT),DPDTQQ(LMAT),DPDRQQ(LMAT),LOC,KPAQQ(LMAT))
      LNKNTR(MATKTR,8)=LNKNTR(MATKTR,8)+1
C
C     KPA=2 IF LIQUID-VAPOR OR SOLID-VAPOR STATE
      IF (KPAQQ(LMAT).EQ.2) GO TO 100
C
C     IS MELT TRANSITION INCLUDED
C
 30   CONTINUE

      ! New code for finding melt-transition state via tables
      IF (MELT_CURVE_TABLE.EQ.1) THEN

         ! Check for LIQUID-SOLID STATE
         CALL ANE2PH_MELT (TQQ(LMAT),RHOQQ(LMAT),MATKTR,PQQ(LMAT),
     $        EQQ(LMAT),SQQ(LMAT),CVQQ(LMAT),DPDTQQ(LMAT),
     $        DPDRQQ(LMAT),LOC,KPAQQ(LMAT))

         IF (KPAQQ(LMAT).EQ.5) GOTO 100  ! Mixed phase

         IF (KPAQQ(LMAT).EQ.4) GOTO 50  ! Solid phase

         IF (KPAQQ(LMAT).EQ.6) GOTO 90  ! Liquid phase

         GOTO 80  ! single phase (must be vapour?)

      ELSEIF (NUMPHASE.GT.1 .AND. MATKTR.EQ.1) THEN  ! New LPP/HPP method not yet compatible with new melt transition method.

         ! Check for LPP-HPP STATE
         CALL ANE2PH_LPPHPP (TQQ(LMAT),RHOQQ(LMAT),MATKTR,PQQ(LMAT),
     $        EQQ(LMAT),SQQ(LMAT),CVQQ(LMAT),DPDTQQ(LMAT),
     $        DPDRQQ(LMAT),LOC,KPAQQ(LMAT))

         GOTO 100   ! This routine defines the state explicitly

      END IF

      IF (ACK46) 80,80,40
   40 KPAQQ(LMAT)=0
C     FATAL ERROR FLAG SET TO STOP IN ANELSM
      CALL ANELSM (TQQ(LMAT),RHOQQ(LMAT),DPDTQQ(LMAT),DPDRQQ(LMAT),LOC,
     &   KPAQQ(LMAT))
      LNKNTR(MATKTR,9)=LNKNTR(MATKTR,9)+1
      IF (KPAQQ(LMAT)-2) 50,70,60
C
C     SOLID STATE (EOS WITH MELT)
C
   50 KPAQQ(LMAT)=4
      CMLT7=-ONE  ! Turn off melt terms
      CALL ANEOS1 (TQQ(LMAT),RHOQQ(LMAT),PQQ(LMAT),EQQ(LMAT),SQQ(LMAT),
     &  CVQQ(LMAT),DPDTQQ(LMAT),DPDRQQ(LMAT),LOC)
      LNKNTR(MATKTR,7)=LNKNTR(MATKTR,7)+1
      CMLT7=ZERO  ! Turn back on melt terms
      GO TO 100
C
C     LIQUID STATE (EOS WITH MELT)
C
   60 KPAQQ(LMAT)=6
      GO TO 90
C
C     LIQUID-SOLID STATE (EOS WITH MELT)
C
   70 KPAQQ(LMAT)=5
      PQQ(LMAT)=PM
      EQQ(LMAT)=EM
      SQQ(LMAT)=SM
      CVQQ(LMAT)=CVM
      DPDTQQ(LMAT)=DPDTM
      DPDRQQ(LMAT)=DPDRM
      GO TO 100
C
C     ONE-PHASE STATE (EOS WITHOUT MELT)
C
   80 KPAQQ(LMAT)=1
   90 CALL ANEOS1 (TQQ(LMAT),RHOQQ(LMAT),PQQ(LMAT),EQQ(LMAT),SQQ(LMAT),
     &   CVQQ(LMAT),DPDTQQ(LMAT),DPDRQQ(LMAT),LOC)
      LNKNTR(MATKTR,7)=LNKNTR(MATKTR,7)+1
      GO TO 100
C
C-----------------------------------------------------------------------
C     error in call parameters
C
  480 PQQ(LMAT)=ZERO
      EQQ(LMAT)=ZERO
      SQQ(LMAT)=ZERO
      CVQQ(LMAT)=ONE
      PHOFG(LMAT)=ZERO
      DPDTQQ(LMAT)=ZERO
      DPDRQQ(LMAT)=ZERO
      FKROQQ(LMAT)=ONE
      R2PHQQ(LMAT)=ZERO
      R2PLQQ(LMAT)=ZERO
      ZBARQQ(LMAT)=ZERO
      CSQQ(LMAT)=ZERO
      IF(TQQ(LMAT).LE.ZERO) THEN
        KPAQQ(LMAT)=-1
      ELSE
        KPAQQ(LMAT)=-2
      END IF
      IF(L8BAD.EQ.0) THEN
        L8BAD=1
        CALL ANMARK
     1  ('ZERO OR NEGATIVE DENSITY OR TEMPERATURE IN ANEOS')
        MESS=
     1   'MATERIAL12345  T=1234567890123   RHO=1234567890123'
        WRITE(MESS(9:13),'(I5)',ERR=694) MATQQ(LMAT)
        WRITE(MESS(18:30),'(1PE13.6)',ERR=694) TQQ(LMAT)
        WRITE(MESS(38:50),'(1PE13.6)',ERR=694) RHOQQ(LMAT)
  694   CONTINUE
        CALL ANMARK(MESS)
      END IF
C     stop code for bad input to eos package
      STOP
C
C----------------------------------------------
C     set arrays for return
C
  100 FKROQQ(LMAT)=FKROS
      R2PHQQ(LMAT)=CRHS2P
      R2PLQQ(LMAT)=CRLS2P
      ZBARQQ(LMAT)=ZBAR
C
  410 CONTINUE
C
      IF(KIDG.EQ.1) THEN
C
C       IDEAL GAS EOS
C
        DO 422 LMAT=1,NUMQQ
        IF(NMATSV(LMAT).EQ.-2) THEN
          CVQQ(LMAT)=ACK(KOC(LMAT)+4)
          EQQ(LMAT)=CVQQ(LMAT)*TQQ(LMAT)
          DPDRQQ(LMAT)=ACK(KOC(LMAT)+5)*TQQ(LMAT)
          DPDTQQ(LMAT)=ACK(KOC(LMAT)+5)*RHOQQ(LMAT)
          PQQ(LMAT)=DPDTQQ(LMAT)*TQQ(LMAT)
          SQQ(LMAT)=CVQQ(LMAT)*(LOG(TQQ(LMAT))-ACK(KOC(LMAT)+15)*
     &              LOG(RHOQQ(LMAT))-ACK(KOC(LMAT)+6))
          FKROQQ(LMAT)=QCB8+ACK(KOC(LMAT)+20)/(ACK(KOC(LMAT)+25)+
     &              TQQ(LMAT))**3
          R2PHQQ(LMAT)=ZERO
          R2PLQQ(LMAT)=ZERO
          ZBARQQ(LMAT)=ZERO
        END IF
  422   CONTINUE
        CRHS2P=ZERO
        CRLS2P=ZERO
        ZBAR=ZERO
      END IF
C
      IF(KUSER.EQ.1) THEN
C
C       USER DEFINED EOS
C
        DO 424 LMAT=1,NUMQQ
        IF(NMATSV(LMAT).EQ.-4) THEN
          CALL ANUEOS(
C         inputs
     &    TQQ(LMAT),RHOQQ(LMAT),MATQQ(LMAT),ACK(KOC(LMAT)+1),
C         outputs
     &    PQQ(LMAT),EQQ(LMAT),SQQ(LMAT),CVQQ(LMAT),DPDTQQ(LMAT),
     &    DPDRQQ(LMAT),FKROQQ(LMAT),KPAQQ(LMAT),
     &    R2PLQQ(LMAT),R2PHQQ(LMAT),ZBARQQ(LMAT))
C
        END IF
  424   CONTINUE
        CRHS2P=R2PHQQ(1)
        CRLS2P=R2PLQQ(1)
        ZBAR=ZBARQQ(1)
        FKROS=FKROQQ(1)
        RHOMPH(IPEOS)=CRHS2P
        RHOMPL(IPEOS)=CRLS2P
        GO TO 428
      END IF
      IF(KKSOL.EQ.1) THEN
C
C     type -1 (Low temperature solid)
C
        CALL ANN1VL (NUMQQ,TQQ,RHOQQ,KOC,NMATSV,
     &  PQQ,EQQ,SQQ,CVQQ,DPDTQQ,DPDRQQ,FKROQQ,KPAQQ,R2PLQQ,R2PHQQ,
     &  ZBARQQ)
        CRHS2P=R2PHQQ(1)
        CRLS2P=R2PLQQ(1)
        ZBAR=ZBARQQ(1)
        FKROS=FKROQQ(1)
        RHOMPH(IPEOS)=CRHS2P
        RHOMPL(IPEOS)=CRLS2P
      END IF
C
C     phonon conduction term
C
      IF(NOTRAD.NE.0) THEN
        DO 320 LMAT=1,NUMQQ
        IF(PHOFG(LMAT).GT.ZERO) THEN
          PHOFG(LMAT)=ACK(KOC(LMAT)+22)*
     &     TQQ(LMAT)**(THREE-ACK(KOC(LMAT)+41))/RHOQQ(LMAT)
          FKROQQ(LMAT)=FKROQQ(LMAT)*PHOFG(LMAT)/
     &     (FKROQQ(LMAT)+PHOFG(LMAT))
        END IF
  320   CONTINUE
      END IF
C
C     old entry variables
C
      IF(NUMQQ.EQ.1) THEN
        RHOMPH(IPEOS)=CRHS2P
        RHOMPL(IPEOS)=CRLS2P
        FKROS=FKROQQ(1)
      END IF
C
C     sound speed, energy-entropy shifts
C
  428   DO 430 LMAT=1,NUMQQ
        IF(CSQQ(LMAT).EQ.CSFLAG) THEN
          CSQQ(LMAT)=SQRT(MAX(QCB1,(DPDRQQ(LMAT)
     &    +(TQQ(LMAT)*DPDTQQ(LMAT)**2)/(CVQQ(LMAT)*RHOQQ(LMAT)**2))))
          EQQ(LMAT)=EQQ(LMAT)+ACK(KOC(LMAT)+64)
          SQQ(LMAT)=SQQ(LMAT)+ACK(KOC(LMAT)+65)
        END IF
  430   CONTINUE
C
      RETURN
      END

max_stars_repo_name: chreinhardt/M-ANEOS

max_issues_repo_name: chreinhardt/M-ANEOS

max_forks_repo_name: chreinhardt/M-ANEOS


text:
module Orbitcombine
  use numeric_kinds
  implicit none
  private

  integer (kind=i4b),private,dimension(:,:),allocatable :: orbinfo
  real (kind=dp)    ,dimension(:),allocatable :: orbwght
  character (len=256) :: outroot
! setup/run/stop the program.
public :: setup,run,stob

private :: read_orbsol
contains

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine setup()
    use integrator         , only : integrator_setup
    use projection         , only : projection_setup,projection_change_direction
    use quadrantgrid       , only : qgrid_setup
    use aperture_routines  , only : aperture_setup
    use histograms         , only : histogram_setup
    use psf                , only : psf_setup
    use output             , only : output_setup
  !----------------------------------------------------------------------
    character(len = 80) :: string
    print *,"  ** Start Orbit combine Setup"
    print *,"  * Give setup version info: [U for unspecified]"
    read *, string
    if (string == "#Orbitcombine_setupfile_version_2" .or. string=="U") then
       print*,"  * Setupfile is Version 1"
       call integrator_setup()
       call projection_setup()
       !call projection_change_direction()
       !call qgrid_setup()
       !call psf_setup()
       !call aperture_setup()
       !call histogram_setup()
       !call output_setup()
	   print*,'bla'
       call read_orbsol()
    else
       print*, "This version is not understood by this program"
       STOP "program terminated in high_level:setup"
    end if

    print *,"  ** Setup Finished"

  end subroutine setup

subroutine read_orbsol ()
  use integrator  , only : integrator_points,integrator_dithering,rcirc
  use initial_parameters, only : totalmass,conversion_factor
  character (len=30) :: infil
  integer (kind=i4b)  :: i,j,norb,nparticles
  real    (kind=dp )  :: precision,t1,t2,rem,t3
  real (kind=dp)    ,dimension(:),allocatable :: massscale

  print*, "  * Give name of file with the orbitsolution (*_orb.out)"
  read (unit=*, fmt="(a256)") infil
  print*,"   ",infil

  open (unit=40,file=infil, status="OLD", action="read",position="rewind")
  read (unit=40, fmt=*) norb
  allocate(orbinfo(6,norb),orbwght(norb),massscale(norb))
  do i=1,norb
     !read (unit=40, fmt="(6i8,es13.5)") j, orbinfo(1,i),orbinfo(2,i),orbinfo(3,i),orbinfo(4,i),orbinfo(5,i),orbwght(i)
     read (unit=40, fmt=*) j, orbinfo(1,i),orbinfo(2,i),orbinfo(3,i),orbinfo(4,i),orbinfo(5,i),orbwght(i)
     !orbint(1:4,i),orbtype(i),orbweight(i)
  end do
  close (unit=40)

   i=LEN_TRIM(infil)
   if (i-4 < 1) stop "output filename too short"
   outroot =      infil(1:i-4)

  print*, "  * max amount of particles that could be generated (int64):", huge(1_i4b)
  print*, "  * Give the number of particles to generate:"
  read *, nparticles
  print*,"   ",nparticles

  if (nparticles <= 1_i4b) stop

  print*,size(rcirc) ,norb/3 ,(norb/3*(integrator_dithering**3))
  !massscale(1:norb/3)= &
  !  rcirc(125/2:(norb/3*(integrator_dithering**3)):integrator_dithering**3) &
  !  /conversion_factor
  !massscale(Norb/3  +1:Norb/3*2)=massscale(1:norb/3)
  !massscale(Norb/3*2+1:Norb)=massscale(1:norb/3)
  !massscale(:)=(massscale(:)+0.01)/(massscale(:)+100)
  !t1=sum(massscale(:))

  print*,'Equal mass particles'
  massscale(:)=sum(orbwght(:)) !1.0_dp
  !t1=1.0_dp
  t1=sum(orbwght(:))

 print*,"  * Total orbweight (~1.00): ",sum(orbwght(:))
  ! Nint rounds to nearest whole number
  !orbinfo(6,:)=nint(orbwght(:)*nparticles/massscale(:)*t1)
  orbinfo(6,:)=nint(orbwght(:)*nparticles/massscale(:)*t1)

  ! Enforce exactly equal mass particles, and add the remainder mass to the
  ! future orbits.
  rem=0.0_dp
  do i=1,norb
	 t3= orbinfo(6,i) * massscale(i) / nparticles /t1
     rem = rem + orbwght(i) - t3
     !if (rem .ge. 4.0_dp/nparticles/t1*massscale(i)) then
	 !  rem = rem - 1.0_dp/nparticles/t1*massscale(i)
	 !  orbinfo(6,i) = orbinfo(6,i) + 1
	 !print*,i,rem,1.0_dp/nparticles
	 !endif
     orbwght(i) = orbinfo(6,i) * massscale(i) / nparticles /t1
  end do
  print*,'remaining mass remainder',rem

  ! ensure all the mass in the particles
  !where ( orbinfo(6,:) .lt. 1 .and. orbwght(:) .gt. 1.0e-20)    orbinfo(6,:) = 1

  t2=0.0_dp
  do i=1,norb
     if ( orbinfo(6,i) .gt. 0 )  t2=t2+orbwght(i)
  end do

  !print*,massscale/conversion_factor
  print*,conversion_factor
  t1=nparticles - sum (orbinfo(6,:))

  print*,"  * Total orbweight (~1.00): ",sum(orbwght(:))
  print*,"  * Total mass in realization: ",t2
  print*,"  * relative Mass percision of the model: ",(t2-sum(orbwght))/sum(orbwght)
  print*,"  * absolute Mass percision of the model: ",(t2-sum(orbwght))*totalmass
  print*,"  * Total number of photons:", sum(dble(orbinfo(6,:)))
  print*,"  * Relative Precision of the model :",(t1)/nparticles
  print*,"  * Absolute number of photons difference:",t1
  print*,"  * symmetry number:",(integrator_dithering**3* 8.0)

  open (unit=33,file=trim(outroot)//"_part_masses.out",action="write", status="replace")
  do i=1, norb
  write (unit=33, fmt=*) orbwght(i), massscale(i), orbinfo(6,i)
  end do
  !stop
end subroutine read_orbsol

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine run()
    use histograms  , only : histogram_reset,hist_thesame, &
         histogram_velbin,histogram_store
    use projection  , only : project,projection_symmetry
    use integrator  , only : integrator_integrate,integrator_points,vyini,&
         integrator_current,vyini,ini_integ,nEner,nI2,nI3,integrator_dithering
    use output      , only : output_write
    use quadrantgrid, only : qgrid_reset,qgrid_store
    use psf         , only : psf_n,psf_gaussian
    use aperture         , only : aperture_n,aperture_psf
    use aperture_routines, only : aperture_find
    use initial_parameters, only : xmbh,grav_const_km, conversion_factor
    !use output , only: out_file
    !----------------------------------------------------------------------
    logical :: done,first,alldone
    real (kind=dp)    ,dimension(:,:),allocatable :: pos,vel
    real (kind=dp)    ,dimension(:,:),allocatable :: proj,vec_gauss
    real (kind=dp)    ,dimension(:  ),allocatable :: losvel
    integer (kind=i4b),dimension(:  ),allocatable :: velb,poly
    integer (kind=i4b)                            :: ap,i,pf,no
    integer (kind=i4b)  :: type
    real    (kind=dp )  :: t1,t2

	print*," writing particles file: ",trim(outroot)//"_particles.out"
	open (unit=33,file=trim(outroot)//"_particles.out",action="write", status="replace")

    ! Define units
    write (unit=33, fmt="(A40,1es20.12)"),"Gravitional constant in km^3/(s^2 Msun):",   grav_const_km
    ! Parsec unit is not needed for N-body calculations. We only use it for MGE->pot
    !write (unit=33, fmt="(A20,1es20.12)"),"Parsec in km",   parsec_km
    write (unit=33, fmt="(A40,1es20.12)"),"To convert from km to arcsec use:",1.0_dp/conversion_factor

    write (unit=33,fmt="(A40)"), "First particle is the black hole"
    write (unit=33,fmt="(A40)"), "Rows: Mass, x, y, z, vx, vy, vz"
    write (unit=33,fmt="(A40)"), "Units: Msun, km, km, km, km/s, km/s, km/s"
    ! Black hole is the first particle
    write (unit=33, fmt="(7es20.12)") xmbh,0.0,0.0,0.0,0.0,0.0,0.0

	  no=(nEner*nI2*nI3/integrator_dithering**3)
    print*,no
    call compute_orbits(1,2,no*2)
    print*,"  * Starting second set of orbits"
    vyini(:)=-vyini(:)
    call compute_orbits(2,2,no*2)

    print*,"  * Starting box orbits"
    call ini_integ()
    call compute_orbits(no*2+1,1,no*3)

    write (unit=33,fmt="(A20)"), "The END_OF_FILE"
    close(unit=33)

    !call output_write()

    print*,"  ** Finished Orbit Calculations"

  end subroutine run

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine compute_orbits(orbstart,orbdelta,orbend)
    use histograms  , only : histogram_reset,hist_thesame, &
         histogram_velbin,histogram_store
    use projection  , only : project,projection_symmetry
    use integrator  , only : integrator_integrate,integrator_points,vyini,&
         integrator_current,nEner,integrator_set_current,integrator_dithering
    use output      , only : output_write
    use quadrantgrid, only : qgrid_reset,qgrid_store
    use psf         , only : psf_n,psf_gaussian
    use aperture         , only : aperture_n,aperture_psf
    use aperture_routines, only : aperture_find
    use initial_parameters, only : totalmass
    !----------------------------------------------------------------------
    integer (kind=i4b),intent(in) :: orbstart,orbdelta,orbend
    logical :: done,first,alldone
    real (kind=dp)    ,dimension(:,:),allocatable :: pos,vel
    real (kind=dp)    ,dimension(:,:),allocatable :: postot,veltot
    real (kind=dp)    ,dimension(:,:),allocatable :: proj,vec_gauss
    real (kind=dp)    ,dimension(:  ),allocatable :: losvel
    integer (kind=i4b),dimension(:  ),allocatable :: velb,poly
    integer (kind=i4b)                            :: ap,i,ri,pf,orbn,nphot ,j
    integer (kind=i4b)  :: type, indtot,index_project
    real    (kind=dp )  :: t1,t2,rand,weight

    orbn=orbstart
    integrator_current=0
    alldone=.false.
    print*,"  ** Starting Orbit Calculations"
    allocate (postot(integrator_points*8*integrator_dithering**3,3))
	allocate (veltot(integrator_points*8*integrator_dithering**3,3))
    allocate (pos(integrator_points,3),vel(integrator_points,3))
    allocate (proj(integrator_points*projection_symmetry,2))
    allocate (vec_gauss(integrator_points*projection_symmetry,2))
    allocate (losvel(integrator_points*projection_symmetry))
    allocate (velb(integrator_points*projection_symmetry))
    allocate (poly(integrator_points*projection_symmetry))

    do  ! for each orbit

       call cpu_time(t1)
       indtot=0
       nphot=orbinfo(6,orbn)
       weight=orbwght(orbn)*totalmass
       orbn=orbn+orbdelta
       integrator_current=integrator_current+1
       if (nphot > 0) then
       	print*," This",(orbn-orbstart)/orbdelta," orbit has ",nphot,"photons"
	    !print*,orbn,integrator_current

	   first=.true.
       call integrator_set_current(integrator_current-1)
	   index_project=1_i4b
		postot=0.0_dp
		veltot=0.0_dp
		j=1
       do ! for all dithers
          call integrator_integrate(pos,vel,type,done,first,alldone)
          first=.false.
          if (done .or. alldone) exit
          if (j .eq. integrator_dithering**3/2+1) then
           do i=1,integrator_points
          write (unit=33, fmt="(7es20.12)") 1.0,pos(i,:),vel(i,:)
           enddo
          endif
		! project overwrite the same vector everytime
!          call project_part(type,pos,vel,postot,veltot,index_project)
           j=j+1
        end do

!	   ! Fisher-Yates shuffle
!       do i=1,nphot
!		   call random_number(rand)
!		   ! random number between 1 and integrator_points*8*integrator_dithering**3- (i-1)
!		   ri=int(rand*(integrator_points*8*integrator_dithering**3- (i-1))) + 1
!		   write (unit=33, fmt="(7es20.12)") weight/nphot,postot(ri,:),veltot(ri,:)
!		   ! Swap elements to make ensure they do not get drawn again.
!		   postot(size(postot,1)-(i-1),:)=postot(ri,:)
!		   veltot(size(postot,1)-(i-1),:)=veltot(ri,:)
!		end do

       call cpu_time(t2)
       print*,"  * Time spent one orbit:",t2-t1," seconds"
       end if
       if (alldone .or. orbn > orbend) exit
    end do
    deallocate (postot,veltot,pos,vel,proj,vec_gauss,losvel,velb,poly)
    print*,"  ** Finished Orbit Calculations"

  end subroutine compute_orbits

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine project_part(type,pos,vel,postot,veltot,ind)
  use initial_parameters, only : theta_view, phi_view
    real (kind=dp),intent(in ),dimension(:        ,:)   :: pos
    real (kind=dp),intent(in ),dimension(size(pos,1),3) :: vel
    real (kind=dp),intent(out),dimension(:,:)           :: postot
    real (kind=dp),intent(out),dimension(:,:)           :: veltot
    integer (kind=i4b),intent(in)                       :: type
	integer (kind=i4b),intent(inout)                    :: ind
  !----------------------------------------------------------------------
    real (kind=dp)               :: t1, t2, t3,theta,phi
    integer (kind=i4b)          :: i ,k ,l
    ! Signs of the (vx,vy,vz) for each Projection and type of Orbit
!    real (kind=dp),dimension(3,8,5),parameter :: vsgn= reshape((/  &
!    ! X tubes
!    1 , 1 , 1    ,-1 , 1 , 1  , -1 , 1 ,-1  ,  1 , 1 , -1 , &
!    1 ,-1 , 1    ,-1 ,-1 , 1  , -1 ,-1 ,-1  ,  1 ,-1 , -1 , &
!    ! Y tubes
!    1 , 1 , 1    , 1 , 1 ,-1  ,  1 ,-1 ,-1  ,  1 ,-1 , 1 , &
!   -1 , 1 , 1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  , -1 ,-1 , 1 , &
!    ! Z tubes
!    1 , 1 , 1    , 1 ,-1 , 1  , -1 ,-1 , 1  , -1 , 1 , 1 , &
!    1 , 1 ,-1    , 1 ,-1 ,-1  , -1 ,-1 ,-1  , -1 , 1 ,-1 , &
!    ! Boxed
!    1 , 1 , 1    ,-1 , 1 , 1  , -1 ,-1 , 1  ,  1 ,-1 , 1 , &
!    1 , 1 ,-1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  ,  1 ,-1 ,-1 , &
!    ! Stochastic
!    1 , 1 , 1    ,-1 , 1 , 1  , -1 ,-1 , 1  ,  1 ,-1 , 1 , &
!    1 , 1 ,-1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  ,  1 ,-1 ,-1 /),(/3,8,5/))
    real (kind=dp),dimension(3,8,5),parameter :: vsgn= reshape((/  &
    ! X tubes
    1 , 1 , 1    ,-1 , 1 , 1  ,  1 , 1 ,-1  , -1 , 1 , -1 , &
   -1, -1 , 1    , 1 ,-1 , 1  , -1 ,-1 ,-1  ,  1 ,-1 , -1 , &
    ! Y tubes
    1 , 1 , 1    , 1 ,-1 ,-1  ,  1 , 1 ,-1  ,  1 ,-1 , 1 , &
    -1,-1 , 1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  , -1 , 1 , 1 , &
    ! Z tubes
    1 , 1 , 1    , 1 ,-1 ,-1  , -1 ,-1 , 1  , -1 , 1 ,-1 , &
    1 , 1 ,-1    , 1 ,-1 , 1  , -1 ,-1 ,-1  , -1 , 1 , 1 , &
    ! Boxed
    1 , 1 , 1    ,-1 , 1 , 1  , -1 ,-1 , 1  ,  1 ,-1 , 1 , &
    1 , 1 ,-1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  ,  1 ,-1 ,-1 , &
    ! Stochastic
    1 , 1 , 1    ,-1 , 1 , 1  , -1 ,-1 , 1  ,  1 ,-1 , 1 , &
    1 , 1 ,-1    ,-1 , 1 ,-1  , -1 ,-1 ,-1  ,  1 ,-1 ,-1 /),(/3,8,5/))

    !Signs of the x,y,z for each projection  :psgn( [x,y,z], project )
    real (kind=dp),dimension(3,8),parameter :: psgn= reshape((/  &
    1 , 1 , 1   , -1 , 1 , 1  , -1 , -1 , 1 ,  1 , -1 , 1 , &
    1 , 1 ,-1   , -1 , 1 ,-1  , -1 , -1 ,-1 ,  1 , -1 ,-1 /),(/3,8/))

    ! check orbit type
    if (type > 5 .or. type < 1 ) stop "project_n: Wrong orbit type"

    ! Use sign matrix for the symmetries.
    ! Using the inverse (transpose) of the projection (eq. 4) of Thesis Ellen.

    do i=1,8
	  do l=1,3
	    postot(ind:ind+size(pos,1)-1,l)=pos(:,l)*psgn(l,i)
	    veltot(ind:ind+size(pos,1)-1,l)=vel(:,l)*vsgn(l,i,type)
      end do
      ind=ind+size(pos,1)
	end do
    ! y'
    !t1 = -cos(theta)*cos(phi) * psgn(1,n)
    !t2 = -cos(theta)*sin(phi) * psgn(2,n)
    !t3 =  sin(theta)          * psgn(3,n)
    !proj(:,2) = t1 * pos(:,1) + t2 * pos(:,2) + t3 * pos(:,3)

    ! v_LOS
    !t1 = sin(theta)*cos(phi)  * vsgn(1,n,type)
    !t2 = sin(theta)*sin(phi)  * vsgn(2,n,type)
    !t3 = cos(theta)           * vsgn(3,n,type)
    !losvel(:) = t1 * vel(:,1) + t2 * vel(:,2) + t3 * vel(:,3)

  end subroutine project_part

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  subroutine stob()
    use integrator        , only : integrator_stop
    use projection        , only : projection_stop
    use histograms        , only : histogram_stop
    use psf               , only : psf_stop
    use aperture_routines , only : aperture_stop
    use output            , only : output_close
  !----------------------------------------------------------------------

    call output_close    ()
    call integrator_stop ()
    call projection_stop ()
    call aperture_stop   ()
    call psf_stop        ()
    call histogram_stop  ()

    deallocate(orbinfo)

  end subroutine stob

end module Orbitcombine

!######################################################################
!######################################################################
!######################################################################

program counterrotation
  use numeric_kinds
  use Orbitcombine, only : setup,run,stob
  implicit none
  real (kind=dp) :: t1,t2
  print*,"  ** Triaxial Orbit library by Remco C.E. van den Bosch <bosch@strw.leidenuniv.nl>"
  print*,"  * Jan 2007 "
  print*,"  * $Id: partgen.f90,v 1.1 2010/03/17 12:56:57 bosch Exp $"

  call setup()
  print*,"Starting orbit integration routines"
  call cpu_time(t1)
  call run()
  call cpu_time(t2)
  call stob()
  print*,"  ** Program Finished"
  print*,"  * Time spent calculating :",t2-t1," seconds"
end program counterrotation

!######################################################################
!######################################################################
!######################################################################

max_stars_repo_name: dynamics-of-stellar-systems/dynamite

max_issues_repo_name: dynamics-of-stellar-systems/dynamite

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


module fields

  implicit none

  private

  public :: init_fields, finish_fields
  public :: read_parameters, wnml_fields, check_fields
  public :: advance, force_maxwell_reinit
  public :: reset_init, set_init_fields
  public :: fields_init_response, set_dump_and_read_response
  public :: dump_response_to_file
  public :: init_fields_level_1, init_fields_level_2
  public :: finish_fields_level_1, finish_fields_level_2

  !> Made public for unit tests
  public :: fields_pre_init
  public :: remove_zonal_flows_switch
  !> Made public for replay
  public :: allocate_arrays

  interface fieldlineavgphi
     module procedure fieldlineavgphi_loc
     module procedure fieldlineavgphi_tot
  end interface

  ! knobs
  integer :: fieldopt_switch
  logical :: remove_zonal_flows_switch
  logical :: force_maxwell_reinit
  integer, parameter :: fieldopt_implicit = 1, fieldopt_test = 2, fieldopt_local = 3
  logical :: dump_response, read_response
  logical :: initialized = .false.
  logical :: exist

contains

  subroutine check_fields(report_unit)
    use fields_local, only: minNrow, do_smart_update
    implicit none
    integer, intent(in) :: report_unit
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       write (report_unit, fmt="('The field equations will be advanced in time implicitly.')")
       if(dump_response) write (report_unit, fmt="('The response matrix will be dumped to file.')")
       if(read_response) write (report_unit, fmt="('The response matrix will be read from file.')")
    case (fieldopt_test)
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('The field equations will only be tested.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    case (fieldopt_local)
       write (report_unit, fmt="('The field equations will be advanced in time implicitly with decomposition respecting g_lo layout.')")
       if(dump_response) write (report_unit, fmt="('The response matrix will be dumped to file.')")
       if(read_response) write (report_unit, fmt="('The response matrix will be read from file.')")
       write(report_unit, fmt="('Using a min block size of ',I0)") minNrow
       if(do_smart_update) write(report_unit, fmt="('Using optimised field update.')")
    end select
  end subroutine check_fields

  subroutine wnml_fields(unit)
    use fields_local, only: minNrow, do_smart_update
    implicit none
    integer, intent(in) :: unit
    if (.not. exist) return
    write (unit, *)
    write (unit, fmt="(' &',a)") "fields_knobs"
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       write (unit, fmt="(' field_option = ',a)") '"implicit"'
    case (fieldopt_test)
       write (unit, fmt="(' field_option = ',a)") '"test"'
    case (fieldopt_local)
       write (unit, fmt="(' field_option = ',a)") '"local"'
       write (unit, fmt="(' minNrow = ',I0)") minNrow
       write (unit, fmt="(' do_smart_update = ',L1)") do_smart_update
    end select
    if(dump_response) write (unit, fmt="(' dump_response = ',L1)") dump_response
    if(read_response) write (unit, fmt="(' read_response = ',L1)") read_response
    write (unit, fmt="(' /')")
  end subroutine wnml_fields

  !> Calls all initialisations required for init_fields_implicit/local,
  !! reads parameters and allocates field arrays
  subroutine fields_pre_init
    use theta_grid, only: init_theta_grid
    use run_parameters, only: init_run_parameters
    use dist_fn, only: init_dist_fn
    use init_g, only: ginit, init_init_g
    use antenna, only: init_antenna
    use unit_tests, only: debug_message
    implicit none
    integer, parameter :: verb=3

    call debug_message(verb, "init_fields: init_theta_grid")
    call init_theta_grid

!CMR,30/3/2009:
! call init_init_g before init_run_parameters to read delt from restart file

    call debug_message(verb, "init_fields: init_init_g")
    call init_init_g
    call debug_message(verb, "init_fields: init_run_parameters")
    call init_run_parameters
    call debug_message(verb, "init_fields: init_dist_fn")
    call init_dist_fn
    !call debug_message(verb, "init_fields: init_parameter_scan")
    !call init_parameter_scan
    call debug_message(verb, "init_fields: init_antenna")
    call init_antenna !Must come before allocate_arrays so we know if we need apar_ext
    !call debug_message(verb, "init_fields: read_parameters")
    !call read_parameters
    !call debug_message(verb, "init_fields: allocate_arrays")
    !call allocate_arrays
  end subroutine fields_pre_init

  subroutine init_fields_level_1
    use unit_tests, only: debug_message
    implicit none
    integer, parameter :: verb=3
    call debug_message(verb, "init_fields: read_parameters")
    call read_parameters
    call debug_message(verb, "init_fields: allocate_arrays")
    call allocate_arrays
  end subroutine init_fields_level_1

  subroutine finish_fields_level_1
    call finish_fields
  end subroutine finish_fields_level_1

  subroutine init_fields_level_2
    call init_fields
  end subroutine init_fields_level_2

  subroutine finish_fields_level_2
    call reset_init
  end subroutine finish_fields_level_2

  subroutine fields_init_response
    use fields_implicit, only: init_fields_implicit
    use fields_test, only: init_fields_test
    use fields_local, only: init_fields_local
    use unit_tests, only: debug_message
    implicit none
    integer, parameter :: verb=3
    logical, parameter :: debug = .false.
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       call debug_message(verb, &
         "fields::fields_init_response init_fields_implicit")
       call init_fields_implicit
    case (fieldopt_test)
       call debug_message(verb, "fields::fields_init_response init_fields_test")
       call init_fields_test
    case (fieldopt_local)
       call debug_message(verb, &
         "fields::fields_init_response init_fields_local")
       call init_fields_local
    case default
       !Silently ignore unsupported field options
    end select
  end subroutine fields_init_response

  subroutine init_fields
!CMR,18/2/2011:
! add optional logical arg "noalloc" to avoid array allocations in ingen
    use theta_grid, only: init_theta_grid
    use run_parameters, only: init_run_parameters
    use dist_fn, only: init_dist_fn
    use init_g, only: ginit, init_init_g
    use nonlinear_terms, only: nl_finish_init => finish_init
    use antenna, only: init_antenna
    use kt_grids, only: gridopt_switch, gridopt_box, kwork_filter
    implicit none
    logical :: restarted
    logical, parameter :: debug=.false.

    if (initialized) return
    initialized = .true.

    call fields_pre_init

    call fields_init_response

! Turn on nonlinear terms.
    if (debug) write(6,*) "init_fields: nl_finish_init"
    call nl_finish_init

    ! EGH Commented out the following lines as they are now
    ! handled by gs2_init

    !if (debug) write(6,*) "init_fields: ginit"
    !call ginit (restarted)
    !if (restarted .and. .not. force_maxwell_reinit) return
    !if (debug) write(6,*) "init_fields: init_antenna"
    !call init_antenna

    !!Set the initial fields
    !call set_init_fields

    !If running in flux tube disable evolution of ky=kx=0 mode
    if(gridopt_switch.eq.gridopt_box) kwork_filter(1,1)=.true.
  end subroutine init_fields

  !>Force the current
  subroutine dump_response_to_file(suffix)
    use fields_implicit, only: dump_response_to_file_imp
    use fields_local, only: dump_response_to_file_local
    implicit none
    character(len=*), intent(in), optional :: suffix
    !Note can pass optional straight through as long as also optional
    !in called routine (and not different routines combined in interface)
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       call dump_response_to_file_imp(suffix)
    case (fieldopt_local)
       call dump_response_to_file_local(suffix)
    case default
       !Silently ignore unsupported field options
    end select
  end subroutine dump_response_to_file

  subroutine set_init_fields
    use fields_implicit, only: init_allfields_implicit
    use fields_test, only: init_phi_test
    use mp, only: proc0
    use fields_local, only: init_allfields_local
    implicit none
    logical, parameter :: debug=.false.
    if(proc0.and.debug) write(6,*) "Syncing fields with g."
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       if (debug) write(6,*) "init_fields: init_allfields_implicit"
       call init_allfields_implicit
    case (fieldopt_test)
       if (debug) write(6,*) "init_fields: init_phi_test"
       call init_phi_test
    case (fieldopt_local)
       if (debug) write(6,*) "init_fields: init_allfields_local"
       call init_allfields_local
    end select
  end subroutine set_init_fields

  subroutine read_parameters
    use file_utils, only: input_unit, error_unit, input_unit_exist
    use text_options, only: text_option, get_option_value
    use mp, only: proc0, broadcast
    use fields_implicit, only: field_subgath
    use fields_local, only: minNrow
    use fields_local, only: do_smart_update, field_local_allreduce, field_local_allreduce_sub
    use fields_arrays, only: response_file
    use file_utils, only: run_name
    implicit none
    type (text_option), dimension (5), parameter :: fieldopts = &
         (/ text_option('default', fieldopt_implicit), &
            text_option('implicit', fieldopt_implicit), &
            text_option('test', fieldopt_test),&
            text_option('local', fieldopt_local),&
            text_option('implicit_local', fieldopt_local)/)
    character(20) :: field_option
    character(len=256) :: response_dir
    namelist /fields_knobs/ field_option, remove_zonal_flows_switch, field_subgath, force_maxwell_reinit,&
         dump_response, read_response, minNrow, do_smart_update, field_local_allreduce, field_local_allreduce_sub, response_dir
    integer :: ierr, in_file

    if (proc0) then
       field_option = 'default'
       remove_zonal_flows_switch = .false.
       field_subgath=.false.
       force_maxwell_reinit=.true.
       dump_response=.false.
       read_response=.false.
       minnrow=64 !Tuning this can influence both init and advance times
       do_smart_update=.false.
       field_local_allreduce=.false.
       field_local_allreduce_sub=.false.
       response_dir=''
       in_file = input_unit_exist ("fields_knobs", exist)
       if (exist) read (unit=in_file, nml=fields_knobs)

       ierr = error_unit()
       call get_option_value &
            (field_option, fieldopts, fieldopt_switch, &
            ierr, "field_option in fields_knobs",.true.)

       if(trim(response_dir).eq.'')then
          write(response_file,'(A)') trim(run_name)
       else
          write(response_file,'(A,"/",A)') trim(response_dir),trim(run_name)
       endif

    end if

    call broadcast (fieldopt_switch)
    call broadcast (remove_zonal_flows_switch)
    call broadcast (field_subgath)
    call broadcast (force_maxwell_reinit)
    call broadcast (dump_response)
    call broadcast (read_response)

    !Setup response file location
    call broadcast(response_dir)
    call broadcast(response_file)

    !Set the solve type specific flags
    call set_dump_and_read_response(dump_response, read_response)
    select case (fieldopt_switch)
    case (fieldopt_implicit)
    case (fieldopt_test)
    case (fieldopt_local)
       call broadcast (minNrow)
       call broadcast (do_smart_update)
       call broadcast (field_local_allreduce)
       call broadcast (field_local_allreduce_sub)
    end select
  end subroutine read_parameters

  subroutine set_dump_and_read_response(dump_flag, read_flag)
    use fields_implicit, only: dump_response_imp => dump_response, read_response_imp=>read_response
    use fields_local, only: dump_response_loc => dump_response, read_response_loc=>read_response
    implicit none
    logical, intent(in) :: dump_flag, read_flag
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       dump_response_imp=dump_flag
       read_response_imp=read_flag
    case (fieldopt_local)
       dump_response_loc=dump_flag
       read_response_loc=read_flag
    case default
       !Silently ignore unsupported field types
    end select
  end subroutine set_dump_and_read_response

  subroutine allocate_arrays
    use theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0
    use antenna, only: no_driver
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew, apar_ext
    use unit_tests, only: debug_message
    implicit none
    integer, parameter :: verb=3

    if (.not. allocated(phi)) then
       call debug_message(verb, 'fields::allocate_arrays allocating')
       allocate (     phi (-ntgrid:ntgrid,ntheta0,naky))
       allocate (    apar (-ntgrid:ntgrid,ntheta0,naky))
       allocate (   bpar (-ntgrid:ntgrid,ntheta0,naky))
       allocate (  phinew (-ntgrid:ntgrid,ntheta0,naky))
       allocate ( aparnew (-ntgrid:ntgrid,ntheta0,naky))
       allocate (bparnew (-ntgrid:ntgrid,ntheta0,naky))
    endif
    phi = 0.; phinew = 0.
    apar = 0.; aparnew = 0.
    bpar = 0.; bparnew = 0.
    if(.not.allocated(apar_ext).and.(.not.no_driver))then
       allocate (apar_ext (-ntgrid:ntgrid,ntheta0,naky))
       apar_ext = 0.
    endif
  end subroutine allocate_arrays

  subroutine advance (istep)
    use fields_implicit, only: advance_implicit
    use fields_test, only: advance_test
    use fields_local, only: advance_local
    implicit none
    integer, intent (in) :: istep

    select case (fieldopt_switch)
    case (fieldopt_implicit)
       call advance_implicit (istep, remove_zonal_flows_switch)
    case (fieldopt_test)
       call advance_test (istep)
    case (fieldopt_local)
       call advance_local (istep, remove_zonal_flows_switch)
    end select
  end subroutine advance

  !This routine has a potentially misleading name and isn't used anywhere
  subroutine kperp (ntgrid_output, akperp)
    use theta_grid, only: delthet
    use kt_grids, only: naky, aky, ntheta0, kperp2
    use run_parameters, only: fphi, fapar, fbpar
    use fields_arrays, only: phinew, aparnew, bparnew
    implicit none
    integer, intent (in) :: ntgrid_output
    real, dimension (:,:), intent (out) :: akperp
    real :: anorm
    integer :: ik, it

    do ik = 1, naky
       do it = 1, ntheta0
          anorm = sum(abs(phinew(-ntgrid_output:ntgrid_output,it,ik)*fphi &
                         + aparnew(-ntgrid_output:ntgrid_output,it,ik)*fapar &
                         + bparnew(-ntgrid_output:ntgrid_output,it,ik)*fbpar)**2 &
                      *delthet(-ntgrid_output:ntgrid_output))
          if (anorm < 2.0*epsilon(0.0) .or. aky(ik) == 0.0) then
             akperp(it,ik) = 0.0
          else
             akperp(it,ik) &
                  = sqrt(sum(kperp2(-ntgrid_output:ntgrid_output,it,ik) &
                     *abs(phinew(-ntgrid_output:ntgrid_output,it,ik)*fphi &
                          + aparnew(-ntgrid_output:ntgrid_output,it,ik)*fapar &
                          + bparnew(-ntgrid_output:ntgrid_output,it,ik)*fbpar)**2 &
                     *delthet(-ntgrid_output:ntgrid_output))/anorm)
          end if
       end do
    end do
  end subroutine kperp

  !This routine isn't used anywhere
  subroutine fieldlineavgphi_loc (ntgrid_output, it, ik, phiavg)
    use theta_grid, only: ntgrid, drhodpsi, gradpar, bmag, delthet
    use fields_arrays, only: phi
    implicit none
    integer, intent (in) :: ntgrid_output, ik, it
    complex, intent (out) :: phiavg
    real, dimension (-ntgrid:ntgrid) :: jac

    jac = 1.0/abs(drhodpsi*gradpar*bmag)
    phiavg = sum(phi(-ntgrid_output:ntgrid_output,it,ik) &
                 *jac(-ntgrid_output:ntgrid_output) &
                 *delthet(-ntgrid_output:ntgrid_output)) &
            /sum(delthet(-ntgrid_output:ntgrid_output) &
                 *jac(-ntgrid_output:ntgrid_output))
  end subroutine fieldlineavgphi_loc

  !This doesn't look like a useful routine and isn't used anywhere
  subroutine fieldlineavgphi_tot (phiavg)
    use theta_grid, only: ntgrid, drhodpsi, gradpar, bmag
!    use theta_grid, only: delthet
    implicit none
    complex, dimension (:,:), intent (out) :: phiavg
    real, dimension (-ntgrid:ntgrid) :: jac
!    integer :: ntg

!    ntg = ntgrid
    jac = 1.0/abs(drhodpsi*gradpar*bmag)
!    phiavg = sum(phi(-ntg:ntg,:,:)*jac(-ntg:ntg)*delthet(-ntg:ntg)) &
!         /sum(delthet(-ntg:ntg)*jac(-ntg:ntg))
    phiavg = 0.
    write(*,*) 'error in fields'
  end subroutine fieldlineavgphi_tot

  !!> This generates a flux surface average of phi.

  !subroutine flux_surface_average_phi (phi_in, phi_average)
    !use theta_grid, only: ntgrid, drhodpsi, gradpar, bmag, delthet
    !use kt_grids, only: ntheta0, naky

    !implicit none
    !complex, intent (in) :: phi_in
    !complex, intent (out) :: phi_average
    !complex, dimension(-ntgrid:ntgrid,1:ntheta0,1:naky) :: phi_fieldline_avg
    !integer it, ig

    !call fieldline_average_phi(phi_in, phi_fieldline_avg)
    !do it = 1,ntheta0
      !do ig = -ntgrid,ntgrid
        !phi_average(ig, it, :) = sum(phi_fieldline_avg(ig, it, :))/real(naky)
      !end do
    !end do

  !end subroutine fieldline_average_phi

  subroutine reset_init
    use fields_implicit, only: fi_reset => reset_init
    use fields_test, only: ft_reset => reset_init
    use fields_local, only: fl_reset => reset_fields_local
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    implicit none
    initialized  = .false.
    phi = 0.
    phinew = 0.
    apar = 0.
    aparnew = 0.
    bpar = .0
    bparnew = 0.
    !What about apar_ext?
    select case (fieldopt_switch)
    case (fieldopt_implicit)
       call fi_reset
    case (fieldopt_test)
       call ft_reset
    case (fieldopt_local)
       call fl_reset
    end select
  end subroutine reset_init

  subroutine finish_fields

    use fields_implicit, only: implicit_reset => reset_init
    use fields_test, only: test_reset => reset_init
    use fields_local, only: finish_fields_local
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use fields_arrays, only: apar_ext
    use unit_tests, only: debug_message

    implicit none
    integer, parameter :: verbosity = 3

    initialized  = .false.
    call debug_message(verbosity, &
      'fields::finish_fields starting')
    phi = 0.
    phinew = 0.
    apar = 0.
    aparnew = 0.
    bpar = .0
    bparnew = 0.
    call debug_message(verbosity, &
      'fields::finish_fields zeroed fields')

    select case (fieldopt_switch)
    case (fieldopt_implicit)
       ! This line is no longer necessary
       ! because fields_implicit::reset_init is
       ! called by finish_fields_level_2
       !call implicit_reset
    case (fieldopt_test)
       ! This line is no longer necessary
       ! because fields_test::reset_init is
       ! called by finish_fields_level_2
       !call test_reset
    case (fieldopt_local)
       call finish_fields_local
    end select

    call debug_message(verbosity, &
      'fields::finish_fields called subroutines')

    if (allocated(phi)) deallocate (phi, apar, bpar, phinew, aparnew, bparnew)
    if (allocated(apar_ext)) deallocate (apar_ext)
    call debug_message(verbosity, &
      'fields::finish_fields deallocated fields')

  end subroutine finish_fields

end module fields


!
! SDF (Self-Describing Format) Fortran Library
! Copyright (c) 2011-2016, SDF Development Team
!
! Distributed under the terms of the BSD 3-clause License.
! See the LICENSE file for details.
!

MODULE sdf_input_r8

  USE sdf_input_ru

  IMPLICIT NONE

CONTAINS

  SUBROUTINE read_header_r8(h, step, time, code_name, code_io_version, &
      string_length, restart_flag, other_domains)

    TYPE(sdf_file_handle) :: h
    INTEGER, INTENT(OUT) :: step
    REAL(r8), INTENT(OUT) :: time
    CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: code_name
    INTEGER, INTENT(OUT), OPTIONAL :: code_io_version, string_length
    LOGICAL, INTENT(OUT), OPTIONAL :: restart_flag, other_domains

    CALL read_header_ru(h, step)

    time = REAL(h%time,r8)
    IF (PRESENT(code_io_version)) code_io_version = h%code_io_version
    IF (PRESENT(restart_flag)) restart_flag = h%restart_flag
    IF (PRESENT(other_domains)) other_domains = h%other_domains
    IF (PRESENT(code_name)) CALL sdf_safe_copy_string(h%code_name, code_name)
    IF (PRESENT(string_length)) string_length = h%string_length

  END SUBROUTINE read_header_r8

  SUBROUTINE read_constant_real_r8(h, value)

    TYPE(sdf_file_handle) :: h
    REAL(r8), INTENT(OUT) :: value
    REAL(r4) :: real4
    REAL(r8) :: real8
    TYPE(sdf_block_type), POINTER :: b

    CALL sdf_read_constant(h)

    b => h%current_block

    IF (b%datatype == c_datatype_real4) THEN
      real4 = TRANSFER(b%const_value, real4)
      value = REAL(real4,r8)
    ELSE IF (b%datatype == c_datatype_real8) THEN
      real8 = TRANSFER(b%const_value, real8)
      value = REAL(real8,r8)
    END IF

  END SUBROUTINE read_constant_real_r8

  SUBROUTINE read_1d_array_real_r8(h, values)

    TYPE(sdf_file_handle) :: h
    REAL(r8), DIMENSION(:), INTENT(OUT) :: values
    INTEGER :: errcode, n1
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    n1 = b%dims(1)

    h%current_location = b%data_location

    IF (h%rank == h%rank_master) THEN
      CALL MPI_FILE_READ_AT(h%filehandle, h%current_location, values, n1, &
          b%mpitype, MPI_STATUS_IGNORE, errcode)
    END IF

    CALL MPI_BCAST(values, n1, b%mpitype, h%rank_master, h%comm, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_1d_array_real_r8

  SUBROUTINE read_1d_array_par_real_spec_r8(h, values, subarray, distribution)

    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:), INTENT(OUT) :: values
    INTEGER, INTENT(IN) :: subarray, distribution
    INTEGER :: errcode
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    h%current_location = b%data_location
    CALL MPI_FILE_SET_VIEW(h%filehandle, h%current_location, MPI_BYTE, &
        distribution, 'native', MPI_INFO_NULL, errcode)
    CALL MPI_FILE_READ_ALL(h%filehandle, values, 1, subarray, &
        MPI_STATUS_IGNORE, errcode)
    CALL MPI_FILE_SET_VIEW(h%filehandle, c_off0, MPI_BYTE, MPI_BYTE, 'native', &
        MPI_INFO_NULL, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_1d_array_par_real_spec_r8

  SUBROUTINE read_1d_array_par_real_r8(h, values, sz, local_starts, &
      local_ghosts, null_proc)

    INTEGER, PARAMETER :: ndims = 1
    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:), INTENT(OUT) :: values
    INTEGER, DIMENSION(ndims), INTENT(IN) :: sz
    INTEGER, DIMENSION(ndims), INTENT(IN) :: local_starts
    INTEGER, DIMENSION(2*ndims), INTENT(IN), OPTIONAL :: local_ghosts
    LOGICAL, INTENT(IN), OPTIONAL :: null_proc
    INTEGER, DIMENSION(2*ndims) :: ghosts
    INTEGER :: distribution, subarray, errcode
    INTEGER, DIMENSION(ndims) :: starts, sizes, subsizes
    LOGICAL :: not_this_processor

    IF (sdf_check_block_header(h)) RETURN

    not_this_processor = .FALSE.
    IF (PRESENT(null_proc)) THEN
      not_this_processor = null_proc
    ELSE IF (ANY(local_starts < 0)) THEN
      not_this_processor = .TRUE.
    END IF

    IF (.NOT. not_this_processor) THEN
      ghosts = 0
      IF (PRESENT(local_ghosts)) ghosts = local_ghosts

      starts = local_starts - 1
      sizes = sz
      subsizes = SHAPE(values) - ghosts(1:ndims) - ghosts(ndims+1:)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)

      ! Subsizes are unchanged
      starts = ghosts(1:ndims)
      sizes = SHAPE(values)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, subarray, errcode)
      CALL MPI_TYPE_COMMIT(subarray, errcode)
    ELSE
      CALL MPI_TYPE_CONTIGUOUS(0, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)
      subarray = distribution
    END IF

    CALL read_1d_array_par_real_spec_r8(h, values, subarray, distribution)

    IF (subarray /= distribution) CALL MPI_TYPE_FREE(subarray, errcode)
    CALL MPI_TYPE_FREE(distribution, errcode)

  END SUBROUTINE read_1d_array_par_real_r8

  SUBROUTINE read_2d_array_real_r8(h, values)

    TYPE(sdf_file_handle) :: h
    REAL(r8), DIMENSION(:,:), INTENT(OUT) :: values
    INTEGER :: errcode, n1
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    n1 = b%dims(1) * b%dims(2)

    h%current_location = b%data_location

    IF (h%rank == h%rank_master) THEN
      CALL MPI_FILE_READ_AT(h%filehandle, h%current_location, values, n1, &
          b%mpitype, MPI_STATUS_IGNORE, errcode)
    END IF

    CALL MPI_BCAST(values, n1, b%mpitype, h%rank_master, h%comm, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_2d_array_real_r8

  SUBROUTINE read_2d_array_par_real_spec_r8(h, values, subarray, distribution)

    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:,:), INTENT(OUT) :: values
    INTEGER, INTENT(IN) :: subarray, distribution
    INTEGER :: errcode
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    h%current_location = b%data_location
    CALL MPI_FILE_SET_VIEW(h%filehandle, h%current_location, MPI_BYTE, &
        distribution, 'native', MPI_INFO_NULL, errcode)
    CALL MPI_FILE_READ_ALL(h%filehandle, values, 1, subarray, &
        MPI_STATUS_IGNORE, errcode)
    CALL MPI_FILE_SET_VIEW(h%filehandle, c_off0, MPI_BYTE, MPI_BYTE, 'native', &
        MPI_INFO_NULL, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_2d_array_par_real_spec_r8

  SUBROUTINE read_2d_array_par_real_r8(h, values, sz, local_starts, &
      local_ghosts, null_proc)

    INTEGER, PARAMETER :: ndims = 2
    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:,:), INTENT(OUT) :: values
    INTEGER, DIMENSION(ndims), INTENT(IN) :: sz
    INTEGER, DIMENSION(ndims), INTENT(IN) :: local_starts
    INTEGER, DIMENSION(2*ndims), INTENT(IN), OPTIONAL :: local_ghosts
    LOGICAL, INTENT(IN), OPTIONAL :: null_proc
    INTEGER, DIMENSION(2*ndims) :: ghosts
    INTEGER :: distribution, subarray, errcode
    INTEGER, DIMENSION(ndims) :: starts, sizes, subsizes
    LOGICAL :: not_this_processor

    IF (sdf_check_block_header(h)) RETURN

    not_this_processor = .FALSE.
    IF (PRESENT(null_proc)) THEN
      not_this_processor = null_proc
    ELSE IF (ANY(local_starts < 0)) THEN
      not_this_processor = .TRUE.
    END IF

    IF (.NOT. not_this_processor) THEN
      ghosts = 0
      IF (PRESENT(local_ghosts)) ghosts = local_ghosts

      starts = local_starts - 1
      sizes = sz
      subsizes = SHAPE(values) - ghosts(1:ndims) - ghosts(ndims+1:)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)

      ! Subsizes are unchanged
      starts = ghosts(1:ndims)
      sizes = SHAPE(values)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, subarray, errcode)
      CALL MPI_TYPE_COMMIT(subarray, errcode)
    ELSE
      CALL MPI_TYPE_CONTIGUOUS(0, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)
      subarray = distribution
    END IF

    CALL read_2d_array_par_real_spec_r8(h, values, subarray, distribution)

    IF (subarray /= distribution) CALL MPI_TYPE_FREE(subarray, errcode)
    CALL MPI_TYPE_FREE(distribution, errcode)

  END SUBROUTINE read_2d_array_par_real_r8

  SUBROUTINE read_3d_array_real_r8(h, values)

    TYPE(sdf_file_handle) :: h
    REAL(r8), DIMENSION(:,:,:), INTENT(OUT) :: values
    INTEGER :: errcode, n1
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    n1 = b%dims(1) * b%dims(2) * b%dims(3)

    h%current_location = b%data_location

    IF (h%rank == h%rank_master) THEN
      CALL MPI_FILE_READ_AT(h%filehandle, h%current_location, values, n1, &
          b%mpitype, MPI_STATUS_IGNORE, errcode)
    END IF

    CALL MPI_BCAST(values, n1, b%mpitype, h%rank_master, h%comm, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_3d_array_real_r8

  SUBROUTINE read_3d_array_par_real_spec_r8(h, values, subarray, distribution)

    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:,:,:), INTENT(OUT) :: values
    INTEGER, INTENT(IN) :: subarray, distribution
    INTEGER :: errcode
    TYPE(sdf_block_type), POINTER :: b

    IF (sdf_check_block_header(h)) RETURN

    b => h%current_block
    IF (.NOT. b%done_info) CALL sdf_read_array_info(h)

    h%current_location = b%data_location
    CALL MPI_FILE_SET_VIEW(h%filehandle, h%current_location, MPI_BYTE, &
        distribution, 'native', MPI_INFO_NULL, errcode)
    CALL MPI_FILE_READ_ALL(h%filehandle, values, 1, subarray, &
        MPI_STATUS_IGNORE, errcode)
    CALL MPI_FILE_SET_VIEW(h%filehandle, c_off0, MPI_BYTE, MPI_BYTE, 'native', &
        MPI_INFO_NULL, errcode)

    h%current_location = b%next_block_location
    b%done_data = .TRUE.

  END SUBROUTINE read_3d_array_par_real_spec_r8

  SUBROUTINE read_3d_array_par_real_r8(h, values, sz, local_starts, &
      local_ghosts, null_proc)

    INTEGER, PARAMETER :: ndims = 3
    TYPE(sdf_file_handle), INTENT(INOUT) :: h
    REAL(r8), DIMENSION(:,:,:), INTENT(OUT) :: values
    INTEGER, DIMENSION(ndims), INTENT(IN) :: sz
    INTEGER, DIMENSION(ndims), INTENT(IN) :: local_starts
    INTEGER, DIMENSION(2*ndims), INTENT(IN), OPTIONAL :: local_ghosts
    LOGICAL, INTENT(IN), OPTIONAL :: null_proc
    INTEGER, DIMENSION(2*ndims) :: ghosts
    INTEGER :: distribution, subarray, errcode
    INTEGER, DIMENSION(ndims) :: starts, sizes, subsizes
    LOGICAL :: not_this_processor

    IF (sdf_check_block_header(h)) RETURN

    not_this_processor = .FALSE.
    IF (PRESENT(null_proc)) THEN
      not_this_processor = null_proc
    ELSE IF (ANY(local_starts < 0)) THEN
      not_this_processor = .TRUE.
    END IF

    IF (.NOT. not_this_processor) THEN
      ghosts = 0
      IF (PRESENT(local_ghosts)) ghosts = local_ghosts

      starts = local_starts - 1
      sizes = sz
      subsizes = SHAPE(values) - ghosts(1:ndims) - ghosts(ndims+1:)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)

      ! Subsizes are unchanged
      starts = ghosts(1:ndims)
      sizes = SHAPE(values)

      CALL MPI_TYPE_CREATE_SUBARRAY(ndims, sizes, subsizes, starts, &
          MPI_ORDER_FORTRAN, MPI_REAL8, subarray, errcode)
      CALL MPI_TYPE_COMMIT(subarray, errcode)
    ELSE
      CALL MPI_TYPE_CONTIGUOUS(0, MPI_REAL8, distribution, errcode)
      CALL MPI_TYPE_COMMIT(distribution, errcode)
      subarray = distribution
    END IF

    CALL read_3d_array_par_real_spec_r8(h, values, subarray, distribution)

    IF (subarray /= distribution) CALL MPI_TYPE_FREE(subarray, errcode)
    CALL MPI_TYPE_FREE(distribution, errcode)

  END SUBROUTINE read_3d_array_par_real_r8

END MODULE sdf_input_r8


text:
C Nga Nguyen: This example shows a bug in regions
C J is modified by the loop but the region does not take into account
C because J is not analysed by semantics analyses.
C
C
      PROGRAM BUGREGION
      DIMENSION ITAB(2000)
      EQUIVALENCE (ITAB(10),J)
      J = 0
      DO I = 1,12
         J = J+1
         ITAB(J) = 1000
      ENDDO
      END

ext: f


!> This is the MPI version of the main program that interprets the
!! input file and decides which methods to execute with mpi implemented.
PROGRAM mqds_mpi
  USE input_output
  USE mpi_variables
  IMPLICIT NONE

  ! Initialize mpi & communication
  CALL MPI_INIT(ierr)
  CALL MPI_Comm_size(MPI_COMM_WORLD, npes, ierr)
  CALL MPI_Comm_rank(MPI_COMM_WORLD, mype, ierr)

  CALL read_input

  IF ( MOD( nbstep, dump ) /= 0 ) THEN
      OPEN(UNIT=10, FILE=ERRORLOG)
      WRITE(10,*) 'Error: nbstep must be divisible by dump!'
      CLOSE(10)
      STOP
  END IF

  ! mapping variable methods
  IF ( method == 'pldm' ) THEN
      CALL pldm_master_mpi
  END IF

  IF ( method == 'twa' ) THEN
      CALL twa_master_mpi
  END IF

  IF ( method == 'sqc' ) THEN
      CALL sqc_master_mpi
  END IF

  IF ( method == 'ipldm' ) THEN
      CALL ipldm_master_mpi
  END IF

  IF ( method == 'equilibrium' ) THEN
      CALL equilibrium_master_mpi
  END IF

  IF ( method == 'redfield' ) THEN
      CALL redfield_master_mpi
  END IF

  CALL MPI_FINALIZE(ierr)

END PROGRAM mqds_mpi


text:
program moist_adiabat
  use scale
  use scale_const, only: &
     PRE00 => CONST_PRE00, &
     Rdry  => CONST_Rdry,  &
     Rvap  => CONST_Rvap,  &
     CPdry => CONST_CPdry
  use scale_atmos_hydrometeor, only: &
     CPvap => CP_VAPOR, &
     CL    => CP_WATER, &
     CI    => CP_ICE
  use scale_io, only: &
     IO_get_available_fid
  use scale_atmos_saturation, only: &
     ATMOS_SATURATION_setup, &
     ATMOS_SATURATION_tdew_liq, &
     ATMOS_SATURATION_pote
  use scale_atmos_adiabat, only: &
     ATMOS_ADIABAT_setup, &
     ATMOS_ADIABAT_cape
  implicit none

  integer :: nx, ny, nz
  real(RP), allocatable :: xaxis(:), yaxis(:)
  real(RP), allocatable :: z(:), zh(:)

  real(RP), allocatable :: TEMP(:), PRES(:)
  real(RP), allocatable :: QV(:), QC(:), QR(:), QI(:), QS(:), QG(:)
  real(RP), allocatable :: POTT(:), DENS(:), QDRY(:)
  real(RP), allocatable :: Rtot(:), CPtot(:)
  real(RP), allocatable :: Tdew(:), POTE(:)
  real(RP), allocatable :: DENS_p(:), TEMP_p(:), BUOY_p(:), QV_p(:)
  real(RP) :: CAPE, CIN, LCL, LFC, LNB

  logical :: converged

  ! conf data
  character(len=H_MID) :: basename_in
  integer :: rankid
  character(len=H_MID) :: filename_out
  real(RP) :: x, y
  integer  :: nstep

  namelist / PARAM_MOIST_ADIABAT / &
       basename_in, &
       rankid, &
       filename_out, &
       x, y, &
       nstep

  integer :: fid
  integer :: k
  integer :: ierr

  call SCALE_init

  call ATMOS_SATURATION_setup
  call ATMOS_ADIABAT_setup

  ! default value
  basename_in  = "history"
  filename_out = "output.dat"
  x = 0.0_RP
  y = 0.0_RP
  nstep = 1

  !--- read namelist
  if ( IO_FID_CONF > 0 ) then
     rewind(IO_FID_CONF)
     read(IO_FID_CONF,nml=PARAM_MOIST_ADIABAT,iostat=ierr)
     if( ierr < 0 ) then !--- missing
        if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
     elseif( ierr > 0 ) then !--- fatal error
        write(*,*) 'xxx Not appropriate names in namelist PARAM_MOIST_ADIABAT. Check!'
        call PRC_abort
     endif
  end if
  if( IO_NML ) write(IO_FID_NML,nml=PARAM_MOIST_ADIABAT)

  call read_data

  ! approximate estimation (ignore other tracers)
  allocate( POTT(nz), DENS(nz), QDRY(nz), Rtot(nz), CPtot(nz) )
  do k = 1, nz
     QDRY(k) = 1.0_RP - QV(k) - QC(k) - QR(k) - QI(k) - QS(k) - QG(k)
     Rtot(k)  =  Rdry * QDRY(k) +  Rvap * QV(k)
     CPtot(k) = CPdry * QDRY(k) + CPvap * QV(k) + CL * ( QC(k) + QR(k) ) + CI * ( QI(k) + QS(k) + QG(k) )
     POTT(k) = TEMP(k) * ( PRE00 / PRES(k) )**( Rtot(k) / CPtot(k) )
     DENS(k) = PRES(k) / ( Rtot(k) * TEMP(k) )
  end do

  ! dew point
  allocate( Tdew(nz) )
  call ATMOS_SATURATION_tdew_liq( nz, 1, nz, &
                                  DENS(:), TEMP(:), QV(:), & ! [IN]
                                  Tdew(:), converged       ) ! [OUT]

  ! equivalent potential temperature
  allocate( POTE(nz) )
  call ATMOS_SATURATION_pote( nz, 1, nz, &
                              DENS(:), POTT(:), TEMP(:), QV(:), & ! [IN]
                              POTE(:)                           ) ! [OUT]

  ! lift parcel
  allocate( DENS_p(nz), TEMP_p(nz), BUOY_p(nz), QV_p(nz) )
  call ATMOS_ADIABAT_cape( nz, 1, nz, 1, &
                           DENS(:), TEMP(:), PRES(:),                & ! [IN]
                           QV(:), QC(:), QDRY(:), Rtot(:), CPtot(:), & ! [IN]
                           z(:), zh(:),                              & ! [IN]
                           CAPE, CIN, LCL, LFC, LNB,                 & ! [OUT]
                           DENS_p(:), TEMP_p(:), BUOY_p(:), QV_p(:), & ! [OUT]
                           converged                                 ) ! [OUT]
  if ( .not. converged ) then
     write(*,*) 'xxx calculation was not converged!'
     call PRC_abort
  end if
  !write(*,*) buoy_p
  !write(*,*) zh

  write(fid,'(5(a,f10.2))') 'CAPE=', CAPE, ', CIN=', CIN, ', LCL=', LCL, ', LFC=', LFC, ', LNB=', LNB

  fid = IO_get_available_fid()
  open( unit=fid, file=filename_out )
  write(fid,*)nz
  write(fid,'(2a8,5a10)') "z", "pres", "temp", "tdew", "pott", "pote", "temp_p"
  do k = 1, nz
     write(fid,'(2f8.1,5f10.2)') z(k), pres(k)/100.0_RP, temp(k), tdew(k), pott(k), pote(k), temp_p(k)
  end do
  close(fid)

  !write(*,*)temp_p(nz)*(PRE00/pres(nz))**(Rdry/CPdry)

  call SCALE_finalize

  stop

contains

  subroutine read_data
    use scale_file, only: &
         FILE_open, &
         FILE_get_shape, &
         FILE_read, &
         FILE_close
    integer :: fid
    integer :: dims(3)
    integer :: start(3), count(3)

    integer :: i, j

    ! file open
    call FILE_open( basename_in,  & ! (in)
                    fid,          & ! (out)
                    rankid=rankid ) ! (in)

    ! get dimension size
    call FILE_get_shape( fid, "height", & ! (in)
                         dims(:)        ) ! (out)
    nx = dims(1)
    ny = dims(2)
    nz = dims(3)

    ! get x and y axis
    allocate( xaxis(nx), yaxis(ny) )
    call FILE_read( fid, "x", xaxis(:) )
    call FILE_read( fid, "y", yaxis(:) )

    ! search index
    do i = 1, nx
       if ( x <= xaxis(i) ) exit
    end do
    do j = 1, ny
       if ( y <= yaxis(j) ) exit
    end do
    start(:) = (/i,j,1/)

    ! get z axis
    allocate( z(nz), zh(0:nz) )
    call FILE_read( fid, "height",     z(:),  start=start(:), count=(/1,1,nz/) )
    call FILE_read( fid, "height_xyw", zh(:), start=start(:), count=(/1,1,nz+1/) )

    ! get variables
    allocate( TEMP(nz), PRES(nz), QV(nz), QC(nz), QR(nz), QI(nz), QS(nz), QG(nz) )
    count(:) = (/1,1,nz/)

    call FILE_read( fid, "T",    TEMP(:), start=start(:), count=count(:), step=nstep )
    call FILE_read( fid, "PRES", PRES(:), start=start(:), count=count(:), step=nstep )
    call FILE_read( fid, "QV",   QV  (:), start=start(:), count=count(:), step=nstep )
    call FILE_read( fid, "QC",   QC  (:), start=start(:), count=count(:), step=nstep, allow_missing=.true. )
    call FILE_read( fid, "QR",   QR  (:), start=start(:), count=count(:), step=nstep, allow_missing=.true. )
    call FILE_read( fid, "QI",   QI  (:), start=start(:), count=count(:), step=nstep, allow_missing=.true. )
    call FILE_read( fid, "QS",   QS  (:), start=start(:), count=count(:), step=nstep, allow_missing=.true. )
    call FILE_read( fid, "QG",   QG  (:), start=start(:), count=count(:), step=nstep, allow_missing=.true. )

    ! close
    call FILE_close( fid )

    return
  end subroutine read_data

end program moist_adiabat

max_stars_repo_name: slayoo/scale

max_issues_repo_name: slayoo/scale

max_forks_repo_name: slayoo/scale


PROGRAM CARBON

    ! ************************************************************************
    !
    !
    !         CARBON ACCUMULATION PER LATITUDE BAND
    !
    !
    ! ************************************************************************

    ! Version 1.0: written by Aitor Aldama-Campino, Autumn 2018

    USE netcdf

    IMPLICIT NONE

    INTEGER, PARAMETER     :: IX = 1440, JYD = 1, JYU = 1021, JY = 1021, KZ = 75
    INTEGER, DIMENSION(4)  :: start2D, count2D, start3D, count3D
    INTEGER                :: ncid, varid, ierr, ii, jj, jp, kk, nindex, mm

    REAL, DIMENSION(IX,JYD:JY)    :: dx, dy, kmt
    REAL, DIMENSION(IX,JYD:JY,KZ) :: dz, ctoti,ctotf, tmask
    REAL, DIMENSION(120,JYD:JYU)  :: accu
    REAL, PARAMETER           :: dt = 1., tmin = -350., trend = 9.11518581e-03

    LOGICAL :: file_exists

    CHARACTER (LEN=400)                :: filepath, commandline
    CHARACTER (LEN=20), DIMENSION(120) :: label1, label2, label3

    ! ***************************************************************

    ! Path to the data input directory
    filepath = '/User/Input/'

    ! READING REGIONS AND MESH
    ! ===============================================================
    start2D = [ 1,JYD,1,1]
    count2D = [IX,JY-JYD+1,1,1]

    start3D = [ 1, JYD,1,1]
    count3D = [IX, JY-JYD+1,KZ,1]

    ! Dx and Dz / Bathimetry
    ierr=NF90_OPEN(trim(filepath)//'mesh_mask.nc',NF90_NOWRITE,ncid)
    IF(ierr.ne.0) STOP 'No mesh mask file'
    ierr=NF90_INQ_VARID(ncid,'e1t',varid) ! the main data fields
    IF(ierr.ne.0) STOP 7
    ierr=NF90_GET_VAR(ncid,varid,dx,start2d,count2d)
    IF(ierr.ne.0) STOP 8
    ierr=NF90_INQ_VARID(ncid,'e2t',varid) ! the main data fields
    IF(ierr.ne.0) STOP 7
    ierr=NF90_GET_VAR(ncid,varid,dy,start2d,count2d)
    IF(ierr.ne.0) STOP 8
    ierr=NF90_INQ_VARID(ncid,'e3t',varid) ! the main data fields
    IF(ierr.ne.0) STOP 9
    ierr=NF90_GET_VAR(ncid,varid,dz,start3d,count3d)
    IF(ierr.ne.0) STOP 10
    ierr=NF90_INQ_VARID(ncid,'tmask',varid) ! the main data fields
    IF(ierr.ne.0) STOP 11
    ierr=NF90_GET_VAR(ncid,varid,tmask,start3d,count3d)
    IF(ierr.ne.0) STOP 12
    ierr=NF90_INQ_VARID(ncid,'bathy',varid) ! the main data fields
    IF(ierr.ne.0) STOP 13
    ierr=NF90_GET_VAR(ncid,varid,kmt,start2d,count2d)
    ierr=NF90_CLOSE(ncid)

    PRINT*, '- Reading mesh and regions done'

    ! ***************************************************************

    ! DATE
    ! ===============================================================
    kk = 1
    DO ii = 2000,2009
        DO jj = 1,12
            WRITE(label1(kk),'(i4,A5,i4,A4)') ii,'0101_',ii,'1230'
            IF (jj<10) THEN
                WRITE(label2(kk),'(i4,A1,i1)') ii,'0',jj
            ELSE
                WRITE(label2(kk),'(i4,i02)') ii,jj
            ENDIF
            WRITE(label3(kk),'(i4)') ii

            kk = kk + 1
        END DO
    END DO

    ! ***************************************************************

    ! BIG LOOP
    ! ===============================================================

    start3D = [ 1, JYD,1,1]
    count3D = [IX, JY-JYD+1,KZ,1]

    accu = 0.
    DO nindex = 1,SIZE(label1)

        IF (nindex == 1) THEN
            PRINT*, '- Time step ::', nindex

            ! Tracer file
            ! Read the DIC
            ierr=NF90_OPEN(trim(filepath)//"/ORCA025-ROAM_1m_" &
                //TRIM(label1(nindex))//"_ptrc_T_"//TRIM(label2(nindex))//'-'//TRIM(label2(nindex))//".nc",NF90_NOWRITE,ncid)
            IF(ierr.ne.0) STOP 'No C1 file'
            ierr=NF90_INQ_VARID(ncid,'DIC',varid) ! the main data fields
            IF(ierr.ne.0) STOP 16
            ierr=NF90_GET_VAR(ncid,varid,ctoti,start3d,count3d)
            IF(ierr.ne.0) STOP 17
            ierr=NF90_CLOSE(ncid)
            WHERE (ctoti > 1e6) ctoti = 0.

            CALL SYSTEM("rm /Users/aitoraldama/Desktop/dumBC1/C1.nc")
        END IF

        PRINT*, '- Time step ::', nindex

        ! Tracer file
        ! Read the DIC
        ierr=NF90_OPEN(trim(filepath)//"/ORCA025-ROAM_1m_" &
            //TRIM(label1(nindex))//"_ptrc_T_"//TRIM(label2(nindex))//'-'//TRIM(label2(nindex))//".nc",NF90_NOWRITE,ncid)
        IF(ierr.ne.0) STOP 'No C1 file'
        ierr=NF90_INQ_VARID(ncid,'DIC',varid) ! the main data fields
        IF(ierr.ne.0) STOP 16
        ierr=NF90_GET_VAR(ncid,varid,ctotf,start3d,count3d)
        IF(ierr.ne.0) STOP 17
        ierr=NF90_CLOSE(ncid)
        WHERE (ctotf > 1e6) ctotf = 0.

        DO jj = 1, JY
            DO ii = 1, IX
                DO kk = 1, KZ
                    accu(nindex,jj) = accu(nindex,jj) + dz(ii,jj,kk)*dx(ii,jj)*dy(ii,jj)*&
                        &(ctotf(ii,jj,kk)-ctoti(ii,jj,kk))*tmask(ii,jj,kk)
                END DO
            END DO
        END DO

    END DO

    ! Save file
    OPEN(UNIT=11,FILE=TRIM(filepath)//"result_accu.bin", STATUS='REPLACE', FORM='UNFORMATTED')
    WRITE(11) accu
    CLOSE(11)

END PROGRAM CARBON


! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program netcdf_to_gitm_blocks

!----------------------------------------------------------------------
! purpose: interface between DART and the GITM model
!
! method: Read DART state netcdf files and overwrite values in a gitm restart file.
!
! NOT FINISHED!  HARDLY STARTED.
! this version assumes that the grid is global and the data needs to be
! blocked into one block per gitm mpi task.  there is a different converter
! for when gitm only needs a single input/output file.
!
!----------------------------------------------------------------------

use        types_mod, only : r8

use    utilities_mod, only : initialize_utilities, finalize_utilities,   &
                             find_namelist_in_file, check_namelist_read, &
                             open_file, close_file, E_MSG, error_handler

use        model_mod, only : netcdf_to_restart_files

use time_manager_mod, only : time_type, print_time, print_date, operator(-), &
                             get_time, get_date

implicit none

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

character(len=*),   parameter :: progname = 'netcdf_to_gitm_blocks'

!-----------------------------------------------------------------------
! namelist parameters with default values.
!-----------------------------------------------------------------------

character (len = 256) :: gitm_restart_input_dirname  = 'none'
character (len = 256) :: gitm_restart_output_dirname = 'none'
character (len = 256) :: netcdf_to_gitm_blocks_input_file = 'filter_restart.nc'

namelist /netcdf_to_gitm_blocks_nml/   &
     gitm_restart_input_dirname,       &
     gitm_restart_output_dirname,      &
     netcdf_to_gitm_blocks_input_file

!----------------------------------------------------------------------
! global storage
!----------------------------------------------------------------------

integer               :: iunit, io
character(len=512)    :: string1, string2, string3

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

call initialize_utilities(progname=progname)

!----------------------------------------------------------------------
! Read the namelist.
!----------------------------------------------------------------------

call find_namelist_in_file("input.nml", "netcdf_to_gitm_blocks_nml", iunit)
read(iunit, nml = netcdf_to_gitm_blocks_nml, iostat = io)
call check_namelist_read(iunit, io, "netcdf_to_gitm_blocks_nml")

call error_handler(E_MSG,progname,'','',revision,revdate)
write(string1,*) 'converting DART file ', "'"//trim(netcdf_to_gitm_blocks_input_file)//"'"
write(string2,*) 'to gitm restart files in directory ', "'"//trim(gitm_restart_output_dirname)//"'"
write(string3,*) 'using the restart files in directory ', "'"//trim(gitm_restart_input_dirname)//"' as a template"
call error_handler(E_MSG,progname,string1,source,revision,revdate,text2=string2,text3=string3)

!----------------------------------------------------------------------
! Reads the valid time, the state, and the target time.
!----------------------------------------------------------------------

call netcdf_to_restart_files(netcdf_to_gitm_blocks_input_file,gitm_restart_output_dirname,&
                             gitm_restart_input_dirname)

!----------------------------------------------------------------------
! Log what we think we're doing, and exit.
!----------------------------------------------------------------------
call error_handler(E_MSG,progname,'','',revision,revdate)
call error_handler(E_MSG,progname,'','',revision,revdate)
write(string1,*) 'Successfully converted to the gitm restart files in directory'
write(string2,*) "'"//trim(gitm_restart_output_dirname)//"'"
call error_handler(E_MSG,progname,string1,source,revision,revdate,text2=string2)

! end - close the log, etc
call finalize_utilities()

!======================================================================
contains
!======================================================================

subroutine write_gitm_time_control(model_time, adv_to_time)
! The idea is to write a text file with the following structure:
!
!#TIMESTART
!2003            year
!06              month
!21              day
!00              hour
!00              minute
!00              second
!
!#TIMEEND
!2003            year
!07              month
!21              day
!00              hour
!00              minute
!00              second
!

type(time_type), intent(in) :: model_time, adv_to_time
integer :: iyear,imonth,iday,ihour,imin,isec

iunit = open_file('DART_GITM_time_control.txt', action='write')
write(iunit,*)

! the end time comes first.

call get_date(adv_to_time,iyear,imonth,iday,ihour,imin,isec)
write(iunit,'(''#TIMEEND'')')
write(iunit,'(i4.4,10x,''year''  )')iyear
write(iunit,'(i2.2,12x,''month'' )')imonth
write(iunit,'(i2.2,12x,''day''   )')iday
write(iunit,'(i2.2,12x,''hour''  )')ihour
write(iunit,'(i2.2,12x,''minute'')')imin
write(iunit,'(i2.2,12x,''second'')')isec
write(iunit,*)

call get_date(model_time,iyear,imonth,iday,ihour,imin,isec)
write(iunit,'(''#TIMESTART'')')
write(iunit,'(i4.4,10x,''year''  )')iyear
write(iunit,'(i2.2,12x,''month'' )')imonth
write(iunit,'(i2.2,12x,''day''   )')iday
write(iunit,'(i2.2,12x,''hour''  )')ihour
write(iunit,'(i2.2,12x,''minute'')')imin
write(iunit,'(i2.2,12x,''second'')')isec
write(iunit,*)

call close_file(iunit)
end subroutine write_gitm_time_control

end program netcdf_to_gitm_blocks

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$


! Run with -np 5 and enter s_size=10.
program scatter
    use mpi_f08
    implicit none
    integer :: my_rank, n_ranks, root, s_size, r_size, i
    type(MPI_Comm) :: comm
    integer, allocatable :: rbuf(:), sbuf(:)

    call MPI_Init()
    comm = MPI_COMM_WORLD
    call MPI_Comm_rank(comm, my_rank)
    call MPI_Comm_size(comm, n_ranks)

    root = 0

    ! Allocate the send buffer `sbuf` of size `s_size` only at rank 0 to scatter
    ! it to the other ranks.
    s_size = 0
    if (my_rank == root) read *, s_size
    allocate(sbuf(s_size))
    if (my_rank == root) sbuf(:) = [ (i, i = 1, s_size) ]

    ! Broadcast the array size and create the receive buffer `sbuf` of size
    ! r_size = s / n.  Initialize it to 0 at all ranks.
    call MPI_Bcast(s_size, 1, MPI_INTEGER, root, comm)
    r_size = s_size / n_ranks
    allocate(rbuf(r_size), source=0)

    call print_barrier()

    call MPI_Scatter(sbuf, r_size, MPI_INTEGER, &
                     rbuf, r_size, MPI_INTEGER, root, comm)

    call print_barrier()

    if (allocated(sbuf)) deallocate(sbuf)
    if (allocated(rbuf)) deallocate(rbuf)

    call MPI_Finalize()

contains

    !---------------------------------------------------------------------------
    subroutine print_barrier()
        integer :: rank

        if (my_rank == 0) print *
        do rank = 0, n_ranks - 1
            if (rank == my_rank) then
                print "(a, i0, a, *(i2, 1x))", "Rank ", rank, ", sbuf(:)=", sbuf
            end if
            call MPI_Barrier(comm)
        end do
        do rank = 0, n_ranks - 1
            if (rank == my_rank) then
                print "(a, i0, a, *(i2, 1x))", "Rank ", rank, ", rbuf(:)=", rbuf
            end if
            call MPI_Barrier(comm)
        end do
    end subroutine print_barrier
    !---------------------------------------------------------------------------

end program scatter


text:
C*********************************************************************
      subroutine STEPG(KK, NN, TSTEP, NT, T1, TN, T, NSTEP)
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  *
C
C   THE ROUTINE GENERATES STEPS IN LINEAR OR LOGARITHMIC INCREMENT
C
C        KK      = SWITCH    KK = 1    STEP INCREASE OF VALUES
C                            KK = 2    LOGARITHMIC INCREASE OF VALUES
C        NN      = NUMBER OF STEPS   OR NUMBER OF VALUES IN EACH 10
C        TSTEP   = LARGEST VALUE IN EACH STEP
C        NT      = NUMBER OF VALUES IN EACH STEP
C        T1      = FIRST VALUE IN LOG-STEP
C        TN      = LAST  VALUE IN LOG-STEP
C        T       = VALUES GENERATED
C        NSTEP   = NUMBER OF VALUES
C
C   CODED PER B SCHNABEL SEPT. 1970
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  *
C
      dimension T(200), TSTEP(27), NT(27)
C
      if (KK .eq. 1) then
        K = 1
        T(K) = 0.
        SAVE = 0.
        do N = 1,NN
          M = NT(N)
          STEP = (TSTEP(N) - SAVE)/FLOAT(M)
          SAVE = TSTEP(N)
          do I = 1,M
            K = K + 1
            T(K) = T(K-1) + STEP
          end do
        end do
        NSTEP = K
        return
      else if (KK .eq. 2) then
        NST = ALOG10(T1)
        if (T1.lt. 1.) NST = NST - 1
        STEP = 1./NN
        K = 1
        TA = 10.**FLOAT(NST)
        T(1) = TA
        do J = 2,NN
          K = K + 1
          T(K) = TA*10.**(STEP*FLOAT(J))
          if (T(K) .gt. T1) go to 221
        end do
  221   TA = T(K-1)
        K = 0
  211   do J = 1,NN
          K = K + 1
          T(K) = TA*10.**(STEP*FLOAT(J))
          if (T(K) .gt. TN)  go to 212
        end do
        TA = TA*10.
        go to 211
  212   NSTEP = K
        return
      end if
      end