!---------------------------------------------------------------------
! README
!
! The attached program will generate an input file for the MAS model.
! This file called 
!      mas_inputfile.dat
! will contain photospheric field values (assumed to be radial)
! on a latitude/longitude grid of size 90x181 defined by the
! user.
! The user is required to modify the program below so as to 
! define the theta grid points in the array x_rdtxt2d, 
! the phi grid points in the array y_rdtxt2d, and the field
! values in the array f_rdtxt2d. Comments direct the user to
! the appropriate places for these modifications.
!
! Note, any structure included in this artificial magnetogram file
! will need to be on a scale of at least the MAS model resolution, if
! it is to appear in the subsequent MAS model output.
!
!---------------------------------------------------------------------



      module rdtxt2d_arrays

! ****** Arrays for writing data in routine WRTXT2D.

      real,allocatable :: f_rdtxt2d(:,:),
     &                    x_rdtxt2d(:),
     &                    y_rdtxt2d(:)
      integer :: ndim,nx,ny,ifscale

      end module


      program mas_magnetogram
      use rdtxt2d_arrays

      implicit none
      character(len=21) :: fname
      integer           :: i,j
      real,parameter    :: pi = 3.14159265358979

      ndim = 2
      ifscale = 1

!----------
! User set
      nx = 90                     ! user modify size of magnetogram grid in
                                  ! latitude(theta) direction
      ny = 181                    ! user modify size of magnetogram grid in
                                  ! azimuthal(phi) direction
      allocate (f_rdtxt2d(nx,ny))
      allocate (x_rdtxt2d(nx))
      allocate (y_rdtxt2d(ny))

      do i=1,nx
!       x_rdtxt2d(i) = pi/real(nx-1) * real(i-1)        ! example
        x_rdtxt2d(i) = ?????                            ! user enter theta grid points
      enddo
      do i=1,ny
!       y_rdtxt2d(i) = 2.*pi/real(ny-1) * real(i-1)     ! example
        y_rdtxt2d(i) = ?????                            ! user enter phi grid points
      enddo
      do j=1,ny
      do i=1,nx
        f_rdtxt2d(i,j) = ?????    ! user enter Br at every grid point on
                                  ! solar surface
      enddo
      enddo
!----------

      fname='mas_inputfile.dat'

      call wrtxt2d (fname)

      deallocate (f_rdtxt2d)
      deallocate (x_rdtxt2d)
      deallocate (y_rdtxt2d)

      stop
      end program mas_magnetogram



      subroutine wrtxt2d (fname)
      use rdtxt2d_arrays
      character*(*) :: fname

      open (20,file=fname,ACTION='WRITE')
c
c ****** Write the number of dimensions.
c
      write (20,*) ndim
c
c
c ****** Write the number of theta and phi coordinates.
c
      write (20,*) nx
      write (20,*) ny
c
c ****** Mark that the scales are present.
c
      write (20,*) ifscale
c
c ****** Write the theta and phi coordinate vectors
c
      write (20,*) x_rdtxt2d
      write (20,*) y_rdtxt2d
c
c ****** Write the array of Br values on the solar surface.
c
      write (20,*) f_rdtxt2d



      close (20)

      return
      end subroutine wrtxt2d