#include "dims.h"

! am_02/02: calculate the coefficient stencil for both hemisphere which is used
!           to calculate the height integrated current densities K_(q,phi)K_(q,lam)

       subroutine nosocoef
!      
      use cons_module,only: pi
      implicit none
C     ****
C     ****     calculate coefficients for dynamo pde for both hemisphere
C     ****
#include "params.h"
#include "coefm.h"
#include "consdyn.h" ! dlatm,dlonm
! Global:
      common/nscoeff/ nscoef(imx0,jmaxm,10),nscrrt(imx0,jmaxm)
      real :: nscoef,nscrrt
! 
     
! Local:
      real :: array(0:imx0+1,jmaxm),cs(jmaxm)
      real :: nszigmc(imaxmp,jmaxm),nszigm2(imaxmp,jmaxm)
      real :: nszigm11(imaxmp,jmaxm),nszigm22(imaxmp,jmaxm)
      real :: nsrhs(imaxmp,jmaxm)
      real :: dfac,dfac1n,dfac1s,dfac2n,dfac2s
      integer :: j,je,jj,jjj,i,n
! Externals:
      real,external :: sddot ! in util.F
!
!   clear arrays
      nszigmc(:,:)  = 0.0 
      nszigm2(:,:)  = 0.0 
      nszigm11(:,:) = 0.0 
      nszigm22(:,:) = 0.0 
      nsrhs(:,:)    = 0.0 
      
!   calculate magnatic latitude cosinus array
!
      dlonm = 2.0*pi/imaxm
      dlatm = pi/(jmaxm-1)
!      
      do j = 1,jmaxm     ! -pi/2 to pi/2
        cs(j) = cos(-pi/2.+(j-1)*dlatm)
      enddo
!
!   reverse sign of ZIGMC to be compatible with Cicely's (richmond)
!   calculate difference
!   for zigmc & zigm2 no sign change because values from the
!                     coresponding hemisphere are used
!   zigmc : sum_{phi lam} = +/(-) (sum_H - sum_C) -> C Ridley
!   zigm2 : sum_{lam phi} = -/(+) (sum_H + sum_C) -> D Ridley
!   zigm11: sum_{phi phi}   -> A Ridley
!   zigm22: sum_{lam lam}   -> B Ridley
!   
!   factors from difference quotients and derivatives
!       4.*dlatm*dlonm for mixed terms zigmc and zigm2
!       dlatm**2 or dlonm**2 for zigm22 and zigm11
!   
!   factor cosin (cs and 1/cs) from pde for zigm22 and zigm11
!
      dfac = 4.*dlatm*dlonm
      do j = 2,jmaxm-1    ! 2,96  not value at the poles 
!	  
	dfac1n = cs(j)/dlatm**2
	dfac2n = cs(j)*dlonm**2
!	  
        do i = 1,imaxmp
	  nszigmc(i,j)  = -zigmc(i,j)
!	  
          nszigmc(i,j)  = (nszigmc(i,j)+zigm2(i,j))/dfac   
          nszigm2(i,j)  = nszigmc(i,j)-2.*zigm2(i,j)/dfac
          nszigm22(i,j) = zigm22(i,j)*dfac1n
          nszigm11(i,j) = zigm11(i,j)/dfac2n
        enddo
!	  
      enddo
!               
      j = (jmaxm+1)/2.0   ! change sign for values at equator maybe not necessary, but
                          ! then the sign for coefficient has to be changed too
      do i = 1,imaxmp
	nszigmc(i,j)  = -nszigmc(i,j)
	nszigm2(i,j)  = -nszigm2(i,j)
      enddo
!
!     values at the poles (1 and 97)
!
      jj  = jmaxm
      jjj = 1
      dfac1n = cs(jj)/dlatm**2
      dfac1s = cs(jjj)/dlatm**2   ! is not necessary cos symmetric
      
      do i = 1,imaxmp
!	  
	nszigmc(i,1)     = -zigmc(i,1)            ! 1
	nszigmc(i,jmaxm) = -zigmc(i,jmaxm)	  ! 97
!	
        nszigmc(i,jj)    = (nszigmc(i,jj)+zigm2(i,jj))/dfac   
        nszigm2(i,jj)    =  nszigmc(i,jj)-2.*zigm2(i,jj)/dfac
        nszigm22(i,jj)   =  zigm22(i,jj)*dfac1n
        nszigmc(i,jjj)   = (nszigmc(i,jjj)+zigm2(i,jjj))/dfac
        nszigm2(i,jjj)   =  nszigmc(i,jjj)-2.*zigm2(i,jjj)/dfac
        nszigm22(i,jjj)  =  zigm22(i,jjj)*dfac1s
!
!       set zigm11 to zero at the magnetic poles (1 and 97) to avoid floating point
!       exception (values at the poles are not used)
!
        nszigm11(i,jj)  = 0.0
        nszigm11(i,jjj) = 0.0
!
      enddo
     
!
!     clear array for difference stencil over north and south hemisphere 
!
      nscoef(:,:,:) = 0.0
  
!     calculate contribution to stencil from each pde coefficient
!     one at a time because of smaller working arrays 
!
      call nsstencil(nszigm11,imx0,jmaxm,nscoef,array,1,jmaxmh)
      call nsstencil(nszigm22,imx0,jmaxm,nscoef,array,4,jmaxmh)
      call nsstencil(nszigmc ,imx0,jmaxm,nscoef,array,2,jmaxmh)
      call nsstencil(nszigm2 ,imx0,jmaxm,nscoef,array,3,jmaxmh) 
!
!     set boundary conditions at pole 
!     value change from 1.0 to 0.5 for each hemisphere
!
      do i = 1,imx0
        do n = 1,8
          nscoef(i,jmaxm,n) = 0.
          nscoef(i,1,n)     = 0.
        enddo
        nscoef(i,jmaxm,9) = 0.5
        nscoef(i,1,9)     = 0.5
      enddo
!
!     divide stencil by cos(theta)
!
      do j = 2,jmaxm-1
        do n = 1,9
	  nscoef(:,j,n) = nscoef(:,j,n)/cs(j)
	enddo
      enddo 
! 
!     calculate right hand side of pde from rim(1) and rim(2)
!
      do 45 j = 2,jmaxmh-1  ! 2,48  south pole-1 to equator-1
        jj = j+jmaxmh-1     ! 50,96 equator-1 to north pole-1 
!  
!       differentiate rim(1) w.r.t lambda
!  
        do i = 2,imaxm-1
          nsrhs(i,j)  = 1.0/(dlonm*cs(j))*
     |                      0.5*(rim(i+1,j,1)-rim(i-1,j,1))
          nsrhs(i,jj) = 1.0/(dlonm*cs(jj))*
     |                      0.5*(rim(i+1,jj,1)-rim(i-1,jj,1))
	enddo
!  
!       values at the poles
!  
        nsrhs(1,j)  = 1.0/(dlonm*cs(j))*
     |                      0.5*(rim(2,j,1)-rim(imaxm,j,1))
        nsrhs(1,jj) = 1.0/(dlonm*cs(jj))*
     |                      0.5*(rim(2,jj,1)-rim(imaxm,jj,1))
        nsrhs(imaxm,j)  = 1.0/(dlonm*cs(j))*
     |                      0.5*(rim(1,j,1)-rim(imaxm-1,j,1))
        nsrhs(imaxm,jj) = 1.0/(dlonm*cs(jj))*
     |                      0.5*(rim(1,jj,1)-rim(imaxm-1,jj,1))
!  
   45 continue
!  
!       differentiate rim(2) w.r.t theta0
!  
      do j = 2,jmaxmh-1  ! 2,48  south pole -1 to equator-1
	jj = j+jmaxmh-1  ! 50,96 equator+1 to north pole-1 
	do i = 1,imaxm
	 nsrhs(i,j)  = nsrhs(i,j) - 1.0/(dlatm*cs(j))*0.5*
     |  	      (rim(i,j+1,2)*cs(j+1)-rim(i,j-1,2)*cs(j-1))
	 nsrhs(i,jj) = nsrhs(i,jj) + 1.0/(dlatm*cs(jj))*0.5*
     |  	      (rim(i,jj+1,2)*cs(jj+1)-rim(i,jj-1,2)*cs(jj-1))
       enddo
      enddo
!    
!     calculate value at the poles by averaging over i:imaxm
!    
      nsrhs(1,jmaxm) = -2./float(imaxm)*
     |                  sddot(imaxm,unit,rim(1,jmaxm-1,2))/cs(jmaxm-1)
      nsrhs(1,1)     = -2./float(imaxm)*
     |                  sddot(imaxm,unit,rim(1,2,2))/cs(2)
!
!     extend over longitude
!
      nsrhs(:,jmaxm) = nsrhs(1,jmaxm)
      nsrhs(:,1)     = nsrhs(1,1)
!
! note: for calculating J_mr values with the stencil at the equator not used
! note: for the test case tstjmrim when both hemisphere are added together
!        get double values at the equator, since no seperate value for south and
!        north of the equator, at the equator jump in rhs, therefore values
!        doubled at equator
! note: nsrhs stencil is the same as rhs in transf.f if average (north & south
!       of equator is taken for derivative in lam direction 
! note: introduced 0.5 to fit coefficient stencil = double nscoef(49)
!       for consistency also for nsrhs: c0(j=49,10) =  double nsrhs(49)
!
      je = jmaxmh
      i = 1
      nsrhs(i,je) = 0.5/dlonm*(rim(i+1,je,1)-rim(imaxm,je,1))
      nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |  		rim(i,je,2)+ cs(je+1)*rim(i,je+1,2))
      nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |  		rim(i,je,2)+ cs(je-1)*rim(i,je-1,2))
      nsrhs(i,je) = 0.5*nsrhs(i,je)
      do i = 2,imaxm-1
	
        nsrhs(i,je) = 0.5/dlonm*(rim(i+1,je,1)-rim(i-1,je,1))
	nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |                    rim(i,je,2)+ cs(je+1)*rim(i,je+1,2))
	nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |                    rim(i,je,2)+ cs(je-1)*rim(i,je-1,2))
        nsrhs(i,je) = 0.5*nsrhs(i,je)
      enddo
      
      i = imaxm
      nsrhs(i,je) = 0.5/dlonm*(rim(1,je,1)-rim(i-1,je,1))
      nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |  		rim(i,je,2)+ cs(je+1)*rim(i,je+1,2))
      nsrhs(i,je) = nsrhs(i,je) + 1./dlatm*(cs(je)*
     |  		rim(i,je,2)+ cs(je-1)*rim(i,je-1,2))
      nsrhs(i,je) = 0.5*nsrhs(i,je)
!
!     periodic points
!
      nsrhs(imaxmp,:) = nsrhs(1,:)
!
!     scale rhs by refernce radius (R_E + H0) in meters dfac = r0*1e-2
!
      dfac = r0*1.0e-2
      nsrhs(:,:) = nsrhs(:,:)*dfac
!
!     insert nsrhs into coefficient : from south to north pole
!     and divide by cos(theta) =1.0 not necessary !
!
      nscoef(:,:,10) = nsrhs(:,:)

!
!     set value of solution to 1.0 at poles ! change control later with equilibrium
!
      nscoef(:,jmaxm,10) = 0.5
      nscoef(:,1,10)	 = 0.5

      return
      end
!
!----------------------------------------------------------------------------
      subroutine nsstencil(zigm,imx0,jmaxm,nscoef,array,ncoef,jmaxmh)
      implicit none
!
! Args:
      integer,intent(in):: imx0,jmaxm,ncoef,jmaxmh
      real,intent(in)   :: zigm(imx0,jmaxm)
      real,intent(out)  :: nscoef(imx0,jmaxm,10),array(0:imx0+1,jmaxm)
!
! Local:
      integer :: i,j
!
!     copy zigm (south to north pole) into array (south to north pole)
!     check: is this change necessary?
! 
      do j = 1,jmaxm    ! 1,97
        do i = 1,imx0
	  array(i,j) = zigm(i,j)
        enddo
      enddo
!
!     extend over 32 grid spaces to alow a total of five grid levels
!     check: is this necessary because no solving with these coefficients?
!     
      i = 1
      do j = 1,jmaxm
        array(i-1,j)	= array(imx0-i,j)
        array(imx0+i,j) = array(1+i,j)
      enddo
      
!     calculate contribution to stencil for each grid point
!  
      call nscnm(array,imx0,jmaxm,nscoef,ncoef,jmaxmh)
!  
      return
      end
!
!----------------------------------------------------------------------------
      subroutine nscnm(array,imx0,jmaxm,nsc,ncoef,jmaxmh)
      implicit none
!      
!      calculate contribution for each zigm
!      the diagonial dominance of the stencil is not checked
!      since the coefficients are not used for solving
!      one reason for difference between added north-south
!      and seperated north south hemisphere
!
!      stencil for southern hemisphere changed, since
!      also for the southern hemisphere latudinal counts from
!      the equator to the south pole as for the northern hemisphere
!      from the equator to the north pole    
!      nsc(i,j,n) n:
!      northern hemisphere stencil n=(1,2,3,4,5,6,7,8,9)
!      southern hemisphere stencil n=(1,-8,-7,-5,-4,-2,9)
!
!      values at the equator (j=49): not separeted for northern and southern
!      hemisphere- only one stencil and later in nscrrt equally  
!      distributed to northern and southern hemisphere
!
! Args:
      integer,intent(in) :: imx0,jmaxm,ncoef,jmaxmh
      real,intent(in)    :: array(0:imx0+1,jmaxm)
      real,intent(out)   :: nsc(imx0,jmaxm,10)
!
! Local:
      integer :: i,j,jj,jjj
!
!     check: what's about boundary equator values for nsc, why not used?
!
!     calculate contribution for zigm11
!
      if(ncoef.eq.1)then	
!      
        do j = 2,jmaxmh-1   ! 2,48  not value at the pols (jj=97 jjj=1)
          jj  = jmaxmh+j-1      ! 50,96 northern henisphere
          jjj = jmaxmh-j+1      ! 48,2  southern henisphere
!
          do i = 1,imx0
	    nsc(i,jj,1) = nsc(i,jj,1)+.5*(array(i,jj)+array(i+1,jj))
	    nsc(i,jj,5) = nsc(i,jj,5)+.5*(array(i,jj)+array(i-1,jj))
	    nsc(i,jj,9) = nsc(i,jj,9)-.5*(array(i+1,jj)+2.*array(i,jj)+
     |                                 array(i-1,jj))
     
	    nsc(i,jjj,1) = nsc(i,jjj,1)+.5*(array(i,jjj)+array(i+1,jjj))
	    nsc(i,jjj,5) = nsc(i,jjj,5)+.5*(array(i,jjj)+array(i-1,jjj))
	    nsc(i,jjj,9) = nsc(i,jjj,9)-.5*(array(i+1,jjj)+
     |                                 2.*array(i,jjj)+array(i-1,jjj))
          enddo
	enddo

! am 2001-6-27 include boundary condition at equator
	
        jj  = jmaxmh      ! 49
!
        do i = 1,imx0
	  nsc(i,jj,1) = nsc(i,jj,1)+.5*(array(i,jj)+array(i+1,jj))
	  nsc(i,jj,5) = nsc(i,jj,5)+.5*(array(i,jj)+array(i-1,jj))
	  nsc(i,jj,9) = nsc(i,jj,9)-.5*(array(i+1,jj)+2.*array(i,jj)+
     |  			     array(i-1,jj))
        enddo
!
!     calculate contribution for zigm12 (=ZIGMC+ZIGM2)
! 
      else if(ncoef.eq.2) then
!	
        do j = 2,jmaxmh-1       ! 2,48  not value at the pols (jj=97 jjj=1)
          jj  = jmaxmh+j-1      ! 50,96 northern henisphere
          jjj = jmaxmh-j+1      ! 48,2  southern henisphere
!
          do i = 1,imx0 
	    nsc(i,jj,2) = nsc(i,jj,2)+.5*(array(i,jj)+array(i+1,jj))
	    nsc(i,jj,4) = nsc(i,jj,4)-.5*(array(i,jj)+array(i-1,jj))
	    nsc(i,jj,6) = nsc(i,jj,6)+.5*(array(i,jj)+array(i-1,jj))
	    nsc(i,jj,8) = nsc(i,jj,8)-.5*(array(i,jj)+array(i+1,jj))
	    nsc(i,jj,3) = nsc(i,jj,3)+.5*(-array(i-1,jj)+array(i+1,jj))
	    nsc(i,jj,7) = nsc(i,jj,7)-.5*(-array(i-1,jj)+array(i+1,jj))
	    
	    nsc(i,jjj,2)=nsc(i,jjj,2)+.5*(array(i,jjj)+array(i+1,jjj))
	    nsc(i,jjj,4)=nsc(i,jjj,4)-.5*(array(i,jjj)+array(i-1,jjj))
	    nsc(i,jjj,6)=nsc(i,jjj,6)+.5*(array(i,jjj)+array(i-1,jjj))
	    nsc(i,jjj,8)=nsc(i,jjj,8)-.5*(array(i,jjj)+array(i+1,jjj))
	    nsc(i,jjj,3)=nsc(i,jjj,3)+.5*(-array(i-1,jjj)+array(i+1,jjj))
	    nsc(i,jjj,7)=nsc(i,jjj,7)-.5*(-array(i-1,jjj)+array(i+1,jjj))
           enddo
	enddo
!
!     calculate contribution for zigm21 (=ZIGMC-ZIGM2)
! 
      else if(ncoef.eq.3) then
!
        do j = 3,jmaxmh-1   ! 3,48  not value at the pols (jj=97 jjj=1)
          jj  = jmaxmh+j-1      ! 51,96 northern hemisphere
          jjj = jmaxmh-j+1      ! 47,2  southern hemisphere
!
          do i = 1,imx0 
!                           
	  nsc(i,jj,2) = nsc(i,jj,2)+.5*(array(i,jj)+array(i,jj+1))
	  nsc(i,jj,4) = nsc(i,jj,4)-.5*(array(i,jj)+array(i,jj+1))
	  nsc(i,jj,6) = nsc(i,jj,6)+.5*(array(i,jj)+array(i,jj-1))
	  nsc(i,jj,8) = nsc(i,jj,8)-.5*(array(i,jj)+array(i,jj-1))
	  nsc(i,jj,1) = nsc(i,jj,1)+.5*(array(i,jj+1)-array(i,jj-1))
	  nsc(i,jj,5) = nsc(i,jj,5)-.5*(array(i,jj+1)-array(i,jj-1))
!
	  nsc(i,jjj,8)=nsc(i,jjj,8)-.5*(array(i,jjj)+array(i,jjj+1))
	  nsc(i,jjj,6)=nsc(i,jjj,6)+.5*(array(i,jjj)+array(i,jjj+1))
	  nsc(i,jjj,4)=nsc(i,jjj,4)-.5*(array(i,jjj)+array(i,jjj-1))
	  nsc(i,jjj,2)=nsc(i,jjj,2)+.5*(array(i,jjj)+array(i,jjj-1))
	  nsc(i,jjj,1)=nsc(i,jjj,1)-.5*(array(i,jjj+1)-array(i,jjj-1))
	  nsc(i,jjj,5)=nsc(i,jjj,5)+.5*(array(i,jjj+1)-array(i,jjj-1))
          enddo
	enddo
	
	j = 2	! 2 change in sign for equatorial values 
	jj  = jmaxmh+j-1      ! 50 northern hemisphere
	jjj = jmaxmh-j+1      ! 48  southern hemisphere
       
       do i = 1,imx0   !am2001-7-3
			
	nsc(i,jj,2) = nsc(i,jj,2)+.5*(array(i,jj)+array(i,jj+1))
	nsc(i,jj,4) = nsc(i,jj,4)-.5*(array(i,jj)+array(i,jj+1))
	nsc(i,jj,6) = nsc(i,jj,6)+.5*(array(i,jj)-array(i,jj-1))
	nsc(i,jj,8) = nsc(i,jj,8)-.5*(array(i,jj)-array(i,jj-1))
	nsc(i,jj,1) = nsc(i,jj,1)+.5*(array(i,jj+1)+array(i,jj-1))
	nsc(i,jj,5) = nsc(i,jj,5)-.5*(array(i,jj+1)+array(i,jj-1))

	nsc(i,jjj,8)=nsc(i,jjj,8)-.5*(array(i,jjj)-array(i,jjj+1))
	nsc(i,jjj,6)=nsc(i,jjj,6)+.5*(array(i,jjj)-array(i,jjj+1))
	nsc(i,jjj,4)=nsc(i,jjj,4)-.5*(array(i,jjj)+array(i,jjj-1))
	nsc(i,jjj,2)=nsc(i,jjj,2)+.5*(array(i,jjj)+array(i,jjj-1))
	nsc(i,jjj,1)=nsc(i,jjj,1)-.5*(-array(i,jjj+1)-array(i,jjj-1))
	nsc(i,jjj,5)=nsc(i,jjj,5)+.5*(-array(i,jjj+1)-array(i,jjj-1))
       
       enddo
!    
!     low latitude boundary conditions: contribution from zigm21(i,j-1/2)=0
!
	j = 1
        jj  = jmaxmh+j-1      ! 49
        jjj = jmaxmh-j+1      ! 49
!
	do i = 1,imx0 
	  nsc(i,jj,2) = nsc(i,jj,2)+.25*(-array(i,jj)+array(i,jj+1))
	  nsc(i,jj,4) = nsc(i,jj,4)-.25*(-array(i,jj)+array(i,jj+1))
	  nsc(i,jj,1) = nsc(i,jj,1)+.25*(-array(i,jj)+array(i,jj+1))
	  nsc(i,jj,5) = nsc(i,jj,5)-.25*(-array(i,jj)+array(i,jj+1))
	
	  nsc(i,jj,2) = nsc(i,jj,2)+.25*(-array(i,jjj)+array(i,jjj-1))
	  nsc(i,jj,4) = nsc(i,jj,4)-.25*(-array(i,jjj)+array(i,jjj-1))
	  nsc(i,jj,1) = nsc(i,jj,1)-.25*(array(i,jjj)-array(i,jjj-1))
	  nsc(i,jj,5) = nsc(i,jj,5)+.25*(array(i,jjj)-array(i,jjj-1))
	enddo
!
!     calculate contribution for zigm22 
! 
      else if(ncoef.eq.4) then
!
        do j = 2,jmaxmh-1   ! 2,48  not value at the pols (jj=97 jjj=1)
          jj  = jmaxmh+j-1      ! 50,96 northern henisphere
          jjj = jmaxmh-j+1      ! 48,2  southern henisphere
!
          do i = 1,imx0 
	    nsc(i,jj,3) = nsc(i,jj,3)+.5*(array(i,jj)+array(i,jj+1))
	    nsc(i,jj,7) = nsc(i,jj,7)+.5*(array(i,jj)+array(i,jj-1))
	    nsc(i,jj,9) = nsc(i,jj,9)-.5*(array(i,jj-1)+2.0*array(i,jj)+
     |                                 array(i,jj+1))
 !    
	    nsc(i,jjj,7)=nsc(i,jjj,7)+.5*(array(i,jjj)+array(i,jjj+1))
	    nsc(i,jjj,3)=nsc(i,jjj,3)+.5*(array(i,jjj)+array(i,jjj-1))
	    nsc(i,jjj,9)=nsc(i,jjj,9)-.5*(array(i,jjj-1)+
     |                                 2.0*array(i,jjj)+ array(i,jjj+1))
           enddo
	enddo
!    
!     low latitude boundary conditions: contribution from zigm22(i,j-1/2)=0
!
	j = 1
        jj  = jmaxmh+j-1      ! 49
        jjj = jmaxmh-j+1      ! 49
	do i = 1,imx0 
	   nsc(i,jj,3) = nsc(i,jj,3)+.25*(array(i,jj)+array(i,jj+1))
	   nsc(i,jj,9) = nsc(i,jj,9)-.25*(array(i,jj)+array(i,jj+1))
!  am 2001-7-02 otherwise double coefficients
	   nsc(i,jj,3)=nsc(i,jj,3)+.25*(array(i,jjj)+array(i,jjj-1))
	   nsc(i,jj,9)=nsc(i,jj,9)-.25*(array(i,jjj-1)+array(i,jjj))
       enddo
!
      endif
!    
      return
      end