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

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

devel : split read_dump_partition into open_local_mesh_file and read_local_mesh
File size: 8.5 KB
Line 
1MODULE init_unstructured_mod
2  USE mpipara, ONLY : is_mpi_master
3  USE data_unstructured_mod
4  USE geometry, ONLY : de
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_local_mesh
131    USE field_mod
132    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e
133    IMPLICIT NONE
134    INTEGER :: ij
135    REAL(rstd), ALLOCATABLE :: angle_e(:)
136    REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat
137         elon(3), elat(3) ! lon/lat unit vectors
138   
139    !status = nf90_get_att(id_nc, nf90_global, "primal_num", primal_num)
140    !status = nf90_get_att(id_nc, nf90_global, "dual_num", dual_num)
141    !status = nf90_get_att(id_nc, nf90_global, "edge_num", edge_num)
142    !print *,primal_num,dual_num,edge_num
143
144    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
145    ALLOCATE(primal_deg, source = Idata_read1)
146    CALL read_from_gridfile(id_nc, 'integer', 'primal_edge')
147    ALLOCATE(primal_edge, source = Idata_read2)
148    CALL read_from_gridfile(id_nc, 'integer', 'primal_ne')
149    ALLOCATE(primal_ne, source = Idata_read2)
150    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
151    ALLOCATE(dual_deg, source = Idata_read1)
152    CALL read_from_gridfile(id_nc, 'integer', 'dual_edge')
153    ALLOCATE(dual_edge, source = Idata_read2)
154    CALL read_from_gridfile(id_nc, 'integer', 'dual_ne')
155    ALLOCATE(dual_ne, source = Idata_read2)
156    CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') 
157    ALLOCATE(primal_vertex, source = Idata_read2)
158    CALL read_from_gridfile(id_nc, 'integer', 'left')
159    ALLOCATE(left, source = Idata_read1)
160    CALL read_from_gridfile(id_nc, 'integer', 'right')
161    ALLOCATE(right, source = Idata_read1)
162    CALL read_from_gridfile(id_nc, 'integer', 'up')
163    ALLOCATE(up, source = Idata_read1)
164    CALL read_from_gridfile(id_nc, 'integer', 'down')
165    ALLOCATE(down, source = Idata_read1)
166    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
167    ALLOCATE(trisk_deg, source = Idata_read1)
168    CALL read_from_gridfile(id_nc, 'integer', 'trisk')
169    ALLOCATE(trisk, source = Idata_read2)
170    CALL read_from_gridfile(id_nc, 'float', 'Ai')
171    ALLOCATE(Ai, source = Ddata_read1)
172    CALL read_from_gridfile(id_nc, 'float', 'Av')
173    ALLOCATE(Av, source = Ddata_read1)
174    CALL read_from_gridfile(id_nc, 'float', 'le_de')
175    ALLOCATE(le_de, source = Ddata_read1)
176    CALL read_from_gridfile(id_nc, 'float', 'le')
177    ALLOCATE(le, source = Ddata_read1)
178    CALL read_from_gridfile(id_nc, 'float', 'de')
179    ALLOCATE(de, source = Ddata_read1)
180    CALL read_from_gridfile(id_nc, 'float', 'Riv2')
181    ALLOCATE(Riv2, source = Ddata_read2)
182    CALL read_from_gridfile(id_nc, 'float', 'wee')
183    ALLOCATE(wee, source = Ddata_read2)
184    CALL read_from_gridfile(id_nc, 'float', 'fv') !
185    ALLOCATE(fv, source = Ddata_read1)
186    CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') 
187    ALLOCATE(dual_vertex, source = Idata_read2)
188
189    CALL read_from_gridfile(id_nc, 'float', 'lon_i')
190    ALLOCATE(lon_i, source = Ddata_read1)
191    CALL read_from_gridfile(id_nc, 'float', 'lat_i')
192    ALLOCATE(lat_i, source = Ddata_read1)
193    CALL read_from_gridfile(id_nc, 'float', 'lon_e')
194    ALLOCATE(lon_e, source = Ddata_read1)
195    CALL read_from_gridfile(id_nc, 'float', 'lat_e')
196    ALLOCATE(lat_e, source = Ddata_read1)
197    CALL read_from_gridfile(id_nc, 'float', 'angle_e')
198    ALLOCATE(angle_e, source = Ddata_read1)
199
200    CALL free_data_read ! free buffers after reading all data from grid file
201
202    edge_num = SIZE(le_de)
203    primal_num = SIZE(Ai)
204    dual_num = SIZE(Av)
205    max_primal_deg = SIZE(primal_edge,1)
206    max_dual_deg = SIZE(dual_edge,1)
207    max_trisk_deg = SIZE(trisk,1)
208
209    ! now post-process some of the data we just read
210    ALLOCATE(ep_e(edge_num,3))
211    DO ij=1,edge_num
212       clon = COS(lon_e(ij))
213       slon = SIN(lon_e(ij))
214       clat = COS(lat_e(ij))
215       slat = SIN(lat_e(ij))
216       ! x,y,z = clat*clon, clat*slon, slat
217       elon = (/ -slon, clon, 0. /)
218       elat = (/ -slat*clon, -slat*slon, clat /)
219       ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
220    END DO
221   
222    DEALLOCATE(angle_e)
223  END SUBROUTINE read_local_mesh
224
225END MODULE init_unstructured_mod
Note: See TracBrowser for help on using the repository browser.