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

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

devel : towards Fortran driver for unstructured/LAM meshes
File size: 8.5 KB
Line 
1MODULE init_unstructured_mod
2  USE mpipara, ONLY : is_mpi_master
3  USE data_unstructured_mod
4  IMPLICIT NONE
5  SAVE
6  PRIVATE
7
8  ! these buffers are used when reading from the grid file
9  REAL(8), ALLOCATABLE :: Ddata_read1(:),Ddata_read2(:,:),Ddata_read3(:,:,:)
10  INTEGER, ALLOCATABLE :: Idata_read1(:),Idata_read2(:,:),Idata_read3(:,:,:)
11
12  PUBLIC :: read_dump_partition
13
14CONTAINS
15
16  SUBROUTINE free_data_read()
17    IF(ALLOCATED(Idata_read1)) DEALLOCATE(Idata_read1) 
18    IF(ALLOCATED(Idata_read2)) DEALLOCATE(Idata_read2) 
19    IF(ALLOCATED(Idata_read3)) DEALLOCATE(Idata_read3) 
20    IF(ALLOCATED(Ddata_read1)) DEALLOCATE(Ddata_read1) 
21    IF(ALLOCATED(Ddata_read2)) DEALLOCATE(Ddata_read2) 
22    IF(ALLOCATED(Ddata_read3)) DEALLOCATE(Ddata_read3) 
23  END SUBROUTINE free_data_read
24
25  SUBROUTINE read_from_gridfile(id_nc, data_type, name)
26    use netcdf_mod
27    CHARACTER(*)    :: data_type, name
28    INTEGER :: id_nc, id_var, status
29    INTEGER :: dim_1, dim_2, dim_3
30    INTEGER :: numDims, dimIDs(3) !max_var_dims
31
32    !-------------------Reading variable-------------------------------
33    status = nf90_inq_varid(id_nc, name,id_var)
34    IF(status /= 0) THEN
35        print *, "Not able to read variable from gridfile :", name
36        STOP "Exit"
37    ENDIF
38    !inquire dimension
39    status = nf90_inquire_variable(id_nc,id_var,dimids=dimIDs,ndims=numDims)
40    status = nf90_inquire_dimension(id_nc, dimIDs(1), len = dim_1)
41    status = nf90_inquire_dimension(id_nc, dimIDs(2), len = dim_2)
42    status = nf90_inquire_dimension(id_nc, dimIDs(3), len = dim_3)
43    SELECT CASE(numDims)
44    CASE(3)
45    print *,"Size of array, ",name,":", dim_1,dim_2,dim_3
46    CASE(2)
47    print *,"Size of array, ",name,":", dim_1,dim_2
48    CASE DEFAULT
49    print *,"Size of array, ",name,":", dim_1
50    END SELECT
51
52    CALL free_data_read
53    SELECT CASE(data_type)
54    CASE('integer')
55       SELECT CASE(numDims)
56       CASE(3)
57          allocate(Idata_read3(dim_1,dim_2,dim_3))
58          status = nf90_get_var(id_nc, id_var,Idata_read3)
59          print *,"First value of array, ",name,":", Idata_read3(1,1,1)
60       CASE(2)
61          allocate(Idata_read2(dim_1,dim_2))
62          status = nf90_get_var(id_nc, id_var,Idata_read2)
63          print *,"First value of array, ",name,":", Idata_read2(1,1)
64       CASE DEFAULT
65          allocate(Idata_read1(dim_1))
66          status = nf90_get_var(id_nc, id_var,Idata_read1)
67          print *,"First value of array, ",name,":", Idata_read1(1)
68       END SELECT
69    CASE DEFAULT
70       SELECT CASE(numDims)
71       CASE(3)
72          allocate(Ddata_read3(dim_1,dim_2,dim_3))
73          status = nf90_get_var(id_nc, id_var,Ddata_read3)
74          print *,"First value of array, ",name,":", Ddata_read3(1,1,1)
75       CASE(2)
76          allocate(Ddata_read2(dim_1,dim_2))
77          status = nf90_get_var(id_nc, id_var,Ddata_read2)
78          print *,"First value of array, ",name,":", Ddata_read2(1,1)
79       CASE DEFAULT
80          allocate(Ddata_read1(dim_1))
81          status = nf90_get_var(id_nc, id_var,Ddata_read1)
82          print *,"First value of array, ",name,":", Ddata_read1(1)
83       END SELECT
84    END SELECT
85   
86    IF(status /= nf90_NoErr) THEN
87        print *, "Error when reading from grid file : ", name
88        STOP "Exit"
89    ENDIF
90
91  END SUBROUTINE read_from_gridfile
92
93
94  SUBROUTINE read_dump_partition
95    use netcdf_mod
96    USE ioipsl
97    USE field_mod
98    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e
99    IMPLICIT NONE
100
101    !!-------------Declare local variables-------------------
102    CHARACTER(LEN=*),PARAMETER :: f="input/mesh_information.nc"
103    INTEGER :: id_nc, ierr, status, descriptionLength, ij
104    CHARACTER(LEN= 80) :: description
105    REAL(rstd), ALLOCATABLE :: angle_e(:)
106    REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat
107         elon(3), elat(3) ! lon/lat unit vectors
108    print *,"------------------ READING FILE " , f, "----------------------- "
109    !open and read the input file
110    ierr = nf90_open(f, NF90_NOWRITE, id_nc)
111    if (ierr /= nf90_noerr) then
112      print *, trim(nf90_strerror(ierr))
113      stop "Error reading file"
114    end if
115
116    status = nf90_inquire_attribute(id_nc, nf90_global, "description", len =descriptionLength)
117    IF(status /= 0 .or. len(description) < descriptionLength) THEN
118        print *, "Not enough space to put NetCDF attribute values."
119        STOP "Error reading file"
120    ENDIF
121
122    !-------------------Reading global attributes-----------------------
123    status = nf90_get_att(id_nc, nf90_global, "description", description)
124    print *,"Data file description :",description
125
126    !status = nf90_get_att(id_nc, nf90_global, "primal_num", primal_num)
127    !status = nf90_get_att(id_nc, nf90_global, "dual_num", dual_num)
128    !status = nf90_get_att(id_nc, nf90_global, "edge_num", edge_num)
129    !print *,primal_num,dual_num,edge_num
130
131    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
132    ALLOCATE(primal_deg, source = Idata_read1)
133    CALL read_from_gridfile(id_nc, 'integer', 'primal_edge')
134    ALLOCATE(primal_edge, source = Idata_read2)
135    CALL read_from_gridfile(id_nc, 'integer', 'primal_ne')
136    ALLOCATE(primal_ne, source = Idata_read2)
137    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
138    ALLOCATE(dual_deg, source = Idata_read1)
139    CALL read_from_gridfile(id_nc, 'integer', 'dual_edge')
140    ALLOCATE(dual_edge, source = Idata_read2)
141    CALL read_from_gridfile(id_nc, 'integer', 'dual_ne')
142    ALLOCATE(dual_ne, source = Idata_read2)
143    CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') 
144    ALLOCATE(primal_vertex, source = Idata_read2)
145    CALL read_from_gridfile(id_nc, 'integer', 'left')
146    ALLOCATE(left, source = Idata_read1)
147    CALL read_from_gridfile(id_nc, 'integer', 'right')
148    ALLOCATE(right, source = Idata_read1)
149    CALL read_from_gridfile(id_nc, 'integer', 'up')
150    ALLOCATE(up, source = Idata_read1)
151    CALL read_from_gridfile(id_nc, 'integer', 'down')
152    ALLOCATE(down, source = Idata_read1)
153    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
154    ALLOCATE(trisk_deg, source = Idata_read1)
155    CALL read_from_gridfile(id_nc, 'integer', 'trisk')
156    ALLOCATE(trisk, source = Idata_read2)
157    CALL read_from_gridfile(id_nc, 'float', 'Ai')
158    ALLOCATE(Ai, source = Ddata_read1)
159    CALL read_from_gridfile(id_nc, 'float', 'Av')
160    ALLOCATE(Av, source = Ddata_read1)
161    CALL read_from_gridfile(id_nc, 'float', 'le_de')
162    ALLOCATE(le_de, source = Ddata_read1)
163    CALL read_from_gridfile(id_nc, 'float', 'Riv2')
164    ALLOCATE(Riv2, source = Ddata_read2)
165    CALL read_from_gridfile(id_nc, 'float', 'wee')
166    ALLOCATE(wee, source = Ddata_read2)
167    CALL read_from_gridfile(id_nc, 'float', 'fv') !
168    ALLOCATE(fv, source = Ddata_read1)
169    CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') 
170    ALLOCATE(dual_vertex, source = Idata_read2)
171
172    CALL read_from_gridfile(id_nc, 'float', 'lon_i')
173    ALLOCATE(lon_i, source = Ddata_read1)
174    CALL read_from_gridfile(id_nc, 'float', 'lat_i')
175    ALLOCATE(lat_i, source = Ddata_read1)
176    CALL read_from_gridfile(id_nc, 'float', 'lon_e')
177    ALLOCATE(lon_e, source = Ddata_read1)
178    CALL read_from_gridfile(id_nc, 'float', 'lat_e')
179    ALLOCATE(lat_e, source = Ddata_read1)
180    CALL read_from_gridfile(id_nc, 'float', 'angle_e')
181    ALLOCATE(angle_e, source = Ddata_read1)
182
183    CALL free_data_read ! free buffers after reading all data from grid file
184
185    edge_num = SIZE(le_de)
186    primal_num = SIZE(Ai)
187    dual_num = SIZE(Av)
188    max_primal_deg = SIZE(primal_edge,1)
189    max_dual_deg = SIZE(dual_edge,1)
190    max_trisk_deg = SIZE(trisk,1)
191
192    ! now post-process some of the data we just read
193    ALLOCATE(ep_e(edge_num,3))
194    DO ij=1,edge_num
195       clon = COS(lon_e(ij))
196       slon = SIN(lon_e(ij))
197       clat = COS(lat_e(ij))
198       slat = SIN(lat_e(ij))
199       ! x,y,z = clat*clon, clat*slon, slat
200       elon = (/ -slon, clon, 0. /)
201       elat = (/ -slat*clon, -slat*slon, clat /)
202       ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
203    END DO
204   
205    DEALLOCATE(angle_e)
206  END SUBROUTINE read_dump_partition
207
208  SUBROUTINE init_grid_type
209    USE grid_param, ONLY : grid_type, grid_unst, grid_ico   
210    USE getin_mod, ONLY : getin
211    CHARACTER(len=255) :: grid_type_var
212    grid_type_var='icosahedral'
213    CALL getin("grid_type",grid_type_var)
214    SELECT CASE(grid_type_var)
215    CASE('unstructured')
216       grid_type = grid_unst
217!       is_omp_level_master=.TRUE.
218!       omp_level_size=1
219       CALL read_dump_partition
220       IF (is_mpi_master) PRINT *,'Using unstructured grid type'
221    CASE DEFAULT
222       grid_type = grid_ico
223       IF (is_mpi_master) PRINT *,'Using default grid type'
224    END SELECT 
225  END SUBROUTINE init_grid_type
226
227END MODULE init_unstructured_mod
Note: See TracBrowser for help on using the repository browser.