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

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

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