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

program formatage_grisli

  implicit none

  character(len=60) :: dir_data    ! le repertoire ou se trouvent les donnees
  character*99 in_surface          ! surface
  character*99 in_socle            ! socle
  character*99 in_thick            ! epaisseur
  character*99 out                 ! fichier de sortie
  character*99 out2                ! fichier de sortie pour grille x y
 ! character*99 out3                ! fichier de sortie pour grille lon lat
  character*99 in_lonlat           ! fichier d entree lon/lat pour la grille

  character*2 resolution           ! la resolution de la maille (string), a changer
  integer pas                      ! la resolution en entier (pas besoin de le changer

! limites "vraies" du domaine, mais pas position du premier "cell center"
  real,parameter :: limx_g = -992.5    ! x a gauche
  real,parameter :: limx_d = 897.5     ! x a droite
  real,parameter :: limy_b = -3517.5   ! y en bas
  real,parameter :: limy_h = -547.5    ! y en haut

  integer nx           ! nb d'elements en x
  integer ny           ! nb d'elements en y
  integer X            ! position en x
  integer Y            ! position en y

  integer i,j         ! entiers de boucle

  real,dimension(800,800) :: S     ! Matrice Surface 
  real,dimension(800,800) :: B     ! Matrice Socle 
  real,dimension(800,800) :: H     ! Matrice Epaisseur 

  real,dimension(800,800) :: lon   ! Matrice des longitudes 
  real,dimension(800,800) :: lat   ! Matrice des latitudes

  real xmin, ymin          ! premier centre
  real xmax, ymax          ! dernier centre

700 format (3(f12.4))
!800 format (4(f10.2),i5)
!801 format (4(f10.2), 2(i5))

  dir_data='/pc96/work250g/aquiquet/DATA/BASE_PROPRE/IN_OUTPUT/TEMP'

! Voici l'unique parametre a changer : la resolution spatiale en km ->
  resolution = '45'                
  read(resolution,*) pas  ! on convertit ce string en entier

  in_surface=trim(dir_data)//'/surface_xyz_'//trim(resolution)//'km.txt'
  in_socle=trim(dir_data)//'/socle_xyz_'//trim(resolution)//'km.txt'
  in_thick=trim(dir_data)//'/thick_xyz_'//trim(resolution)//'km.txt'
  out=trim(dir_data)//'/groen_grisli_'//trim(resolution)//'km.txt'
  out2=trim(dir_data)//'/groen_grisli_grid_xylonlat_'//trim(resolution)//'km.txt'
!  out3=trim(dir_data)//'/groen_grisli_grid_lonlat_'//trim(resolution)//'km.txt'
  in_lonlat=trim(dir_data)//'/grid_'//trim(resolution)//'km_lonlat.txt'

  print*, pas
  print*, in_thick

  select case (pas)
  case(5)
     nx=378
     ny=594
  case(15)
     nx=126
     ny=198
  case(45)
     nx=42
     ny=66
  case default
     print*, 'erreur, resolution non valide'
     nx=1
     ny=1
  end select

  print*, nx,ny

! on definit la position des limites du domaine (grid registration)
  xmin=limx_g+pas/2
  ymin=limy_b+pas/2
  xmax=limx_d-pas/2
  ymax=limy_h-pas/2
  
  open(71, file=in_surface,form='formatted',status='old')
  open(72, file=in_socle,form='formatted',status='old')
  open(73, file=in_thick,form='formatted',status='old')

  do j=1,ny    ! attention du nord au sud dans fichier initial
     do i=1,nx
        read(71,*) X, Y, S(i,ny-j+1)
        read(72,*) X, Y, B(i,ny-j+1)
        read(73,*) X, Y, H(i,ny-j+1)
     end do
  end do

  close(71)
  close(72)
  close(73)

  open(74, file=out,form='formatted', status='unknown')
  write(74,*) 'Groenland : Topo actuelle (ETOPO + Bamber)'
  write(74,*) 'nx, ny et resolution :'
  write(74,*) nx, ny, pas
  write(74,*) '     S                H               B'
  do j=1,ny   
     do i=1,nx  
        write(74,700) S(i,j),H(i,j), B(i,j)
     end do
  end do
  close(74)

  open(76, file=in_lonlat,form='formatted', status='old')     
  do j=1,ny            ! attention j=1 au nord 
     do i=1,nx
        read(76,*) lon(i,ny-j+1), lat(i,ny-j+1)
     end do
  end do   
  close(76)

  open(75, file=out2,form='formatted', status='unknown')     
  write(75,*) 'Groenland : Grille x y'
  write(75,*) 'nx, ny et resolution :'
  write(75,*) nx, ny, pas
  write(75,*) '          I           J           X           Y        LON            LAT'
  Y=ymin
  do j=1,ny 
     X=xmin
     do i=1,nx  
        write(75,*) i,j, X,Y,lon(i,ny-j+1), lat(i,ny-j+1)
        X=X+pas
     end do
     Y=Y+pas
  end do
  close(75)



!!$  open(77, file=out3,form='formatted', status='unknown')     
!!$  write(77,*) 'Groenland : Grille lon lat'
!!$  write(77,*) 'nx, ny et resolution :'
!!$  write(77,*) nx, ny, pas
!!$  write(77,*) '    I      J        X       Y'
!!$  do j=1,ny 
!!$     do i=1,nx  
!!$        write(77,*) i,j, lon(i,j), lat(i,j)
!!$     end do
!!$  end do
!!$  close(77)


end program formatage_grisli