[4] | 1 | program formatage_grisli |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | character(len=60) :: dir_data ! le repertoire ou se trouvent les donnees |
---|
| 6 | character*99 in_surface ! surface |
---|
| 7 | character*99 in_socle ! socle |
---|
| 8 | character*99 in_thick ! epaisseur |
---|
| 9 | character*99 out ! fichier de sortie |
---|
| 10 | character*99 out2 ! fichier de sortie pour grille x y |
---|
| 11 | ! character*99 out3 ! fichier de sortie pour grille lon lat |
---|
| 12 | character*99 in_lonlat ! fichier d entree lon/lat pour la grille |
---|
| 13 | |
---|
| 14 | character*2 resolution ! la resolution de la maille (string), a changer |
---|
| 15 | integer pas ! la resolution en entier (pas besoin de le changer |
---|
| 16 | |
---|
| 17 | ! limites "vraies" du domaine, mais pas position du premier "cell center" |
---|
| 18 | real,parameter :: limx_g = -992.5 ! x a gauche |
---|
| 19 | real,parameter :: limx_d = 897.5 ! x a droite |
---|
| 20 | real,parameter :: limy_b = -3517.5 ! y en bas |
---|
| 21 | real,parameter :: limy_h = -547.5 ! y en haut |
---|
| 22 | |
---|
| 23 | integer nx ! nb d'elements en x |
---|
| 24 | integer ny ! nb d'elements en y |
---|
| 25 | integer X ! position en x |
---|
| 26 | integer Y ! position en y |
---|
| 27 | |
---|
| 28 | integer i,j ! entiers de boucle |
---|
| 29 | |
---|
| 30 | real,dimension(800,800) :: S ! Matrice Surface |
---|
| 31 | real,dimension(800,800) :: B ! Matrice Socle |
---|
| 32 | real,dimension(800,800) :: H ! Matrice Epaisseur |
---|
| 33 | |
---|
| 34 | real,dimension(800,800) :: lon ! Matrice des longitudes |
---|
| 35 | real,dimension(800,800) :: lat ! Matrice des latitudes |
---|
| 36 | |
---|
| 37 | real xmin, ymin ! premier centre |
---|
| 38 | real xmax, ymax ! dernier centre |
---|
| 39 | |
---|
| 40 | 700 format (3(f12.4)) |
---|
| 41 | !800 format (4(f10.2),i5) |
---|
| 42 | !801 format (4(f10.2), 2(i5)) |
---|
| 43 | |
---|
| 44 | dir_data='/pc96/work250g/aquiquet/DATA/BASE_PROPRE/IN_OUTPUT/TEMP' |
---|
| 45 | |
---|
| 46 | ! Voici l'unique parametre a changer : la resolution spatiale en km -> |
---|
| 47 | resolution = '45' |
---|
| 48 | read(resolution,*) pas ! on convertit ce string en entier |
---|
| 49 | |
---|
| 50 | in_surface=trim(dir_data)//'/surface_xyz_'//trim(resolution)//'km.txt' |
---|
| 51 | in_socle=trim(dir_data)//'/socle_xyz_'//trim(resolution)//'km.txt' |
---|
| 52 | in_thick=trim(dir_data)//'/thick_xyz_'//trim(resolution)//'km.txt' |
---|
| 53 | out=trim(dir_data)//'/groen_grisli_'//trim(resolution)//'km.txt' |
---|
| 54 | out2=trim(dir_data)//'/groen_grisli_grid_xylonlat_'//trim(resolution)//'km.txt' |
---|
| 55 | ! out3=trim(dir_data)//'/groen_grisli_grid_lonlat_'//trim(resolution)//'km.txt' |
---|
| 56 | in_lonlat=trim(dir_data)//'/grid_'//trim(resolution)//'km_lonlat.txt' |
---|
| 57 | |
---|
| 58 | print*, pas |
---|
| 59 | print*, in_thick |
---|
| 60 | |
---|
| 61 | select case (pas) |
---|
| 62 | case(5) |
---|
| 63 | nx=378 |
---|
| 64 | ny=594 |
---|
| 65 | case(15) |
---|
| 66 | nx=126 |
---|
| 67 | ny=198 |
---|
| 68 | case(45) |
---|
| 69 | nx=42 |
---|
| 70 | ny=66 |
---|
| 71 | case default |
---|
| 72 | print*, 'erreur, resolution non valide' |
---|
| 73 | nx=1 |
---|
| 74 | ny=1 |
---|
| 75 | end select |
---|
| 76 | |
---|
| 77 | print*, nx,ny |
---|
| 78 | |
---|
| 79 | ! on definit la position des limites du domaine (grid registration) |
---|
| 80 | xmin=limx_g+pas/2 |
---|
| 81 | ymin=limy_b+pas/2 |
---|
| 82 | xmax=limx_d-pas/2 |
---|
| 83 | ymax=limy_h-pas/2 |
---|
| 84 | |
---|
| 85 | open(71, file=in_surface,form='formatted',status='old') |
---|
| 86 | open(72, file=in_socle,form='formatted',status='old') |
---|
| 87 | open(73, file=in_thick,form='formatted',status='old') |
---|
| 88 | |
---|
| 89 | do j=1,ny ! attention du nord au sud dans fichier initial |
---|
| 90 | do i=1,nx |
---|
| 91 | read(71,*) X, Y, S(i,ny-j+1) |
---|
| 92 | read(72,*) X, Y, B(i,ny-j+1) |
---|
| 93 | read(73,*) X, Y, H(i,ny-j+1) |
---|
| 94 | end do |
---|
| 95 | end do |
---|
| 96 | |
---|
| 97 | close(71) |
---|
| 98 | close(72) |
---|
| 99 | close(73) |
---|
| 100 | |
---|
| 101 | open(74, file=out,form='formatted', status='unknown') |
---|
| 102 | write(74,*) 'Groenland : Topo actuelle (ETOPO + Bamber)' |
---|
| 103 | write(74,*) 'nx, ny et resolution :' |
---|
| 104 | write(74,*) nx, ny, pas |
---|
| 105 | write(74,*) ' S H B' |
---|
| 106 | do j=1,ny |
---|
| 107 | do i=1,nx |
---|
| 108 | write(74,700) S(i,j),H(i,j), B(i,j) |
---|
| 109 | end do |
---|
| 110 | end do |
---|
| 111 | close(74) |
---|
| 112 | |
---|
| 113 | open(76, file=in_lonlat,form='formatted', status='old') |
---|
| 114 | do j=1,ny ! attention j=1 au nord |
---|
| 115 | do i=1,nx |
---|
| 116 | read(76,*) lon(i,ny-j+1), lat(i,ny-j+1) |
---|
| 117 | end do |
---|
| 118 | end do |
---|
| 119 | close(76) |
---|
| 120 | |
---|
| 121 | open(75, file=out2,form='formatted', status='unknown') |
---|
| 122 | write(75,*) 'Groenland : Grille x y' |
---|
| 123 | write(75,*) 'nx, ny et resolution :' |
---|
| 124 | write(75,*) nx, ny, pas |
---|
| 125 | write(75,*) ' I J X Y LON LAT' |
---|
| 126 | Y=ymin |
---|
| 127 | do j=1,ny |
---|
| 128 | X=xmin |
---|
| 129 | do i=1,nx |
---|
| 130 | write(75,*) i,j, X,Y,lon(i,ny-j+1), lat(i,ny-j+1) |
---|
| 131 | X=X+pas |
---|
| 132 | end do |
---|
| 133 | Y=Y+pas |
---|
| 134 | end do |
---|
| 135 | close(75) |
---|
| 136 | |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | !!$ open(77, file=out3,form='formatted', status='unknown') |
---|
| 140 | !!$ write(77,*) 'Groenland : Grille lon lat' |
---|
| 141 | !!$ write(77,*) 'nx, ny et resolution :' |
---|
| 142 | !!$ write(77,*) nx, ny, pas |
---|
| 143 | !!$ write(77,*) ' I J X Y' |
---|
| 144 | !!$ do j=1,ny |
---|
| 145 | !!$ do i=1,nx |
---|
| 146 | !!$ write(77,*) i,j, lon(i,j), lat(i,j) |
---|
| 147 | !!$ end do |
---|
| 148 | !!$ end do |
---|
| 149 | !!$ close(77) |
---|
| 150 | |
---|
| 151 | |
---|
| 152 | end program formatage_grisli |
---|