source: trunk/LMDZ.PLUTO/util/startarchive2icosa/rearrange_startphy.f90 @ 3545

Last change on this file since 3545 was 3545, checked in by afalco, 36 hours ago

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