source: codes/icosagcm/devel/src/unstructured/init_unstructured.f90 @ 879

Last change on this file since 879 was 879, checked in by dubos, 5 years ago

devel : introduced derived type to store cell bounds
File size: 10.1 KB
Line 
1MODULE init_unstructured_mod
2  USE mpipara, ONLY : is_mpi_master
3  USE data_unstructured_mod
4  USE geometry, ONLY : swap_geometry
5  IMPLICIT NONE
6  SAVE
7  PRIVATE
8
9  ! these buffers are used when reading from the grid file
10  REAL(8), ALLOCATABLE :: Ddata_read1(:),Ddata_read2(:,:),Ddata_read3(:,:,:)
11  INTEGER, ALLOCATABLE :: Idata_read1(:),Idata_read2(:,:),Idata_read3(:,:,:)
12
13  CHARACTER(LEN=*),PARAMETER :: meshfile="input/mesh_information.nc"
14  INTEGER :: id_nc ! NetCDF id of mesh file open by open_local_mesh_file
15
16  PUBLIC :: open_local_mesh_file, read_local_mesh
17
18CONTAINS
19
20  SUBROUTINE free_data_read()
21    IF(ALLOCATED(Idata_read1)) DEALLOCATE(Idata_read1) 
22    IF(ALLOCATED(Idata_read2)) DEALLOCATE(Idata_read2) 
23    IF(ALLOCATED(Idata_read3)) DEALLOCATE(Idata_read3) 
24    IF(ALLOCATED(Ddata_read1)) DEALLOCATE(Ddata_read1) 
25    IF(ALLOCATED(Ddata_read2)) DEALLOCATE(Ddata_read2) 
26    IF(ALLOCATED(Ddata_read3)) DEALLOCATE(Ddata_read3) 
27  END SUBROUTINE free_data_read
28
29  SUBROUTINE read_from_gridfile(id_nc, data_type, name)
30    use netcdf_mod
31    CHARACTER(*)    :: data_type, name
32    INTEGER :: id_nc, id_var, status
33    INTEGER :: dim_1, dim_2, dim_3
34    INTEGER :: numDims, dimIDs(3) !max_var_dims
35
36    !-------------------Reading variable-------------------------------
37    status = nf90_inq_varid(id_nc, name,id_var)
38    IF(status /= 0) THEN
39        print *, "Not able to read variable from gridfile :", name
40        STOP "Exit"
41    ENDIF
42    !inquire dimension
43    status = nf90_inquire_variable(id_nc,id_var,dimids=dimIDs,ndims=numDims)
44    status = nf90_inquire_dimension(id_nc, dimIDs(1), len = dim_1)
45    status = nf90_inquire_dimension(id_nc, dimIDs(2), len = dim_2)
46    status = nf90_inquire_dimension(id_nc, dimIDs(3), len = dim_3)
47    SELECT CASE(numDims)
48    CASE(3)
49    print *,"Size of array, ",name,":", dim_1,dim_2,dim_3
50    CASE(2)
51    print *,"Size of array, ",name,":", dim_1,dim_2
52    CASE DEFAULT
53    print *,"Size of array, ",name,":", dim_1
54    END SELECT
55
56    CALL free_data_read
57    SELECT CASE(data_type)
58    CASE('integer')
59       SELECT CASE(numDims)
60       CASE(3)
61          allocate(Idata_read3(dim_1,dim_2,dim_3))
62          status = nf90_get_var(id_nc, id_var,Idata_read3)
63          print *,"First value of array, ",name,":", Idata_read3(1,1,1)
64       CASE(2)
65          allocate(Idata_read2(dim_1,dim_2))
66          status = nf90_get_var(id_nc, id_var,Idata_read2)
67          print *,"First value of array, ",name,":", Idata_read2(1,1)
68       CASE DEFAULT
69          allocate(Idata_read1(dim_1))
70          status = nf90_get_var(id_nc, id_var,Idata_read1)
71          print *,"First value of array, ",name,":", Idata_read1(1)
72       END SELECT
73    CASE DEFAULT
74       SELECT CASE(numDims)
75       CASE(3)
76          allocate(Ddata_read3(dim_1,dim_2,dim_3))
77          status = nf90_get_var(id_nc, id_var,Ddata_read3)
78          print *,"First value of array, ",name,":", Ddata_read3(1,1,1)
79       CASE(2)
80          allocate(Ddata_read2(dim_1,dim_2))
81          status = nf90_get_var(id_nc, id_var,Ddata_read2)
82          print *,"First value of array, ",name,":", Ddata_read2(1,1)
83       CASE DEFAULT
84          allocate(Ddata_read1(dim_1))
85          status = nf90_get_var(id_nc, id_var,Ddata_read1)
86          print *,"First value of array, ",name,":", Ddata_read1(1)
87       END SELECT
88    END SELECT
89   
90    IF(status /= nf90_NoErr) THEN
91        print *, "Error when reading from grid file : ", name
92        STOP "Exit"
93    ENDIF
94
95  END SUBROUTINE read_from_gridfile
96
97
98  SUBROUTINE open_local_mesh_file
99    USE netcdf_mod
100    CHARACTER(LEN= 80) :: description
101    INTEGER :: ierr, status, descriptionLength
102
103    PRINT *,"------------------ READING FILE " , meshfile, "----------------------- "
104    !open and read the input file
105    ierr = nf90_open(meshfile, NF90_NOWRITE, id_nc)
106    if (ierr /= nf90_noerr) then
107      print *, trim(nf90_strerror(ierr))
108      stop "Error reading file"
109    end if
110
111    status = nf90_inquire_attribute(id_nc, nf90_global, "description", len=descriptionLength)
112    IF(status /= 0 .or. len(description) < descriptionLength) THEN
113        print *, "Not enough space to put NetCDF attribute values."
114        STOP "Error reading file"
115    ENDIF
116
117    !-------------------Reading global attributes-----------------------
118    status = nf90_get_att(id_nc, nf90_global, "description", description)
119    print *,"Data file description :",description
120
121    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
122    primal_num = SIZE(Idata_read1)
123    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
124    dual_num = SIZE(Idata_read1)
125    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
126    edge_num = SIZE(Idata_read1)
127  END SUBROUTINE open_local_mesh_file
128
129
130  SUBROUTINE read_field(id_nc, field)
131    USE field_mod, ONLY : t_field, type_integer, type_real
132    INTEGER :: id_nc
133    TYPE(t_field), POINTER :: field(:), fld
134    fld=>field(1)
135    SELECT CASE(fld%data_type)
136    CASE(type_integer)
137       CALL read_from_gridfile(id_nc, 'integer', TRIM(fld%name))
138       SELECT CASE(field(1)%ndim)
139       CASE(2)
140          fld%ival2d = Idata_read1
141       CASE(3)
142          fld%ival3d = Idata_read2
143       END SELECT
144    CASE(type_real)
145       CALL read_from_gridfile(id_nc, 'float', TRIM(fld%name))
146       SELECT CASE(field(1)%ndim)
147       CASE(2)
148          fld%rval2d = Ddata_read1
149       CASE(3)
150          fld%rval3d = Ddata_read2
151       END SELECT
152    END SELECT
153  END SUBROUTINE read_field
154
155  SUBROUTINE read_local_mesh
156    USE field_mod
157    USE geometry, ONLY : geom, lon_i, lat_i, lon_e, lat_e, ep_e
158    USE domain_mod, ONLY : domain, domain_glo, t_domain, swap_needed
159    IMPLICIT NONE
160    INTEGER :: ij
161    INTEGER, ALLOCATABLE :: cell_ij(:)
162    REAL(rstd), ALLOCATABLE :: angle_e(:)
163    REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat
164         elon(3), elat(3) ! lon/lat unit vectors
165   
166    ! if possible data is read into fields allocated in geom(1) using read_field
167    CALL read_field(id_nc, geom%Ai)
168    CALL read_field(id_nc, geom%Av)
169    CALL read_field(id_nc, geom%fv)
170    CALL read_field(id_nc, geom%le_de)
171    CALL read_field(id_nc, geom%le)
172    CALL read_field(id_nc, geom%de)
173    CALL read_field(id_nc, geom%lon_i)
174    CALL read_field(id_nc, geom%lat_i)
175    CALL read_field(id_nc, geom%lon_e)
176    CALL read_field(id_nc, geom%lat_e)
177
178    swap_needed = .TRUE.
179    CALL swap_geometry(1)
180
181    ! cell centers and bounds are read into domain(1)%primal_own
182    domain(1)%primal_own%ncell = primal_num
183    domain_glo(1)%primal_own%ncell = primal_num
184    ALLOCATE(domain(1)%primal_own%lon,     source=lon_i)
185    ALLOCATE(domain_glo(1)%primal_own%lon, source=lon_i)
186    ALLOCATE(domain(1)%primal_own%lat,     source=lat_i)
187    ALLOCATE(domain_glo(1)%primal_own%lat, source=lat_i)
188
189    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lon')
190    ALLOCATE(domain(1)%primal_own%bnds_lon,     source=Ddata_read2)
191    ALLOCATE(domain_glo(1)%primal_own%bnds_lon, source=Ddata_read2)
192    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lat')
193    ALLOCATE(domain(1)%primal_own%bnds_lat,     source=Ddata_read2)
194    ALLOCATE(domain_glo(1)%primal_own%bnds_lat, source=Ddata_read2)
195    PRINT *, 'read_local_mesh : primal_num =', primal_num, domain_glo(1)%primal_own%ncell
196
197    ! other fields are defined only in the case of an unstructured mesh
198    ! those fields are read into simple arrays
199    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
200    ALLOCATE(primal_deg, source = Idata_read1)
201    CALL read_from_gridfile(id_nc, 'integer', 'primal_edge')
202    ALLOCATE(primal_edge, source = Idata_read2)
203    CALL read_from_gridfile(id_nc, 'integer', 'primal_ne')
204    ALLOCATE(primal_ne, source = Idata_read2)
205    CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') 
206    ALLOCATE(primal_vertex, source = Idata_read2)
207
208    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
209    ALLOCATE(dual_deg, source = Idata_read1)
210    CALL read_from_gridfile(id_nc, 'integer', 'dual_edge')
211    ALLOCATE(dual_edge, source = Idata_read2)
212    CALL read_from_gridfile(id_nc, 'integer', 'dual_ne')
213    ALLOCATE(dual_ne, source = Idata_read2)
214    CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') 
215    ALLOCATE(dual_vertex, source = Idata_read2)
216
217    CALL read_from_gridfile(id_nc, 'integer', 'left')
218    ALLOCATE(left, source = Idata_read1)
219    CALL read_from_gridfile(id_nc, 'integer', 'right')
220    ALLOCATE(right, source = Idata_read1)
221    CALL read_from_gridfile(id_nc, 'integer', 'up')
222    ALLOCATE(up, source = Idata_read1)
223    CALL read_from_gridfile(id_nc, 'integer', 'down')
224    ALLOCATE(down, source = Idata_read1)
225    CALL read_from_gridfile(id_nc, 'float', 'angle_e')
226    ALLOCATE(angle_e, source = Ddata_read1)
227    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
228    ALLOCATE(trisk_deg, source = Idata_read1)
229    CALL read_from_gridfile(id_nc, 'integer', 'trisk')
230    ALLOCATE(trisk, source = Idata_read2)
231    CALL read_from_gridfile(id_nc, 'float', 'Riv2')
232    ALLOCATE(Riv2, source = Ddata_read2)
233    CALL read_from_gridfile(id_nc, 'float', 'wee')
234    ALLOCATE(wee, source = Ddata_read2)
235
236    CALL free_data_read ! free buffers after reading all data from grid file
237
238
239!    edge_num = SIZE(le_de)
240!    primal_num = SIZE(Ai)
241!    dual_num = SIZE(Av)
242    max_primal_deg = SIZE(primal_edge,1)
243    max_dual_deg = SIZE(dual_edge,1)
244    max_trisk_deg = SIZE(trisk,1)
245
246    ! now post-process some of the data we just read
247    ALLOCATE(ep_e(edge_num,3))
248    DO ij=1,edge_num
249       clon = COS(lon_e(ij))
250       slon = SIN(lon_e(ij))
251       clat = COS(lat_e(ij))
252       slat = SIN(lat_e(ij))
253       ! x,y,z = clat*clon, clat*slon, slat
254       elon = (/ -slon, clon, 0. /)
255       elat = (/ -slat*clon, -slat*slon, clat /)
256       ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
257    END DO
258   
259    DEALLOCATE(angle_e)
260
261    ALLOCATE(domain(1)%primal_own%ij(primal_num))
262    domain(1)%primal_own%ij(:) = (/(ij,ij=1,primal_num)/)
263    PRINT *, 'read_local_mesh :', SHAPE(domain), primal_num, &
264         SHAPE(domain(1)%primal_own%ij)
265
266    ALLOCATE(cell_ij(dual_num))
267    cell_ij(:) = (/(ij,ij=1,dual_num)/)
268    ALLOCATE(domain(1)%dual_own%ij, source=cell_ij)
269    DEALLOCATE(cell_ij)
270
271    CALL swap_geometry(1)
272    swap_needed = .FALSE.
273  END SUBROUTINE read_local_mesh
274
275END MODULE init_unstructured_mod
Note: See TracBrowser for help on using the repository browser.