  program crea_grid_B
  
! -----------------------------------------------------------------------------
! Frame   : Procedures and program used in 
!
! Program : crea_grid_B.f90
! Object  : Compute B 3D grid from FGMPOS_database
! Input   : file containig all FGM_POS.dat files
! Output  : B 3D grid file
! Author  : P. Robert, LPP-ScientiDev, Sept 2025
! Mail    : patrick.robert@lpp.polytechnique.fr
! -----------------------------------------------------------------------------


  character(len=255) :: file_fgmpos,datafileB,name_dir
  character(len=24)  :: datiso
  character(len=2)   :: CLU
  character(len=3)   :: cdta,Kp_on,sIGRF,bid
  real(kind=4), dimension(:,:,:,:,:,:), allocatable :: gridB
  integer,      dimension(:,:,:,:,:),   allocatable :: nvalB
  real(kind=4), dimension(9) :: sw
  
  character(len=13) :: datim(80000),curdate
  integer           :: ikp(80000)
  
  print*,'Kp_on ?(Yes or No) if yes, one grid by Kp, need memory'
  read(*,'(a)') Kp_on
  print*, Kp_on
  
  print*, 'IGRF substraction ? (Yes or No)'
  read(*,'(a)') sIGRF
  print*, sIGRF

  if(Kp_on == 'Yes') then
     a=18.
     idimkp=5
     dd=0.30
                     else
     a=20.
     idimkp=1
     dd=0.25
  endif
  
  x1=-a
  x2= a
  y1=-a
  y2= a
  z1=-a
  z2= a
  
  dx= dd
  dy= dd
  dz= dd

  nx= int((x2-x1)/dx +0.5)
  ny= int((y2-y1)/dy +0.5)
  nz= int((z2-z1)/dz +0.5)

  print*, '-----------------------------------------------------------------'
  print*, 'nx,ny,nz=',nx,ny,nz

  deltax=float(nx)*dx
  deltay=float(ny)*dy
  deltaz=float(nz)*dz

  difx=(x2-x1)-deltax
  if(abs(difx) > 1.e-10) then
     print*, 'x2-x1 =',x2-x1
     print*, 'deltax=',deltax
     print*, 'diff  =',difx
     print*, '***incoherence! aborted!!'
     stop '***incoherence! aborted!!'
  endif

  dify=(y2-y1)-deltay
  if(abs(dify) > 1.e-10) then
     print*, 'y2-y1 =',y2-y1
     print*, 'deltay=',deltay
     print*, 'diff  =',dify
     print*, '***incoherence! aborted!!'
     stop '***incoherence! aborted!!'
  endif

  difz=(z2-z1)-deltaz
  if(abs(difz) > 1.e-10) then
     print*, 'z2-z1 =',z2-z1
     print*, 'deltaz=',deltaz
     print*, 'diff  =',difz
     print*, '***incoherence! aborted!!'
     stop '***incoherence! aborted!!'
  endif

  allocate(gridB(4,nx,ny,nz,9,idimkp))

  allocate(nvalB(nx,ny,nz,9,idimkp))

  gridB(:,:,:,:,:,:)=0.

  nvalB(:,:,:,:,:)=0
  
! read Kp file
  if(Kp_on == 'Yes') call read_Kp_data(datim,ikp,nrec)
  
! ---------------
! read files fgmpos*.dat
! ---------------

  iosta=0
  igx=0
  igy=0
  igz=0
  sw(:)=0.
  nf=0
  ntot=0
  
! reading all files
! =================  

  DO while (iosta == 0)
     print*, '------------------------------------------------------------------------'
     print*, 'path of file to read ? (ex: FGMPOS.dat)'
     read (*,'(a)', iostat=iosta) file_fgmpos
     print*,'reading file: ',trim(file_fgmpos)
 
! probleme geopack sur un fichier 
     if(sIGRF == 'Yes' .and. file_fgmpos == 'FGMPOS_asc_database/2004/2004_09/C1_FGM_SPIN_20040904.txt') then
       print*, '*** anomalie geopack_2008 on this file, ignored'
       cycle
     endif
     
     if(iosta /= 0) then
        print*,'*** file does not exist, ignored'
        cycle
     endif
     if(file_fgmpos(1:5) == '     ') then
        print*,'*** file name is blank, ignored'
        cycle
     endif
     print*, trim(file_fgmpos)
   
     close(1)
     open(1, file=file_fgmpos,iostat=iosta3)
     if(iosta3 /= 0) then
        print*,'*** file cannot be openned, ignored'
        cycle
                     else
        print*,'file opened'
     endif
     nf=nf+1
     nt=0
     nerrR=0
     nerrD=0
   
   ! reading current file
   ! --------------------
     iosta2=0
       
     DO while(iosta2 == 0)  ! reading each record of the file
   
         read(1,10,iostat=iosta2, end=20) datiso,bx,by,bz,px,py,pz,CLU
         if(iosta2 /= 0) then
                 print*,'iosta2=',iosta2
                 print*,'record #',nt,' cannot be read'
                 nerrR=nerrR+1
                 go to 20
         endif
      10 format(a24,1x,6E14.6,1x,a2)
       
       ! passage en gsm
       
         read(datiso,'(i4,5(1x,i2))') iyear,imon,iday,ih,im,is
       
         call ctimpar(iyear,imon,iday,ih,im,is)
         call t_gse_to_gsm(px,py,pz,pxg,pyg,pzg)
         call t_gse_to_gsm(bx,by,bz,bxg,byg,bzg)
         call g_gsm_dipole_tilt_angle(diptan)
         diptan=diptan*180./3.14159
       
       ! km -> RE
       
         RE=6372.2
       
         pxg=pxg/RE
         pyg=pyg/RE
         pzg=pzg/RE
         
         r=sqrt(pxg**2 +pyg**2 +pzg**2)
  
        if(r > 25. .or. r < 1.1) cycle    
         
        nt=nt+1
        ntot=ntot+1
        
! on retire l'IGRF si demande    

        IF(sIGRF == 'Yes') THEN
        
!   (3) IF NO   INFORMATION IS AVAILABLE ON THE SOLAR WIND SPEED, THEN SET VGSEX=-400.0
!       AND VGSEY=VGSEZ=0. IN THAT CASE,    THE GSW COORDINATE SYSTEM BECOMES
!       IDENTICAL TO    THE STANDARD ONE. 

!       call DIP_08 (pxg,pyg,pzg,HXGSW,HYGSW,HZGSW) 
  
        VGSEX=-400.0 
        VGSEY=0.
        VGSEZ=0.
        call cdoyear(iyear,imon,iday,idoty)
        call RECALC_08 (IYEAR,idoty,ih,im,is,VGSEX,VGSEY,VGSEZ)
        call GDIPTIL(Tdiptan)
        Tdiptan=Tdiptan*180./3.14159 
  
        if(abs(diptan-Tdiptan) > 1.) then
           nblorej=nblorej+1
           print*, '*** dipta, Tdiptan,diff=',diptan,Tdiptan,diptan-Tdiptan
           print*, '*** block rejected :',nblorej,' /',nbvec
           cycle
        endif
  
        call IGRF_GSW_08 (pxg,pyg,pzg,HXGSW,HYGSW,HZGSW)
        
        bxg=bxg -HXGSW
        byg=byg -HYGSW
        bzg=bzg -HZGSW
      
        ENDIF ! on IGRF request
        
       ! compute Kp
       
         if(Kp_on == 'Yes') then 
            write(curdate,'(i4,"-",i2.2,"-",i2.2,1x,i2.2)') iyear,imon,iday,ih

            do k=1,nrec
               if(curdate(1:10) == datim(k)(1:10)) then
                  read(datim(k),'(11x,i2)') ihk 
                  if(ih == ihk-1 .or. ih == ihk .or. ih == ihk+1) then
                     kp=ikp(k)
                     exit
                   endif
               endif
            enddo
         
            kp4=kp/4 +1
            if(Kp4 > 5) Kp4=5
                            else
             Kp4=1
         endif
         
       ! calcul des grilles
       ! ------------------
       
         if(diptan > 0) then
            kdta= int(( diptan+5.)/10.)
                        else
            kdta=-int((-diptan+5.)/10.)
         endif
         idta=kdta+5  ! varie entre 1 et 9 pour dta entre -35 et 35
         nbdta=9
       
         if(idta < 1 .or. idta > nbdta) then
              print*,'**** idta=',idta,'  anormal!!!'
              nerrD=nerrD+1
              cycle
         endif
       
         dtak= float(kdta)*10.
       
         igx=int((pxg-x1)/dx) +1
         igy=int((pyg-y1)/dy) +1
         igz=int((pzg-z1)/dz) +1
       
         if(nt < 6 .and. nf < 5) print*,'record #',nt, ' igx,igy,igz=',igx,igy,igz,'  idta,Kp,kp4=',idta,kp,kp4
       
         if(igx < 1 .or. igx > nx) cycle
         if(igy < 1 .or. igy > ny) cycle
         if(igz < 1 .or. igz > nz) cycle
       
       ! B ----------------------------
         nvalB(igx,igy,igz,idta,kp4)=nvalB(igx,igy,igz,idta,kp4) +1
       
         if(nt == 7) print*, 'Reading file...'
       
         gridB(1,igx,igy,igz,idta,kp4)=gridB(1,igx,igy,igz,idta,kp4) +bxg
         gridB(2,igx,igy,igz,idta,kp4)=gridB(2,igx,igy,igz,idta,kp4) +byg
         gridB(3,igx,igy,igz,idta,kp4)=gridB(3,igx,igy,igz,idta,kp4) +bzg
       
         B0=sqrt(bxg**2 +byg**2 +bzg**2)
       
         gridB(4,igx,igy,igz,idta,kp4)=gridB(4,igx,igy,igz,idta,kp4) +B0
   
     ENDDO ! end of reading file
   
    20 continue  ! pour EOF sur fichier en cours
                 ! -----------------------------
   
     if(nt == 0) then
        print*, trim(file_fgmpos),' nf, nt=',nf,nt, '*** empty file'
                 else
        print*, trim(file_fgmpos),' nf, nt=',nf,nt
     endif
   
 ENDDO ! sur lecture de tous les fichiers
       ! ================================

  print*, 'nb read files = ',nf
  print*, 'nb all records=',ntot
  print*, 'nb errors Read=',nerrR
  print*, 'nb errors Dip =',nerrD

! definition noms repertoire et fichiers
! --------------------------------------

if(Kp_on == 'Yes') then
   write(bid,'(a2,i1)') 'Kp',kp4
   name_dir= 'grid_B_Kp'
   nca=9
                   else
   write(bid,'(a3)') 'AKp'
   name_dir= 'grid_B'  
   nca=6
endif

   
if(sIGRF == 'Yes') then
   name_dir=name_dir(1:nca)//'_mIGRF'
endif


! creation des fichiers grilles pour chaque valeur du dta (et eventuellement du Kp)
! -------------------------------------------------------
  print*
  print*,'grid creation, wait...'
  
do idta=1,nbdta
do kp4=1,idimkp

  dta=float(idta-5)*10.
  kdta=INT(dta)
  write(cdta,'(i3.2)') abs(kdta)
  if(kdta >= 0) then
          cdta(1:1)='p'
                else
          cdta(1:1)='m'
  endif
  print*, '---------------------'
  print*, 'idta,cdta,kp4=',idta,cdta,kp4

   write(datafileB, '(a,3i3.3,5a)') trim(name_dir)//'/grid_Bxyz_',nx,ny,nz,'_',cdta,'_',bid,'.dat'

  close(11)
  open(11,file=datafileB)

  write(11,'(i8,1x,i8,1x,i8,1x,i1)') nx,ny,nz,4

  nb=0

! normalisation

  do i=1,nx
     do j=1,ny
        do k=1,nz
          if(nvalB(i,j,k,idta,kp4) /= 0) then
             Bxx=gridB(1,i,j,k,idta,kp4)/float(nvalB(i,j,k,idta,kp4))
             Byy=gridB(2,i,j,k,idta,kp4)/float(nvalB(i,j,k,idta,kp4))
             Bzz=gridB(3,i,j,k,idta,kp4)/float(nvalB(i,j,k,idta,kp4))
             B00=gridB(4,i,j,k,idta,kp4)/float(nvalB(i,j,k,idta,kp4))
             nb=nb+1
                              else
             Bxx=1.e30
             Byy=1.e30
             Bzz=1.e30
             B00=1.e30
          endif
          write(11,'(4E13.5,i8)') Bxx,Byy,Bzz,B00, nvalB(i,j,k,idta,kp4)
        enddo
     enddo
  enddo

  print*,trim(datafileB)

  print*,'grid dimension      b :',nx,ny,nz

  close(11)

enddo ! idta
enddo ! kp4

  print*, ('-',i=1,72)

  print*, "crea_grid_B.exe  : NORMAL TERMINATION"
     stop "crea_grid_B.exe  : NORMAL TERMINATION"

  end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
       SUBROUTINE GDIPTIL(PS)
!
!-----INPUT PARAMETERS:  NONE
!-----OUTPUT PARAMETERS: PS - GEODIPOLE TILT ANGLE IN RADIANS
!
!     ATTENTION:  SUBROUTINE  RECALC  MUST BE INVOKED BEFORE GDIPTIL IN
!     TWO CASES:
!     /A/  BEFORE THE FIRST INVOCATION OF GDIPTIL
!     /B/  IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
!
!     P. Robert, June 2002.
!
!                      ----------------------
!
      COMMON /GEOPACK1/ A(12),DS3,BB(2),PSI,CC(18)
!
      PS=PSI
!
      RETURN
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

