subroutine readModel (A, B, lmax, kmax, choiceD)

c  This subroutine reads a standard model Harvard 3D and 2D model
c  coefficients (and converts 2D coefficients to real spherical 
c  harmonics).

c  April 26, 1999.
c  Last Modified:  September 07, 1999.

cccccccccc
c  Output:
c	A(l,m,k) = cosine coefficients
c	B(l,m,k) = sine coefficients
c	lmax = maximum angular degree
c	kmax = maximum order of the radial expansion
c	choiceD = 1 for 2D model and 2 for 3D model
cccccccccc

c=======Declaration===========

      integer lmax, kmax
      real A(0:lmax,0:lmax,0:kmax), B(0:lmax,0:lmax,0:kmax)

      character*100, inputfile
c--full filename from which the mode is read.

      integer choice_model
c  a value determining which model file is going to be read.

      integer choiceD
c  this is 1 for 2D model and 2 for 3D model.

      integer tmpl, tmpm 
c  temporary values of l and m to check.

      integer l, m, k

c======PROGRAM===============

      do l = 0, lmax
	 do m = 0, lmax
	    do k = 0, kmax
	       A(l,m,k) = 0.0
	       B(l,m,k) = 0.0  ! cleaning up
	    end do  ! end k loop
	 end do  ! end m loop
      end do  ! end l loop

 35   if (choiceD.eq.2) then  ! 3D model

         print*, 'Choose one of the following models.'
         print*, '	1. Shear Velocity Model'
         print*, '	2. Compressional Velocity Model'
         print*, '	3. Density Model'
         read*, choice_model
     
         if (choice_model.eq.1) then
	    inputfile = 'sprd6_s'
         elseif (choice_model.eq.2) then
	    inputfile = 'sprd6_p'
         elseif (choice_model.eq.3) then
	    inputfile = 'sprd6_r'
         else  ! if a choice not allowed
	    print*, choice_model, ' is not one of the options.'
	    go to 35
         end if  ! if (choice_model.eq.1)...

      elseif (choiceD.eq.1) then  ! 2D model

         print*, 'Choose one of the following models.'
         print*, '	1. Free-Surface Dynamic Topography Model'
         print*, '	2. 660-km Discontinuity Model'
         print*, '	3. Core-Mantle Boundary Model'
         read*, choice_model
     
         if (choice_model.eq.1) then
	    inputfile = 'sprd6_fs.XCS'
         elseif (choice_model.eq.2) then
	    inputfile = 'sprd6_d660.XCS'
         elseif (choice_model.eq.3) then
	    inputfile = 'sprd6_dcmb.XCS'
         else  ! if a choice not allowed
	    print*, choice_model, ' is not one of the options.'
	    go to 35
         end if  ! if (choice_model.eq.1)...

      end if  ! if (choiceD.eq.2)...

      open(10, file = inputfile) ! opening the model file
      rewind(10)

      if (choiceD.eq.2) then  ! if 3D model

         do k = 0, kmax
	    do l = 0, lmax
	    read(10,*), A(l,0,k), (A(l,m,k), B(l,m,k), m = 1, l)

	       if (B(l,0,k).ne.0.0) then  ! check
	          print*, '!!Error:  B(l,0,k).ne.0.0 for k = ', k, ', l = ', l,
     +			' (sub_readModel.f)'
	          B(l,0,k) = 0.0
	       end if  ! if (B(l,0,k).ne.0.0)

	    end do  ! end l loop
         end do  ! end k loop

      elseif (choiceD.eq.1) then  ! if 2D model

	 do l = 0, lmax
	    do m = 0, l
	   
	       read(10,*), tmpl, tmpm, A(l,m,0), B(l,m,0)

	       if (tmpl.ne.l) then  ! check
		  print*, '!!Error(sub_readModel): tmpl = ', tmpl,
     +			  '                        l = ', l
	       end if  ! if (tmpl.ne.l)
	       if (tmpm.ne.m) then  ! check
		  print*, '!!Error(sub_readModel): tmpm = ', tmpm,
     +			  '                        m = ', m
	       end if  ! if (tmpm.ne.m)

	       if (m.ne.0) then  ! converting to real coefficients
	          A(l,m,0) = 2.0*A(l,m,0)
		  B(l,m,0) = -2.0*B(l,m,0)
	       end if  ! if (m.ne.0)

	    end do  ! end m loop
	 end do  ! end l loop

      end if  ! if (choiceD.eq.2)...

      close(10)

      RETURN 
      END