fortran source that causes ICE.

readmajorcube.f (text/plain), 6.78 KB, created by Tom Fogal on 2005-08-19 11:26:21 UTC

c **********************************************************************
      subroutine ReadMajorCube(iounit,Imode,
     %           HotCubeIndex,MajorCubeNumber)
c
c Purpose :: Read in the values of the interpolation cubes from files
c
c Programmer: Douglas J. Larson
c Creation Date : Friday, May 30, 1997
c Affiliation : University of New Hampshire -  Space Plasma Physics
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      common/cpu_info/cputype ! Generic across all applications
      integer cputype       ! Computer type (0-G:1-V:2-I:3-R4.4k:4-R10k)
      include 'INTERP.FOR'
c
c Local Logical Variables :
c -------------------------
      logical firstpass ! TRUE-Write header, FALSE-No header needed
c
c Local Character Variables :
c ---------------------------
      character*1 ans
      character*80 question
      character*80 header(7)
      character*128 filename
c
c Integer Arguments :
c -------------------
      integer iounit
      integer Imode  ! 1-Electric Field, 2-magnetic field, 3-potential
      integer MajorCubeNumber ! Index to the MajorCube parameter arrays
      integer HotCubeIndex    ! Where to put the MC values
c
c Local Integer Variables :
c -------------------------
      integer debug_level
      integer rwflag     ! 0-Read, 1-Write
      integer typecount(4)    ! 1-Bfield, 2-Birkland, 3-Xtail, 4-Efield
      integer HeaderLineCount
      integer NumberofColumns
      integer filecount  ! Number of Major Cubes to process
      integer i,j,k   ! Loop indices
      integer ldim       ! Actual number of data points
      integer FLtype     ! Type of field line system (B, Xtail, etc.)
      integer SLtype     ! 1-E, 2-B, 3-j, 4-Divj, 5-K, 6-JxB, etc.
      integer inputstyle ! 0-Interactive, 1-File, 2-Grid setup
      integer coorsys    ! 1-Cartesian, 2-Cylindrical, 3-Spherical
      integer range_total(3) ! Number of boxes in the range
      integer currentpoint
c
c Local Real Variables :
c ----------------------
      real*8 range(6)  ! Range start and end for each axes
      real*8 range_delta(3)     ! Range delta for each axes
      real*8 x,y,z              ! position
      real*8 array(maxvalues,7) ! (x,y,z,Vecx,Vecy,Vecz)
c
c **********************************************************************
c *****                                                            *****
c *****     R E A D  I N  M A J O R C U B E  D A T A F I L E S     *****
c *****                                                            *****
c **********************************************************************
c
      debug_level=0
      inputstyle=3
      HeaderLineCount=7
      NumberofColumns=7

      if(debug_level.gt.0)then
         print *,'ReadMajorCube: Imode=',Imode
         print *,'ReadMajorCube: filecount=',filecount
         print *,'ReadMajorCube: NumberofColumns=',NumberofColumns
         print *,'ReadMajorCube: HeaderLineCount=',HeaderLineCount
      end if

      rwflag=0
      if(Imode.eq.1)then
         filename=EI_filenames(MajorCubeNumber)

         if(debug_level.gt.0)then
            print *,'ReadMajorCube: filename=',filename
            print *,'ReadMajorCube: MajorCubeNumber=',MajorCubeNumber
            print *,'ReadMajorCube: EI_TimeDep=',EI_TimeDep
         end if

         call flsetup(iounit,filename)

         call rw_mcdata(iounit,0,"NA",filename,HeaderLineCount,
     %        NumberofColumns,header,maxvalues,ldim,array)

         currentpoint=1
         do i=1,EI_rangetotal(1,MajorCubeNumber)
            do j=1,EI_rangetotal(2,MajorCubeNumber)
               do k=1,EI_rangetotal(3,MajorCubeNumber)
                  EI_ex(i,j,k,HotCubeIndex)=array(currentpoint,4)
                  EI_ey(i,j,k,HotCubeIndex)=array(currentpoint,5)
                  EI_ez(i,j,k,HotCubeIndex)=array(currentpoint,6)
                  if(debug_level.gt.0)then
                     print 100,currentpoint,i,j,k,
     %                     EI_ex(i,j,k,HotCubeIndex),
     %                     EI_ey(i,j,k,HotCubeIndex),
     %                     EI_ez(i,j,k,HotCubeIndex) 
                  end if
100               format(4(i6),1x,3(1pd23.16,1x))
                  currentpoint=currentpoint+1
               end do
            end do
         end do
      else if(Imode.eq.2)then
         filename=BI_filenames(MajorCubeNumber)

         if(debug_level.gt.0)then
            print *,'ReadMajorCube: filename=',filename
            print *,'ReadMajorCube: MajorCubeNumber=',MajorCubeNumber
            print *,'ReadMajorCube: BI_TimeDep=',BI_TimeDep
         end if

         call flsetup(iounit,filename)

         call rw_mcdata(iounit,0,"NA",filename,HeaderLineCount,
     %        NumberofColumns,header,maxvalues,ldim,array)

         currentpoint=1
         do i=1,BI_rangetotal(1,MajorCubeNumber)
            do j=1,BI_rangetotal(2,MajorCubeNumber)
               do k=1,BI_rangetotal(3,MajorCubeNumber)
                  BI_bx(i,j,k,HotCubeIndex)=array(currentpoint,4)
                  BI_by(i,j,k,HotCubeIndex)=array(currentpoint,5)
                  BI_bz(i,j,k,HotCubeIndex)=array(currentpoint,6)
                  if(debug_level.gt.0)then
                     print 100,currentpoint,i,j,k,
     %                        BI_bx(i,j,k,HotCubeIndex),
     %                        BI_by(i,j,k,HotCubeIndex),
     %                        BI_bz(i,j,k,HotCubeIndex) 
                  end if
                  currentpoint=currentpoint+1
               end do
            end do
         end do
      else if(Imode.eq.3)then
         NumberofColumns=4
         filename=PI_filenames(MajorCubeNumber)

         if(debug_level.gt.0)then
            print *,'ReadMajorCube: filename=',filename
            print *,'ReadMajorCube: MajorCubeNumber=',MajorCubeNumber
            print *,'ReadMajorCube: PI_TimeDep=',PI_TimeDep
         end if

         call flsetup(iounit,filename)

         call rw_mcdata(iounit,0,"NA",filename,HeaderLineCount,
     %        NumberofColumns,header,maxvalues,ldim,array)

         currentpoint=1
         do i=1,PI_rangetotal(1,MajorCubeNumber)
            do j=1,PI_rangetotal(2,MajorCubeNumber)
               do k=1,PI_rangetotal(3,MajorCubeNumber)
                  PI_data(i,j,k,HotCubeIndex)=array(currentpoint,4)
                  if(debug_level.gt.0)then
                     print 200,currentpoint,i,j,k,
     %                     PI_data(i,j,k,HotCubeIndex)
                  end if
200               format(4(i4),1(f10.4))
                  currentpoint=currentpoint+1
               end do
            end do
         end do
cc         print *,'setup_Interp: Imode=P is not implemented yet'
cc         if(SLtype.ne.10)then
cc            print *,'ReadMajorCube: You requested an E-potential'
cc            print *,'*** file, this file is SLtype=',SLtype 
cc             print *,'*** Sorry, This is a FATAL error'
cc             stop
c          end if
      end if
      return
      end