PROGRAM Convert_rcm_to_tecplot
!     USE Rcm_variables
!     USE Rcm_io
USE Rcm_mod_subs
      IMPLICIT NONE
INTEGER, PARAMETER :: N_GC=0
!     include 'mpif90.h'
!
      INTEGER :: itime = -999, time_min = -999, record_number = -1
!
      CHARACTER (LEN=24) :: filename_out_1d = "rcm-output-1D-xxxxxx.dat", &
                            filename_out_2d = "rcm-output-2D-xxxxxx.dat", &
                            filename_out_3d = "rcm-output-3D-xxxxxx.dat"
      CHARACTER (LEN=10) :: time_string = "xxxxxxxxxx", run_id_string
      CHARACTER (LEN=12) :: char_date, char_time
      CHARACTER (LEN=05) :: char_kvalue='xxxxx'
!
      INTEGER :: i,j,k, k_zone, rec_1, rec_2, k_selected=0, k_start, k_stop
      REAL (rprec), DIMENSION (1-n_gc:isize+n_gc, 1-n_gc:jsize+n_gc):: RCM_p,&
          RCM_n, RCM_T, RCM_PVgamma, MHD_p, veff, &
          temperature=0.0, density=0.0, birk_mhd=0.0, &
          RCm_p_40kv = 0.0  
      LOGICAL :: do_3d
!      
!
      WRITE (*,'(A)',ADVANCE='NO') 'ENTER RUN LABEL (10 CHARS MAX):__'
      READ (*,*) run_id_string
      WRITE (*,'(A)',ADVANCE='NO') 'ENTER RCM RECORDS TO START AND STOP, do_3d: '
      READ (*,*) rec_1, rec_2, do_3d
      IF (do_3d) THEN
          WRITE (*,'(A)',ADVANCE='NO') 'ENTER K-CHANNEL TO WRITE (-1 FOR ALL): '
          READ (*,*) k_selected
      END IF
!
!
      CALL Read_grid ()
      CALL Read_plasma ()
      CALL Read_qtcond ()
!
!
   DO record_number = rec_1, rec_2
!
      CALL Read_array ('rcmbndloc', record_number, label, ARRAY_1D = bndloc); imin_j = CEILING(bndloc)
!
      CALL Read_array ('rcmxmin'  , record_number, label, ARRAY_2D = xmin  )
      CALL Read_array ('rcmymin'  , record_number, label, ARRAY_2D = ymin  )
      CALL Read_array ('rcmvm'    , record_number, label, ARRAY_2D = vm    )
      CALL Read_array ('rcmbmin'  , record_number, label, ARRAY_2D = bmin  )
!
      CALL Read_array ('rcmv'     , record_number, label, ARRAY_2D = v     )
      CALL Read_array ('rcmbirk'  , record_number, label, ARRAY_2D = birk  )
      CALL Read_array ('rcmpedlam', record_number, label, ARRAY_2D = pedlam)
      CALL Read_array ('rcmpedpsi', record_number, label, ARRAY_2D = pedpsi)
      CALL Read_array ('rcmhall'  , record_number, label, ARRAY_2D = hall  )
      CALL Read_array ('rcmssequat', record_number, label, ARRAY_1D = ss  )
      CALL Read_array ('rcmvpar'  , record_number, label, ARRAY_2D = vpar  )
!!    CALL Read_array ('rcmbirkmhd',record_number, label, ARRAY_2D = birk_mhd)
!!      CALL Read_array ('rcmtemperature',record_number, label, ARRAY_2D = temperature)
!!      CALL Read_array ('rcmdensity',record_number, label, ARRAY_2D = density)
!
      CALL Read_array ('rcmeeta'  , record_number, label, ARRAY_3D = eeta  )
!      print*,'tecploteta(10,27,2)=',eeta(1,27,32), eeta(1,27,60), eeta(1,27,90)
      CALL Read_array ('rcmeflux' , record_number, label, ARRAY_3D = eflux )
      CALL Read_array ('rcmeavg'  , record_number, label, ARRAY_3D = eavg  )
!
      itime = label%intg(6)
      time_min = label%intg(3)*60 + label%intg(4) + NINT(REAL(label%intg(5))/60)
      WRITE (time_string,'(I2.2,A1,I2.2,A1,I2.2)') label%intg(3),":",label%intg(4),":",label%intg(5)
      WRITE (*,'(A,I6.6,A,A10,A)',ADVANCE='NO') 'ITIME=', itime, ' T=', time_string
!
!
      WRITE (filename_out_2d(15:16),'(I2.2)') label%intg(3)
      WRITE (filename_out_2d(17:18),'(I2.2)') label%intg(4)
      WRITE (filename_out_2d(19:20),'(I2.2)') label%intg(5)

      WRITE (filename_out_3d(15:16),'(I2.2)') label%intg(3)
      WRITE (filename_out_3d(17:18),'(I2.2)') label%intg(4)
      WRITE (filename_out_3d(19:20),'(I2.2)') label%intg(5)
!
      RCM_p = 0.
      RCM_p_40kv = 0.0
      do i=1,isize
       do j=1,jsize
         if( i<imin_j(j) ) then
            RCM_p(i,j) = -1.
            RCM_p_40kV (i,j) = -1.0
         else    
            do k=1,kcsize
	      
               RCM_p(i,j) = RCM_p(i,j) + vm(i,j)**2.5*eeta(i,j,k)*ABS(alamc(k))
               
               IF (alamc(k)*vm(i,j) >= 20.E+3 .AND. alamc(k)*vm(i,j) <= 60.E+3) THEN
                 RCM_p_40kv(i,j) = RCM_p_40kv(i,j) + &
                                   vm(i,j)**2.5*eeta(i,j,k)*ABS(alamc(k))
               END IF
	      
          end do  
         end if 
!	 end do
	 if (i==10) then 
!	 print*, 'eta=',j,eeta(10,j,31:40)
	 end if
      end do
!      do j=1,jsize
!      print*,'eta=',j, eeta(i,:,2) 
      end do
      
      RCM_p = RCM_p * 1.67E-35/ 1.0E-9
      RCM_p_40kv = RCM_p_40kv * 1.67E-35/ 1.0E-9
!      print*,'p=',RCM_p_40kv(10,27)

      RCM_pvgamma = 0.
      do i=1,isize; do j=1,jsize
         if( i<imin_j(j)-1 ) then
            RCM_pvgamma(i,j) = -1.
         else    
            do k=1,kcsize
!            do k=31,115
               RCM_pvgamma(i,j) = RCM_pvgamma(i,j) + eeta(i,j,k)*ABS(alamc(k))
            end do  
!	print*, 'pvgamma=',RCM_pvgamma(:,27)
         end if  
      end do; end do
      RCM_pvgamma = RCM_pvgamma * 2.0/3.0
      RCM_pvgamma = RCM_pvgamma /3.98E+25 ! convert to [nPa*(Re/nT)**(5/3)]
!
      RCM_n = 0.
      RCM_T = 0.
      do i=1,isize
       do j=1,jsize
             if( i<imin_j(j)-1 ) then
            RCM_T(i,j) = -1.
            RCM_n(i,j) = -1.
         else if (rcm_pvgamma(i,j) == 0.0) THEN
!	 if (rcm_pvgamma(i,j) == 0.0) THEN
	 
            RCM_n(i,j) = 0.0
            RCM_t(i,j) = 0.0
         else    
            do k=1,kcsize
               IF (alamc(k) <= 0.0) CYCLE
               RCM_n(i,j) = RCM_n(i,j) + eeta(i,j,k)*vm(i,j)**1.5/6.37E+21
            end do  
            RCM_T(i,j) = vm(i,j)*RCM_pvgamma(i,j)/SUM(eeta(i,j,:))
!            RCM_T(i,j) = vm(i,j)*RCM_pvgamma(i,j)/SUM(eeta(i,j,31:115))
         end if  
         MHD_P(i,j) = density(i,j)*1.0E+6*temperature(i,j)*1.6E-19 / 1.0E-9 ![nPa]
         IF (MHD_p(i,j) < 0.0) MHD_p(i,j) = 0.0
      end do
      end do
!      print*,'J=27', RCM_n(30:31, 27)
!      print*,'J=26',RCM_n(20:30,26)
!      print*,'density=',RCM_n(:,27)
!      print*,'temp=',RCM_T(:,27)
!
!!    CALL Wrap_around_ghostcells (RCM_n, isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (RCM_T, isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (RCM_p, isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (birk, isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (birk_mhd, isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (RCM_pvgamma,isize,jsize,n_gc)
!!    CALL Wrap_around_ghostcells (MHD_p,isize,jsize,n_gc)
!

      CALL Date_and_time (char_date, char_time)


      WRITE (*,'(A24,A2,A23)',ADVANCE='NO') filename_out_2d
      OPEN (UNIT=LUN, FILE = filename_out_2d, STATUS='REPLACE')

         WRITE (LUN,*)'TITLE="RCM-2D, RUN='//run_id_string//': T='//time_string//'"'
         WRITE (LUN,'(A)')'VARIABLES="I","J","COLAT","ALOCT","MLT",&
                     &"BNDLOC","XMIN","YMIN","VM", "BMIN",&
                     &"V", "BIRK(NH)", "PEDLAM", "PEDPSI","HALL","EFLUX",&
                     &"EAVG","N(RCM),cm-3","T(RCM),eV","P(RCM),nPa",&
                     &"PV_gamma","Birk_mhd(NH)","N(mhd),cm-3",&
                     &"T(MHD),eV","P(MHD),nPa","VPAR,V","VCOROT,V","ETA_PS","P_40kV,nPa","SS_eq"'
!
         WRITE (LUN,*) 'ZONE T="RCM-2D-'//time_string//'" I=', isize, ', J=',jsize+0, ',&
                        &F=POINT AUXDATA T="'//time_string//'"'
         WRITE (LUN,'(A,F5.1,A)') 'AUXDATA PCP='//'"',label%real(16),'"'
         WRITE (LUN,'(A,A,A)')    'AUXDATA RUN='//'"',run_id_string,'"'
         WRITE (LUN,'(A,F5.1,A)') 'AUXDATA IMF_Bz='//'"',label%real(19),'"'
         WRITE (LUN,'(A,F5.1,A)') 'AUXDATA Kp='//'"',label%real(17),'"'
         WRITE (LUN,'(A,F6.1,A)') 'AUXDATA Dst='//'"',label%real(14),'"'
         WRITE (LUN,'(A,F5.1,A)') 'AUXDATA P_dyn='//'"',label%real(11),'"'
         WRITE (LUN,'(A,A4,A1,A2,A1,A2,A1,A2,A1,A2,A1,A2,A1)')    &
               'AUXDATA Time_created='//'"', &
               char_date(1:4), '/', char_date(5:6), '/', CHAR_DATE(7:8), ' ', &
               char_time(1:2), ':', char_time(3:4), ':', char_time(5:6), '"'

         DO j = 1, jsize+0
         DO i = 1, isize
             WRITE (LUN,*) i, j, colat(i,j)*rtd, aloct(i,j), MODULO(aloct(i,j)*rth+12.0,24.01),&
                           bndloc(j), xmin(i,j), ymin(i,j), vm(i,j), bmin(i,j), &
                           v(i,j), 0.5*birk(i,j), pedlam(i,j), pedpsi(i,j),&
                           hall(i,j), eflux(i,j,1), eavg(i,j,1), &
                           RCM_n(i,j), RCM_T(i,j), RCM_P(i,j), RCM_pvgamma(i,j),&
                           -birk_mhd(i,j)/1.0E-6, density(i,j), temperature(i,j),&
                           MHD_p(i,j),vpar(i,j),vcorot(i,j) , eeta(i,j,1), RCM_p_40kv(i,j),&
                           ss(j)
!         print*, RCM_n(i,j)
!	 print*, temperature(i,j)			   
         END DO
         END DO
	 print*,'RCM_n=',rcm_n(30:31,27),rcm_T(30:31,27)*6.4*1e6*6.24*1.e15,rcm_P(30:31,27)
!
      CLOSE (LUN)
      WRITE (*,'(A)') '... DONE' 
!
!
IF (do_3d) THEN

      WRITE (*,'(A,I6.6,A,A10,A)',ADVANCE='NO') 'ITIME=', itime, ' T=', time_string
      WRITE (*,'(A24,A2,A23)',ADVANCE='NO') filename_out_3d
      OPEN (UNIT=LUN, FILE = filename_out_3d, STATUS='REPLACE')
         WRITE (LUN,*) 'TITLE="RCM-ENERGY-FILE: T='//time_string//'"'

         IF (K_selected  < 0) THEN
            k_start = 1
            k_stop  = ksize
         ELSE
            k_start = k_selected
            k_stop  = k_selected
         END IF

         WRITE (LUN,'(A)') 'VARIABLES="I","J","COLAT","ALOCT","BNDLOC","XMIN","YMIN","VM", "BMIN",&
                              &"V", "BIRK(N.H)", "EETA", "VEFF","W,eV"'
!
         DO k_zone = k_start, k_stop

            WRITE (char_kvalue,'(A2,I3.3)') 'K=', k_zone

            IF (k_zone == k_start) THEN
               WRITE (LUN,*) 'ZONE T="RCM-'//char_kvalue//'-'//time_string//&
                                 '" I=', isize, ', J=',jsize+0, ', F=POINT'
            ELSE
               WRITE (LUN,*) 'ZONE T="RCM-'//char_kvalue//'-'//time_string//&
                                 &'" I=', isize, ', J=',jsize+0, ', F=POINT, VARSHARELIST=([1-11]=1)'
            END IF

            WRITE (LUN,'(A,F5.1,A)') 'AUXDATA PCP='//'"',label%real(16),'"'
            WRITE (LUN,'(A,A,A)')    'AUXDATA RUN='//'"',run_id_string,'"'
            WRITE (LUN,'(A,F5.1,A)') 'AUXDATA IMF_Bz='//'"',label%real(19),'"'
            WRITE (LUN,'(A,F5.1,A)') 'AUXDATA Kp='//'"',label%real(17),'"'
            WRITE (LUN,'(A,F6.1,A)') 'AUXDATA Dst='//'"',label%real(14),'"'
            WRITE (LUN,'(A,F5.1,A)') 'AUXDATA P_dyn='//'"',label%real(11),'"'
            WRITE (LUN,'(A,I5.5,A)')   'AUXDATA Kvalue='//'"', k_zone,'"'
            WRITE (LUN,'(A,F12.3,A)') 'AUXDATA Alam='//'"', alamc(k_zone),'"'
            WRITE (LUN,'(A,I2.2,A)') 'AUXDATA IKFLAV='//'"', ikflavc(k_zone),'"'
            WRITE (LUN,'(A,ES12.3,A)') 'AUXDATA Xmass='//'"', xmass(ikflav(k_zone)),'"'
            WRITE (LUN,'(A,A4,A1,A2,A1,A2,A1,A2,A1,A2,A1,A2,A1)')    &
               'AUXDATA Time_created='//'"', &
               char_date(1:4), '/', char_date(5:6), '/', CHAR_DATE(7:8), ' ', &
               char_time(1:2), ':', char_time(3:4), ':', char_time(5:6), '"'

            veff = v + vcorot + vm*alamc(k_zone)
!!          CALL Wrap_around_ghostcells (veff, isize,jsize,n_gc)

            IF (k_zone == K_start) THEN
               DO j = 1, jsize+0
               DO i = 1, isize
                WRITE (LUN,*) i, j, colat(i,j)*rtd, aloct(i,j),&
                           bndloc(j), xmin(i,j), ymin(i,j), vm(i,j), bmin(i,j), &
                           v(i,j), birk(i,j)/2.0, eeta(i,j,k_zone), veff(i,j), alamc(k_zone)*vm(i,j)
               END DO
               END DO
            ELSE
               DO j = 1, jsize+0
               DO i = 1, isize
                WRITE (LUN,*) eeta(i,j,k_zone), veff(i,j), alamc(k_zone)*vm(i,j)
               END DO
               END DO
            END IF

         END DO

      CLOSE (LUN)
      WRITE (*,'(A)') '... DONE' 
END IF

   END DO


      STOP
      END PROGRAM Convert_rcm_to_tecplot