!----------------------------------------------------------------------------
!
! This fortran subroutine outputs model solution data cube with the correct 
! format for the SHINE validation study
!
!  Written : Peter MacNeice                March 2012


!  The steps to using the DATACUBE model are:
!
!    1 - Make sure that the output from your simulation is in spherical
!        coordinates with cgs units. The grid must be regular and ordered,
!        but does not need to be uniform. At minimum, your inital simulation 
!        must include r, theta, phi, br, btheta, bphi. You can also include 
!        density, temperature, and pressure, but they are optional. If your 
!        cube is not axisymmetric, bear in mind that for_drive expects models 
!        to be centered on phi = 0, theta=90 (equator). However, the default 
!        POS plot will be from the viewer position of central meridian -90, 
!        so that your cube will be seen centered at the West (right) limb. 
!        Check FORWARD/MODELS/DATACUBE/make_my_cube.pro for more detailed 
!        information about what is needed from your initial simulation.
!
! The arguments to this routine are
 
!   nr                integer                                   number of radial grid points
!   nth               integer                                   number of latitudinal grid points
!   nph               integer                                   number of longitudinal grid points
!   r                 float vector    dimension nr              radial grid point coords 
!   th                float vector    dimension nth             latitudinal grid point coords 
!   ph                float vector    dimension nph             longitudinal grid point coords 
!   br                float vector    dimension nr*nth*nph      vector storing radial component of magnetic field
!   bth               float vector    dimension nr*nth*nph      vector storing latitudinal component of magnetic field
!   bph                float vector   dimension nr*nth*nph      vector storing longitudinal component of magnetic field
!   dens              float vector    dimension nr*nth*nph      vector storing density at grid points
!   temp              float vector    dimension nr*nth*nph      vector storing temperature at grid points
!   pres              float vector    dimension nr*nth*nph      vector storing pressure at grid points
!   units_descriptor  character string                          string describing units of variables
!   filename          character string                          string containing the output filename
                             
                             
!   The 3D data in br,bth,bph,dens,temp and pres is passed into this routine in a flar 1D vector with
!   the radial coordinate index changing fastest, followed by the latitudinal index and finally the longitudinal
!   index changing slowest.
                                   
                                   
!
!----------------------------------------------------------------------------


         subroutine shine_valid_format_fortran( nr, nth, nph,
     .                                          r, th, phi,
     .                                          br, bth, bph, 
     .                                          dens, temp, pres,
     .                                          units_descriptor,
     .                                          filename
     .                                        )

         implicit none

         integer, intent(in)  :: nr, nth, nph
         real,intent(in)      :: r(nr), th(nth), phi(nph)
         real,intent(in)      :: br(nr,nth,nph)
         real,intent(in)      :: bth(nr,nth,nph)
         real,intent(in)      :: bph(nr,nth,nph)
         real,intent(in)      :: dens(nr,nth,nph)
         real,intent(in)      :: temp(nr,nth,nph)
         real,intent(in)      :: pres(nr,nth,nph)
         character(len=80), intent(in) :: units_descriptor
         character(len=*), intent(in)  :: filename

         integer              :: i, j, k



! Open file
         open(unit=100,status='unknown',file=filename)

! Write units description string (eg. 'Units are cgs')
         write(100,'(a80)') units_descriptor

! Write grid dimensions
         write(100,'(i4,1x,i4,1x,i4)')  nr, nth, nph

! Write radial grid points
         do i=1,nr
           write(100,'(1pe13.6)') r(i)
         enddo

! Write latitudinal grid points
         do i=1,nth
           write(100,'(1pe13.6)') th(i)
         enddo

! Write azimuthal grid points
         do i=1,nph
           write(100,'(1pe13.6)') phi(i)
         enddo

! Write magnetic field values
         do k=1,nph
         do j=1,nth
         do i=1,nr
           write(100,'(1pe13.6,1x,1pe13.6,1x,1pe13.6)') 
     .                   br(i,j,k), bth(i,j,k), bph(i,j,k)
         enddo
         enddo
         enddo

! Write thermodynamic quantities
         do k=1,nph
         do j=1,nth
         do i=1,nr
           write(100,'(1pe13.6,1x,1pe13.6,1x,1pe13.6)') 
     .                   dens(i,j,k), temp(i,j,k), pres(i,j,k)
         enddo
         enddo
         enddo


         close(unit=100)


         return
         end subroutine shine_valid_format_fortran