!
      module nc_module
!      
      use param_module,only: mlat,mlon,rmlat,rmlon,rglat,rglon
!      
      implicit none
!      
      integer,parameter :: count3(3)= (/mlon,mlat,1/),
     |     rm_count3(3)= (/rmlon,rmlat,1/),
     |     rg_count3(3)= (/rglon,rglat,1/)
      integer ::  start3(3),rm_start3(3)
      integer :: noid,	! id of output netcdf-file
     |   rm_dim3(3),	! id of  defining 3dim-arrays on regular mag. grid
     |   rmlt_dim3(3),	! id of  defining 3dim-arrays on regular mag.local time grid
     |   rg_dim3(3),	! id of  defining 3dim-arrays on regular geo. grid
     |   dim3(3)	! id of  defining 3dim-arrays on irregular grid
!      
      contains
!- -------------------------------------------------------------------------
      subroutine ncdefpdim

! nc-output: define plot dimensions & attributes
      use param_module,only: rmlat,rmlon,rglat,rglon
      use eqvcur_module,only: vallat,vallon,
     |     reglat,reglon,regmlt,reg_glat,reg_glon
!
      implicit none
!
#ifdef SUN
#include "/opt/local/include/netcdf.inc"
#else
#include "netcdf.inc"
#endif
!
      integer :: istat,ilon,ilat,imlt,idtime,idlat,idlon,idmlt,i,j,
     |   iglon,iglat,idglat,idglon
      real :: delta

!  put the longitudes & latitudes in output nc-file 
      istat = nf_redef(noid)
      if (istat /= NF_NOERR) call check_err(istat,'redef ')

! magnetic grid
      istat = nf_def_dim(noid,'rmlon',rmlon,idlon)
      if (istat /= NF_NOERR) call check_err(istat,'defdim rmlon ')

      istat = nf_def_dim(noid,'rmlat',rmlat,idlat)
      if (istat /= NF_NOERR) call check_err(istat,'defdim rmlat ')
      
      istat = nf_def_dim(noid,'rmlt',rmlon,idmlt)
      if (istat /= NF_NOERR) call check_err(istat,'defdim rmlt ')

      istat = nf_def_var(noid,'rmlon',NF_REAL,1,idlon,ilon)
      if (istat /= NF_NOERR) call check_err(istat,'defvar rmlon ')

      istat = nf_def_var(noid,'rmlat',NF_REAL,1,idlat,ilat)
      if (istat /= NF_NOERR) call check_err(istat,'defvar rmlat ')

      istat = nf_def_var(noid,'rmlt',NF_REAL,1,idmlt,imlt)
      if (istat /= NF_NOERR) call check_err(istat,'defvar rmlt ')

      istat = nf_put_att_text(noid,ilon,'long_name',15,
     |   'magn. longitude')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlon name')

      istat = nf_put_att_text(noid,ilat,'long_name',14,
     |   'magn. latitude')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlat name')
      
      istat = nf_put_att_text(noid,imlt,'long_name',16,
     |   'magn. local time')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlt name')

      istat = nf_put_att_text(noid,ilon,'units',7,'degrees')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlon unit')

      istat = nf_put_att_text(noid,ilat,'units',7,'degrees')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlat unit')

      istat = nf_put_att_text(noid,imlt,'units',3,'hrs')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rmlt unit')
      
! geographic grid
      istat = nf_def_dim(noid,'rglon',rglon,idglon)
      if (istat /= NF_NOERR) call check_err(istat,'defdim rglon ')

      istat = nf_def_dim(noid,'rglat',rglat,idglat)
      if (istat /= NF_NOERR) call check_err(istat,'defdim rglat ')

      istat = nf_def_var(noid,'rglon',NF_REAL,1,idglon,iglon)
      if (istat /= NF_NOERR) call check_err(istat,'defvar rglon ')

      istat = nf_def_var(noid,'rglat',NF_REAL,1,idglat,iglat)
      if (istat /= NF_NOERR) call check_err(istat,'defvar rglat ')

      istat = nf_put_att_text(noid,iglon,'long_name',14,
     |   'geo. longitude')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rglon unit')

      istat = nf_put_att_text(noid,iglat,'long_name',13,
     |   'geo. latitude')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rglat unit')

      istat = nf_put_att_text(noid,iglon,'units',7,'degrees')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rglon unit')

      istat = nf_put_att_text(noid,iglat,'units',7,'degrees')
      if (istat /= NF_NOERR) call check_err(istat,'put att_rglat unit')

      istat = nf_enddef(noid)
      if (istat /= NF_NOERR) call check_err(istat,'enddef ') 

! put magnetic coordinates to nc-file
! set rmlat regular latitude grid 
      delta = 360./float(rmlon-1)
      do j = 1,rmlon
	 reglon(j) = float(j-1)*delta -180.
      enddo
      
!      reglon = vallon   ! set rmlon regular longitude grid (= vallon)
      
! put mag longitude in new file
      istat = nf_put_var_double(noid,ilon,reglon)
      if (istat /= NF_NOERR) call check_err(istat,'putvar reglon ')  

! set rmlat regular latitude grid 
      delta = 180./float(rmlat-1)
      do j = 1,rmlat
	 reglat(j) = float(j-1)*delta -90.
      enddo
!      reglat = vallat ! ptest
      
! put mag. latitude in new file
      istat = nf_put_var_double(noid,ilat,reglat)
      if (istat /= NF_NOERR) call check_err(istat,'putvar reglat ')

! set rmlt regular magnetic local time
      delta = 24./float(rmlon-1)
      do j = 1,rmlon
	 regmlt(j) = float(j-1)*delta 
      enddo
      
! put mag. local time in new file
      istat = nf_put_var_double(noid,imlt,regmlt)
      if (istat /= NF_NOERR) call check_err(istat,'putvar rmltt ')

! put geographic coordinates to nc-file

! set rglon regular longitude grid 
      delta = 360./float(rglon-1)
      do i = 1,rglon
        reg_glon(i) = delta*float(i-1) - 180.	! set regular geo. grid
      enddo
      
! put geo. longitude in new file
      istat = nf_put_var_double(noid,iglon,reg_glon)
      if (istat /= NF_NOERR) call check_err(istat,'putvar reg_glon ')  

! set rglat regular latitude grid (= vallat)
      delta = 180./float(rglat-1)
      do j = 1,rglat
	  reg_glat(j) = delta*float(j-1) -90.
      enddo
      
! put mag. latitude in new file
      istat = nf_put_var_double(noid,iglat,reg_glat)
      if (istat /= NF_NOERR) call check_err(istat,'putvar reg_glat ')

! get time id      
      istat = nf_inq_varid(noid,'time',idtime)
      if (istat /= NF_NOERR) call check_err(istat,'id time ')

      rm_dim3(1) = idlon
      rm_dim3(2) = idlat
      rm_dim3(3) = idtime

      rmlt_dim3(1) = idmlt
      rmlt_dim3(2) = idlat
      rmlt_dim3(3) = idtime
      
      rg_dim3(1) = idglon
      rg_dim3(2) = idglat
      rg_dim3(3) = idtime

      end subroutine ncdefpdim
!----------------------------------------------------------------------------
      subroutine output
      
      use param_module,only: mlon,mlat,igeoplot,
     |    rmlon,rmlat,rglon,rglat
      use eqvcur_module
      use map_module
      use coef_module,only: nm,mm,c_fld,c_tor,c_tot,c_grd,
     |   c_mptext,c_mpttotgrd ,c_mptextgrd,c_mptintgrd
      
      implicit none
      
      integer :: nw,nw_u,i,j
      character :: fname*13,labl*56,units*12
      real :: cvt_geo(rglon,rglat),psi(rmlon,rmlat)
      
!   on irregular grid
   	
      fname = 'kqlam_org'
      labl = 'K_qlam(org)'
      units = 'A/m'
      call ncplot(noid,fname,labl,start3,count3,
     |  	  dim3,mlon,mlat,kqlam_org,11,units,3)
     
      fname = 'kqphi_org'
      labl = 'K_qphi(org)'
      units = 'A/m'
      call ncplot(noid,fname,labl,start3,count3,
     |  	  dim3,mlon,mlat,kqphi_org,11,units,3)
     
!   on regular grid
     
      fname = 'kqlam'
      labl = 'K_qlam'
      units = 'A/m'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	  rm_dim3,rmlon,rmlat,kqlam,6,units,3)
     
      fname = 'kqphi'
      labl = 'K_qphi'
      units = 'A/m'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	  rm_dim3,rmlon,rmlat,kqphi,6,units,3)

! equivalent current function (toroidal part)
      fname = 'eqcf_tor'
      labl = 'eqcf toroidal ionospheric part'
      call plt(c_tor,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,40,units,3)
      
      fname = 'eqcf_tor_asym'
      labl = 'eqcf antisymmetric toroidal part'
      call plt(c_tor,1,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,32,units,3)
     
      fname = 'eqcf_tor_sym'
      labl = 'eqcf symmetric toroidal part'
      call plt(c_tor,-1,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,29,units,3)
     
      fname = 'eqc_tot'
      labl = 'eqcf total'
      call plt(c_tot,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,10,units,3)
     
      fname = 'eqc_fldalg'
      labl = 'eqcf field-aligned part'
      call plt(c_fld,-1,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,23,units,3)
     
      fname = 'eqc_tot_sym'
      labl = 'eqcf total symmetric'
      call plt(c_tot,-1,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,20,units,3)
     
      fname = 'eqc_grd'
      labl = 'eqcf ground'
      call plt(c_grd,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,11,units,3)
     
      fname = 'magpot_ext'
      labl = 'magn. potential external'
      call plt(c_mptext,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,24,units,3)
     
      fname = 'magpot_int'
      labl = 'magn. potential internal @ground'
      call plt(c_mptintgrd,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	rmlt_dim3,rmlon,rmlat,psi,24,units,3)
     
      fname = 'magpot_tot'
      labl = 'magpot total @ground'
      call plt(c_mpttotgrd,0,psi)
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rmlt_dim3,rmlon,rmlat,psi,12,units,3)
     
      fname = 'dB_lam'
      labl = 'dB_lam (northward)'
      units = 'nT'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	rm_dim3,rmlon,rmlat,dB_lam,18,units,2)
     
      fname = 'dB_phi'
      labl = 'dB_phi (eastward)'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	rm_dim3,rmlon,rmlat,dB_phi,17,units,2)
!     
      fname = 'dB_dwn'
      labl = 'dB_dwn (downward)'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rm_dim3,rmlon,rmlat,dB_dwn,17,units,2)
     
      fname = 'B_abs_mag'
      labl = '|dB(mag)|'
      call ncplot(noid,fname,labl,rm_start3,rm_count3,
     |  	 rm_dim3,rmlon,rmlat,babs_mag,9,units,2)
     
 !  
      if(igeoplot == 1) then
      
        fname = 'eqc_tot_geo'
        labl = 'eqcf (geo.)'
        units = 'A/m'
        call plt_mlon(c_tot,0,psi)
        call mag2geo(psi(1,1),cvt_geo(1,1),im(1,1),jm(1,1),
     |     dim(1,1),djm(1,1),rglon,rmlon,rglon,rglat,rmlon,rmlat)
        call ncplot(noid,fname,labl,rm_start3,rg_count3,
     |  	   rg_dim3,rglon,rglat,cvt_geo,11,units,1)
      
        fname = 'dB_lam_geo'
        labl = 'dB_lam (northward geo.)'
        units = 'nT'
        call mag2geo(dB_lam_geo(1,1),cvt_geo(1,1),im(1,1),jm(1,1),
     |     dim(1,1),djm(1,1),rglon,rmlon,rglon,rglat,rmlon,rmlat)
        call ncplot(noid,fname,labl,rm_start3,rg_count3,
     |  	   rg_dim3,rglon,rglat,cvt_geo,23,units,2)
     
        fname = 'dB_phi_geo'
        labl = 'dB_phi (eastward geo.)'
        call mag2geo(dB_phi_geo(1,1),cvt_geo(1,1),im(1,1),jm(1,1),
     |     dim(1,1),djm(1,1),rglon,rmlon,rglon,rglat,rmlon,rmlat)
        call ncplot(noid,fname,labl,rm_start3,rg_count3,
     |  	   rg_dim3,rglon,rglat,cvt_geo,23,units,2)
!
        fname = 'dB_dwn_geo'
        labl = 'dB_dwn (downward geo.)'
        call mag2geo(dB_dwn_geo(1,1),cvt_geo(1,1),im(1,1),jm(1,1),
     |     dim(1,1),djm(1,1),rglon,rmlon,rglon,rglat,rmlon,rmlat)
        call ncplot(noid,fname,labl,rm_start3,rg_count3,
     |  	    rg_dim3,rglon,rglat,cvt_geo,23,units,2)
     
        fname = 'B_abs_geo'
        labl = 'b_0*dB(geo)'
        call mag2geo(babs_geo(1,1),cvt_geo(1,1),im(1,1),jm(1,1),
     |     dim(1,1),djm(1,1),rglon,rmlon,rglon,rglat,rmlon,rmlat)
        call ncplot(noid,fname,labl,rm_start3,rg_count3,
     |  	 rg_dim3,rglon,rglat,cvt_geo,11,units,2)
 
      endif
     
      
      return
      end subroutine output
!- -------------------------------------------------------------------------
      end module nc_module
      
!- -------------------------------------------------------------------------
      subroutine ncplot(noid,fname,labl,start3,count3,dim3,nlmb,
     |                  nthplt,psi,nw,units,nw_u)

! nc-output: define plot variable & attributes
!
      implicit none
!
#ifdef SUN
#include "/opt/local/include/netcdf.inc"
#else
#include "netcdf.inc"
#endif
!      
      integer, intent(in) :: noid,nlmb,nthplt,start3(3),count3(3),nw
      integer, intent(in) :: dim3(3),nw_u
      real, intent(in) :: psi
      character, intent(in) :: fname*13,labl*56,units*12

      integer :: istat,id
      character :: msgerr*19

      
! check if field is already defined
      istat = nf_inq_varid(noid,fname,id)
      if (istat /= NF_NOERR) then
	
        istat = nf_redef(noid)
        if (istat /= NF_NOERR) call check_err(istat,'redef ')

        istat = nf_def_var(noid,fname,NF_REAL,3,dim3,id)
        msgerr = 'defvar '//fname
        if (istat /= NF_NOERR) call check_err(istat,msgerr)

        istat = nf_put_att_text(noid,id,'long_name',nw,labl)
        msgerr = 'put att_'//fname//' name'
        if (istat /= NF_NOERR) call check_err(istat,msgerr)
	
        istat = nf_put_att_text(noid,id,'units',nw_u,units)
        msgerr = 'put att_'//fname//' units'
        if (istat /= NF_NOERR) call check_err(istat,msgerr)

        istat = nf_enddef(noid)
        if (istat /= NF_NOERR) call check_err(istat,'enddef ')
	 
      endif
      
! put values in new file
      msgerr = 'putvar '//fname
      istat = nf_put_vara_double(noid,id,start3,count3,psi)
      if (istat /= NF_NOERR) call check_err(istat,msgerr) 

      end subroutine ncplot
!----------------------------------------------------------------------------
      SUBROUTINE PLT(C,ISM,psi)                                               
C Constructs gridded function PSI (in COMMON/PSIARR/) from coefficients C
C  ISM=-1 GIVES SYMMETRIC COMPONENT                                             
C  ISM=0  GIVES TOTAL                                                           
C  ISM=1  GIVES ANTISYMMETRIC COMPONENT                                         
C
C This subroutine originally called another routine to produce a contour plot 
C  of PSI.  Now it is set up for testing, merely calling OUTPSI to print out
C  a few values from the array.
C To activate this subroutine, remove the RETURN that immediately follows
C  the DIMENSION statements.

      use coef_module,only: fcmp,tmpo,mmp,nm,mm,nmp,
     |   p,
     |   r,pmopmmo,	! common/eqcr/
     |   f		! common/far/
      use param_module,only: rad,rmlon,rmlat
      use eqvcur_module,only: reglat,regmlt,reglon
!      
      implicit none
!      
      integer, intent(in) :: ism
      real, intent(in)  :: c(0:nm,-mm:mm) 
      real, intent(out) :: psi(rmlon,rmlat)     

      integer :: i,ii,j,jmx,jcon,l,n,m,mp,np,rmlat2
      real :: ph,th,ct,st,pm2,x                            
                     
      DO 100 J=1,rmlat 	
         jcon = rmlat + 1 -j						   
	 TH = (90. + reglat(jcon))/RAD 	! colatitude
         CT = COS(TH)   
         ST = SIN(TH)  
         P(0,0) = 1.
         DO 80 MP=1,MMP    
            M = MP - 1  		       
            IF (MP.NE.1) P(M,M) = PMOPMMO(M)*P(M-1,M-1)*ST
            PM2 = 0.	 
            DO 80 N=MP,NM	      
               NP = N + 1
               P(N,M) = (CT*P(N-1,M) -R(N-1,M)*PM2)/R(N,M)
               PM2 = P(N-1,M)
   80 	 continue 
                                                            
         DO 100 L=1,rmlon               
            PH = regmlt(l)*15./rad	! mag.longitude [rad] 
            CALL FCMP(PH,MM)
	    psi(l,j)    = 0.
	    
            DO 90 m = -mm,mm		!Loop over order M  
               MP = 1 + IABS(M)                                                       
               DO 90 NP=MP,NMP 		!Loop over degree N     
                  N = NP - 1                                                            
                  psi(L,j) = psi(L,j) + c(n,m)*P(N,MP-1)*F(M)
   90       CONTINUE   
  100 continue
  
      rmlat2 = (rmlat+1)/2
      if(ism == -1) then	! symmetric
        do i = 1,rmlon
	  do j = 1,rmlat2
            jcon = rmlat + 1 -j	
            x = psi(i,j) + psi(i,jcon)
	    x = x /2.
	    psi(i,j)    = x
	    psi(i,jcon) = x
	  enddo
	enddo
      elseif(ism == 1) then	! antisymmetric
        do i = 1,rmlon
	  do j = 1,rmlat2
            jcon = rmlat + 1 -j	
            x = psi(i,j) + psi(i,jcon)
	    x = x /2.
	    psi(i,j)    = psi(i,j) - x
	    psi(i,jcon) = psi(i,jcon)- x
	  enddo
	enddo
      endif

      return                                                                   
      end subroutine plt
 
!----------------------------------------------------------------------------
      SUBROUTINE PLT_mlon(C,ISM,psi)                                               
C Constructs gridded function PSI (in COMMON/PSIARR/) from coefficients C
C  ISM=-1 GIVES SYMMETRIC COMPONENT                                             
C  ISM=0  GIVES TOTAL                                                           
C  ISM=1  GIVES ANTISYMMETRIC COMPONENT                                         
C
C This subroutine originally called another routine to produce a contour plot 
C  of PSI.  Now it is set up for testing, merely calling OUTPSI to print out
C  a few values from the array.
C To activate this subroutine, remove the RETURN that immediately follows
C  the DIMENSION statements.

      use coef_module,only: fcmp,tmpo,mmp,nm,mm,nmp,
     |   p,
     |   r,pmopmmo,	! common/eqcr/
     |   f		! common/far/
      use param_module,only: rad,rmlon,rmlat
      use eqvcur_module,only: reglat,regmlt,reglon
!      
      implicit none
      
      real :: xmlon
      common/ssolar/ xmlon
!      
      integer, intent(in) :: ism
      real, intent(in)  :: c(0:nm,-mm:mm) 
      real, intent(out) :: psi(rmlon,rmlat)     

      integer :: i,ii,j,jmx,jcon,l,n,m,mp,np,rmlat2
      real :: ph,th,ct,st,pm2,x                            
                     
      DO 100 J=1,rmlat 	
         jcon = rmlat + 1 -j						   
	 TH = (90. + reglat(jcon))/RAD 	! colatitude
         CT = COS(TH)   
         ST = SIN(TH)  
         P(0,0) = 1.
         DO 80 MP=1,MMP    
            M = MP - 1  		       
            IF (MP.NE.1) P(M,M) = PMOPMMO(M)*P(M-1,M-1)*ST
            PM2 = 0.	 
            DO 80 N=MP,NM	      
               NP = N + 1
               P(N,M) = (CT*P(N-1,M) -R(N-1,M)*PM2)/R(N,M)
               PM2 = P(N-1,M)
   80 	 continue 
                                                            
         DO 100 L=1,rmlon      
            PH = (reglon(l)+180.-xmlon)/RAD        ! use magn. longitude             
!            PH = reglon(l)/rad	! mag.longitude [rad] 
            CALL FCMP(PH,MM)
	    psi(l,j)    = 0.
	    
            DO 90 m = -mm,mm		!Loop over order M  
               MP = 1 + IABS(M)                                                       
               DO 90 NP=MP,NMP 		!Loop over degree N     
                  N = NP - 1                                                            
                  psi(L,j) = psi(L,j) + c(n,m)*P(N,MP-1)*F(M)
   90       CONTINUE   
  100 continue
  
      rmlat2 = (rmlat+1)/2
      if(ism == -1) then	! symmetric
        do i = 1,rmlon
	  do j = 1,rmlat2
            jcon = rmlat + 1 -j	
            x = psi(i,j) + psi(i,jcon)
	    x = x /2.
	    psi(i,j)    = x
	    psi(i,jcon) = x
	  enddo
	enddo
      elseif(ism == 1) then	! antisymmetric
        do i = 1,rmlon
	  do j = 1,rmlat2
            jcon = rmlat + 1 -j	
            x = psi(i,j) + psi(i,jcon)
	    x = x /2.
	    psi(i,j)    = psi(i,j) - x
	    psi(i,jcon) = psi(i,jcon)- x
	  enddo
	enddo
      endif

      return                                                                   
      end subroutine plt_mlon                                                                     
