source: trunk/LMDZ.MARS/util/startarchive2icosa/rearrange_startphy.f90 @ 2745

Last change on this file since 2745 was 2532, checked in by adelavois, 4 years ago

update of start_archive2icosa tool
File size: 8.0 KB
Line 
1PROGRAM rearrange_startphy
2
3USE netcdf
4
5IMPLICIT NONE
6
7INTEGER :: ncid
8INTEGER :: dimids(4)
9INTEGER :: varid
10INTEGER :: ierr
11INTEGER :: nvar
12INTEGER :: ndim
13CHARACTER(LEN=100) :: varname
14INTEGER :: varname_dimids(5)
15
16INTEGER :: i,j,k,nb_cells,nlev,nsoil,nvertex,nlev_p1
17INTEGER,ALLOCATABLE :: cell_index(:)
18REAL,ALLOCATABLE :: ref_lon(:)
19REAL,ALLOCATABLE :: ref_lat(:)
20REAL,ALLOCATABLE :: lon(:)
21REAL,ALLOCATABLE :: lat(:)
22REAL,ALLOCATABLE :: ref_field(:),field(:)
23REAL,ALLOCATABLE :: ref_field_3D(:,:),field_3D(:,:)
24REAL,ALLOCATABLE :: ref_field_3D_p1(:,:),field_3D_p1(:,:)
25REAL,ALLOCATABLE :: ref_field_soil(:,:),field_soil(:,:)
26REAL,ALLOCATABLE :: ref_field_vertex(:,:),field_vertex(:,:)
27REAL :: lati,loni
28REAL :: diff_lat,diff_lon
29CHARACTER(LEN=*),PARAMETER :: ref_file="startphy_icosa_ref.nc"
30CHARACTER(LEN=*),PARAMETER :: file="startfi.nc"
31
32DOUBLE PRECISION,PARAMETER :: pi=acos(-1.d0)
33REAL :: reflatj,reflonj
34
35! load coordinates from files
36  ierr=NF90_OPEN(ref_file, NF90_NOWRITE, ncid)
37  ierr=NF90_INQ_DIMID(ncid,"physical_points",dimids(1))
38  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=nb_cells)
39  write(*,*) "nb_cells=",nb_cells,"dimids(1)=",dimids(1)
40  allocate(ref_lat(nb_cells),ref_lon(nb_cells))
41 
42  ierr=NF90_INQ_VARID(ncid,"latitude",varid)
43  ierr=NF90_GET_VAR(ncid,varid,ref_lat)
44  if (ierr /= nf90_noerr) then
45    write(*,*) "cannot load latitude from ",trim(ref_file)
46  else
47    write(*,*) "ref_lat(1:5)=",ref_lat(1:5)
48  endif
49  ierr=NF90_INQ_VARID(ncid,"longitude",varid)
50  ierr=NF90_GET_VAR(ncid,varid,ref_lon)
51  if (ierr /= nf90_noerr) then
52    write(*,*) "cannot load longitude from ",trim(ref_file)
53  else
54    write(*,*) "ref_lon(1:5)=",ref_lon(1:5)
55  endif
56  ierr=NF90_CLOSE(ncid)
57 
58  ierr=NF90_OPEN(file, NF90_WRITE, ncid)
59  allocate(lat(nb_cells),lon(nb_cells))
60  ierr=NF90_INQ_VARID(ncid,"latitude",varid)
61  ierr=NF90_GET_VAR(ncid,varid,lat)
62  if (ierr /= nf90_noerr) then
63    write(*,*) "cannot load lat from ",trim(file)
64  else
65    write(*,*) "lat(1:5)=",lat(1:5)
66  endif
67  ierr=NF90_INQ_VARID(ncid,"longitude",varid)
68  ierr=NF90_GET_VAR(ncid,varid,lon)
69  if (ierr /= nf90_noerr) then
70    write(*,*) "cannot load lat from ",trim(file)
71  else
72    write(*,*) "lon(1:5)=",lon(1:5)
73  endif
74 
75  ! find correspondances between lon/ref_lon & lat/ref_lat
76  allocate(cell_index(nb_cells))
77  cell_index(:)=0
78  do i=1,nb_cells !in case lon and lat are expressed in rad in startphy_icosa_ref /180*pi
79    lat(i)=lat(i) /(180.d0/pi) ; lon(i)=lon(i)/(180.d0/pi)
80    do j=1,nb_cells
81     reflatj=ref_lat(j)
82     reflonj=ref_lon(j)
83
84     if (lat(i) ==0.) then
85       diff_lat=abs((lat(i)-ref_lat(j))/1.e-8)
86     else
87       diff_lat=abs((lat(i)-ref_lat(j))/lat(i))
88     endif
89
90     if (lon(i) ==0.) then
91      diff_lon=abs((lon(i)-ref_lon(j))/1.e-8)
92     else
93      diff_lon=abs((lon(i)-ref_lon(j))/lon(i))
94     endif
95
96     if ((diff_lat <= 1.e-5).and.(diff_lon <= 1.e-5)) then
97      cell_index(i)=j
98      write(*,*)"j=",j," lat(i)=",lat(i)," ref_lat(j)=",reflatj
99      write(*,*)"j=",j," lon(i)=",lon(i)," ref_lon(j)=",reflonj
100     endif
101    enddo ! of do j=1,nb_cells
102    write(*,*) ">> i=",i,"cell_index(i)=",cell_index(i)
103    ! sanity check:
104    if (cell_index(i)==0) then
105      write(*,*) "Error, could not find lon-lat match for i=",i , lat(i), lon(i)
106      stop
107    endif
108    ! writing lat and lon in rad rather than deg
109    ierr=NF90_INQ_VARID(ncid,"longitude",varid)
110    ierr=NF90_PUT_VAR(ncid,varid,lon)
111
112    ierr=NF90_INQ_VARID(ncid,"latitude",varid)
113    ierr=NF90_PUT_VAR(ncid,varid,lat)
114  enddo ! of do i=1,nb_cells
115 
116!  do i=1,nb_cells
117!    write(*,*) "i=",i," cell_index(i)=",cell_index(i)
118!    write(*,*) "  lat(i)=",lat(i),"ref_lat(cell_index(i))=",ref_lat(cell_index(i))
119!  enddo
120 
121  ! load, rearrange and write variables
122  ierr=NF90_INQ_DIMID(ncid,"physical_points",dimids(1))
123  ierr=NF90_INQ_DIMID(ncid,"lev",dimids(2))
124  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=nlev)
125  ierr=NF90_INQ_DIMID(ncid,"subsurface_layers",dimids(3))
126  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(3), len=nsoil)
127  ierr=NF90_INQ_DIMID(ncid,"nvertex",dimids(4))
128  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(4), len=nvertex)
129  ierr=NF90_INQ_DIMID(ncid,"lev_p1",dimids(5))
130  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=nlev_p1)
131  ierr=NF90_INQUIRE(ncid,nVariables=nvar)
132  write(*,*) "nvar=",nvar
133  allocate(ref_field(nb_cells),field(nb_cells))
134  allocate(ref_field_3D(nb_cells,nlev),field_3D(nb_cells,nlev))
135  allocate(ref_field_3D_p1(nb_cells,nlev_p1),field_3D_p1(nb_cells,nlev_p1))
136  allocate(ref_field_soil(nb_cells,nsoil),field_soil(nb_cells,nsoil))
137  allocate(ref_field_vertex(nvertex,nb_cells),field_vertex(nvertex,nb_cells))
138  write(*,*) "dimids :", dimids
139  ! loop over variables:
140  do k=1,nvar
141    varname_dimids(:)=0
142    ierr=NF90_INQUIRE_VARIABLE(ncid,k,name=varname,ndims=ndim,dimids=varname_dimids)
143      write(*,*) "name: ",trim(varname),"; dimensions:",varname_dimids
144    if (ierr /= nf90_noerr) then
145      write(*,*) "error for variable k=",k
146      write(*,*) nf90_strerror(ierr)
147    endif
148    ! Processing 2D variables
149    if ((ndim==1).and.(varname_dimids(1)==dimids(1))) then
150      write(*,*) "processing ",trim(varname)
151      ! load field
152      ierr=NF90_INQ_VARID(ncid,varname,varid)
153      ierr=NF90_GET_VAR(ncid,varid,ref_field)
154      ! rearrange field
155      do i=1,nb_cells
156        field(cell_index(i))=ref_field(i)
157      enddo
158      ! write field
159      ierr=NF90_PUT_VAR(ncid,varid,field)
160
161    ! Processing 3D variables : bounds_lon and bounds_lat
162    else if ((ndim==2).and.(varname_dimids(1)==dimids(4))) then
163      write(*,*) "processing ",trim(varname)
164      ! load field_3D
165      ierr=NF90_INQ_VARID(ncid,varname,varid)
166      if (ierr /= nf90_noerr) print*,"error inqvarid ",trim(varname)
167      ierr=NF90_GET_VAR(ncid,varid,ref_field_vertex)
168      if (ierr /= nf90_noerr) print*,"error get_var ",trim(varname)
169      ! rearrange field_3D
170      do i=1,nb_cells
171        field_vertex(:,cell_index(i))=ref_field_vertex(:,i)
172      enddo
173      ! write field_3D
174      ierr=NF90_PUT_VAR(ncid,varid,field_vertex)
175      if (ierr /= nf90_noerr) print*,"error putvar ",trim(varname)
176 
177    ! Processing 3D variables : altitude
178    else if ((ndim==2).and.(varname_dimids(2)==dimids(2))) then
179      write(*,*) "processing ",trim(varname)
180      ! load field_3D
181      ierr=NF90_INQ_VARID(ncid,varname,varid)
182      if (ierr /= nf90_noerr) print*,"error inqvarid ",trim(varname)
183      ierr=NF90_GET_VAR(ncid,varid,ref_field_3D)
184      if (ierr /= nf90_noerr) print*,"error get_var ",trim(varname)
185      ! rearrange field_3D
186      do i=1,nb_cells
187        field_3D(cell_index(i),:)=ref_field_3D(i,:)
188      enddo
189      ! write field_3D
190      ierr=NF90_PUT_VAR(ncid,varid,field_3D)
191      if (ierr /= nf90_noerr) print*,"error putvar ",trim(varname)
192
193    ! Processing 3D variables : soil
194    else if ((ndim==2).and.(varname_dimids(2)==dimids(3))) then
195      write(*,*) "processing ",trim(varname)
196      ! load field_soil
197      ierr=NF90_INQ_VARID(ncid,varname,varid)
198      if (ierr /= nf90_noerr) print*,"error inqvarid ",trim(varname)
199      ierr=NF90_GET_VAR(ncid,varid,ref_field_soil)
200      if (ierr /= nf90_noerr) print*,"error get_var ",trim(varname)
201      ! rearrange field_soil
202      do i=1,nb_cells
203        field_soil(cell_index(i),:)=ref_field_soil(i,:)
204      enddo
205      ! write field_soil
206      ierr=NF90_PUT_VAR(ncid,varid,field_soil)
207      if (ierr /= nf90_noerr) print*,"error putvar ",trim(varname)
208
209    ! Processing 3D variables : q2 with lev_p1 as vertical axis
210    else if ((ndim==2).and.(varname_dimids(2)==dimids(5))) then
211      write(*,*) "processing ",trim(varname)
212      ! load field_3D
213      ierr=NF90_INQ_VARID(ncid,varname,varid)
214      if (ierr /= nf90_noerr) print*,"error inqvarid ",trim(varname)
215      ierr=NF90_GET_VAR(ncid,varid,ref_field_3D_p1)
216      if (ierr /= nf90_noerr) print*,"error get_var ",trim(varname)
217      ! rearrange field_3D
218      do i=1,nb_cells
219        field_3D_p1(cell_index(i),:)=ref_field_3D_p1(i,:)
220      enddo
221      ! write field_3D
222      ierr=NF90_PUT_VAR(ncid,varid,field_3D_p1)
223      if (ierr /= nf90_noerr) print*,"error putvar ",trim(varname)
224    endif
225  enddo
226 
227END PROGRAM rearrange_startphy
Note: See TracBrowser for help on using the repository browser.