source: trunk/INPUT/Greeneemxy/changement_grille.f90 @ 308

program changement_grille
  implicit none

  integer i,j,ii,jj                 ! entiers de boucle
  integer aux, numx, numy           ! un entier auxiliaire

  integer,parameter :: Xmin = -990  ! les bornes du domaine (en km)
  integer,parameter :: Xmax = 895   ! dans la grille de Bamber
  integer,parameter :: Ymin = -3515
  integer,parameter :: Ymax = -550

  integer,parameter :: pasx_init = 5      ! pas initial selon x : 5km
  integer,parameter :: pasy_init = 5      ! pas initial selon y : 5km

!  integer,parameter :: pasx_new = 15      ! nouveau pas selon x : 15km
!  integer,parameter :: pasy_new = 15      ! nouveau pas selon y : 15km

  integer,parameter :: pasx_new = 45      ! nouveau pas selon x : 45km
  integer,parameter :: pasy_new = 45      ! nouveau pas selon y : 45km

! nx_init=((Xmax-Xmin)+pasx_init)/pasx_init
! nx_new=((Xmax-Xmin)+pasx_init)/pasx_new
  integer,parameter :: nx_init = 378    !  nb initial d'elements en x
  integer,parameter :: ny_init = 594    !  nb initial d'elements en y 
! nb d'elmt en x et y, pour 15km :
!  integer,parameter :: nx_new = 126     ! nouveau nb elements en x
!  integer,parameter :: ny_new = 198     ! nouveau nb elements en y
! nb d'elmt en x et y, pour 45km : 
  integer,parameter :: nx_new = 42
  integer,parameter :: ny_new = 66


  integer indicex  ! un indice pour la matrice initiale, x
  integer indicey  ! un indice pour la matrice initiale, y

  integer pt_courantx
  integer pt_couranty

  character*99 out_surface
  character*99 out_socle
  character*99 out_thick
  character*99 out_grid

  real,dimension(nx_init,ny_init) :: Surf_init  ! Matrice Surface initiale
  real,dimension(nx_new,ny_new) :: Surf_new     ! Nouvelle matrice Surface

  real,dimension(nx_init,ny_init) :: Socle_init ! Matrice Socle initiale
  real,dimension(nx_new,ny_new) :: Socle_new    ! Nouvelle matrice Socle

  real,dimension(nx_init,ny_init) :: Thick_init ! Matrice Epaisseur initiale
  real,dimension(nx_new,ny_new) :: Thick_new    ! Nouvelle matrice Epaisseur

  character(len=60) :: dir_data
  dir_data='/pc96/work250g/aquiquet/DATA/BASE_PROPRE/IN_OUTPUT/TEMP'

! Pour passer à 15 km :
!  out_surface=trim(dir_data)//'/surface_xyz_15km.txt'
!  out_socle=trim(dir_data)//'/socle_xyz_15km.txt'
!  out_thick=trim(dir_data)//'/thick_xyz_15km.txt'
!  out_grid=trim(dir_data)//'/grid_xy_15km.txt'

! Pour passer à 45 km :
  out_surface=trim(dir_data)//'/surface_xyz_45km.txt'
  out_socle=trim(dir_data)//'/socle_xyz_45km.txt'
  out_thick=trim(dir_data)//'/thick_xyz_45km.txt'
  out_grid=trim(dir_data)//'/grid_xy_45km.txt'

!__________________________________________________________________________
! Lecture des donnees au pas initial
  open(71, file=trim(dir_data)//'/surface_xyz_5km.txt',   &
                 form='formatted',status='old')
  open(72, file=trim(dir_data)//'/socle_xyz_5km.txt',     &
                 form='formatted',status='old')
  open(73, file=trim(dir_data)//'/thick_xyz_5km.txt',     &
                 form='formatted',status='old')

  do j=1,ny_init          ! Attention : j=1 au nord, j=ny_init au sud
     do i=1,nx_init       ! i=1 à l'ouest 
        read(71,*) aux, aux, Surf_init(i,j)
        read(72,*) aux, aux, Socle_init(i,j)
        read(73,*) aux, aux, Thick_init(i,j)
     end do
  end do
  close(71)
  close(72)
  close(73)
! Fin de lecture
!--------------------------------------------------------------------------




!__________________________________________________________________________
! Calcul des donnees au pas plus large
  numx = pasx_new/pasx_init
  numy = pasy_new/pasy_init

  do j=1,ny_new
     do i=1,nx_new
        Surf_new(i,j)=0
        Socle_new(i,j)=0
        Thick_new(i,j)=0
        indicex=((numx-1)/2 + 1) + (i-1)*numx
        indicey=((numy-1)/2 + 1) + (j-1)*numy
        do jj=1,numy        ! on va sommer tous les elmt de la petite grille
           do ii=1,numx
              Surf_new(i,j)=Surf_new(i,j) +             &
                   Surf_init(indicex+(ii-1)-(numx-1)/2,   &
                                indicey+(jj-1)-(numy-1)/2)
              Socle_new(i,j)=Socle_new(i,j) +                 &
                   Socle_init(indicex+(ii-1)-(numx-1)/2,     &
                              indicey+(jj-1)-(numy-1)/2)
              Thick_new(i,j)=Thick_new(i,j) +                 &
                   Thick_init(indicex+(ii-1)-(numx-1)/2,     &
                              indicey+(jj-1)-(numy-1)/2)
              if(i.eq.50) then           ! un test
                 print*, 'i,j,ii,jj,socle(i,j)', i,j,ii,jj, Socle_new(i,j)
              end if
           end do
        end do
        Surf_new(i,j)=Surf_new(i,j)/(numx*numy)    ! on veut la moyenne
        Socle_new(i,j)=Socle_new(i,j)/(numx*numy)
        Thick_new(i,j)=Thick_new(i,j)/(numx*numy)
     end do
  end do

! Fin de calcul
!--------------------------------------------------------------------------

!__________________________________________________________________________
! Ecriture des donnees au pas de temps plus large
  open(81, file=out_surface,form='formatted',status='unknown')
  open(82, file=out_socle,form='formatted',status='unknown')
  open(83, file=out_thick,form='formatted',status='unknown')
  pt_couranty=Ymax-(numy-1)*pasy_init/2    ! /I\ on avait j=1 au nord.
  do j=1,ny_new
     pt_courantx=Xmin+(numx-1)*pasx_init/2
     do i=1,nx_new
        write(81,*) pt_courantx, pt_couranty, Surf_new(i,j)
        write(82,*) pt_courantx, pt_couranty, Socle_new(i,j)
        write(83,*) pt_courantx, pt_couranty, Thick_new(i,j)
        pt_courantx = pt_courantx +  pasx_new      
     enddo
     pt_couranty = pt_couranty -  pasy_new
  enddo
  close(81)
  close(82)
  close(83)
! Fin d'ecriture
!--------------------------------------------------------------------------

!--------------------------------------------------------------------------
! Ecriture d'un simple fichier grille x y :
  open(91, file=out_grid,form='formatted',status='unknown')
  pt_couranty=Ymax-(numy-1)*pasy_init/2    ! /I\ on avait j=1 au nord.
  do j=1,ny_new
     pt_courantx=Xmin+(numx-1)*pasx_init/2
     do i=1,nx_new
        write(91,*) pt_courantx, pt_couranty        
        pt_courantx = pt_courantx +  pasx_new      
     end do     
     pt_couranty = pt_couranty -  pasy_new
  end do
  close(91)

end program changement_grille