[813] | 1 | MODULE init_unstructured_mod |
---|
[883] | 2 | USE prec, ONLY : rstd |
---|
| 3 | USE math_const, ONLY : radian_to_degree |
---|
[813] | 4 | USE mpipara, ONLY : is_mpi_master |
---|
| 5 | USE data_unstructured_mod |
---|
[879] | 6 | USE geometry, ONLY : swap_geometry |
---|
[813] | 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 | |
---|
[863] | 15 | CHARACTER(LEN=*),PARAMETER :: meshfile="input/mesh_information.nc" |
---|
| 16 | INTEGER :: id_nc ! NetCDF id of mesh file open by open_local_mesh_file |
---|
[813] | 17 | |
---|
[863] | 18 | PUBLIC :: open_local_mesh_file, read_local_mesh |
---|
| 19 | |
---|
[813] | 20 | CONTAINS |
---|
| 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 | |
---|
[863] | 99 | SUBROUTINE open_local_mesh_file |
---|
| 100 | USE netcdf_mod |
---|
| 101 | CHARACTER(LEN= 80) :: description |
---|
| 102 | INTEGER :: ierr, status, descriptionLength |
---|
[813] | 103 | |
---|
[863] | 104 | PRINT *,"------------------ READING FILE " , meshfile, "----------------------- " |
---|
[813] | 105 | !open and read the input file |
---|
[863] | 106 | ierr = nf90_open(meshfile, NF90_NOWRITE, id_nc) |
---|
[813] | 107 | if (ierr /= nf90_noerr) then |
---|
| 108 | print *, trim(nf90_strerror(ierr)) |
---|
| 109 | stop "Error reading file" |
---|
| 110 | end if |
---|
| 111 | |
---|
[863] | 112 | status = nf90_inquire_attribute(id_nc, nf90_global, "description", len=descriptionLength) |
---|
[813] | 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 | |
---|
[863] | 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 | |
---|
[879] | 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 | |
---|
[883] | 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 | |
---|
[863] | 218 | SUBROUTINE read_local_mesh |
---|
| 219 | USE field_mod |
---|
[883] | 220 | USE domain_mod, ONLY : swap_needed, domain, domain_glo |
---|
[879] | 221 | USE geometry, ONLY : geom, lon_i, lat_i, lon_e, lat_e, ep_e |
---|
[883] | 222 | USE netcdf_mod, ONLY : nf90_close |
---|
[863] | 223 | IMPLICIT NONE |
---|
[883] | 224 | INTEGER :: ij, ierr |
---|
[879] | 225 | INTEGER, ALLOCATABLE :: cell_ij(:) |
---|
[863] | 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 | |
---|
[879] | 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) |
---|
[813] | 241 | |
---|
[879] | 242 | swap_needed = .TRUE. |
---|
| 243 | CALL swap_geometry(1) |
---|
[883] | 244 | |
---|
[879] | 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 |
---|
[813] | 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) |
---|
[879] | 255 | CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') |
---|
| 256 | ALLOCATE(primal_vertex, source = Idata_read2) |
---|
| 257 | |
---|
[813] | 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) |
---|
[879] | 264 | CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') |
---|
| 265 | ALLOCATE(dual_vertex, source = Idata_read2) |
---|
| 266 | |
---|
[813] | 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) |
---|
[879] | 275 | CALL read_from_gridfile(id_nc, 'float', 'angle_e') |
---|
| 276 | ALLOCATE(angle_e, source = Ddata_read1) |
---|
[813] | 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, source = Ddata_read2) |
---|
| 285 | |
---|
[883] | 286 | ! read cell centers and bounds |
---|
| 287 | CALL read_local_mesh_bounds(domain(1)) |
---|
| 288 | CALL read_local_mesh_bounds(domain_glo(1)) |
---|
| 289 | |
---|
[813] | 290 | CALL free_data_read ! free buffers after reading all data from grid file |
---|
[883] | 291 | ierr = nf90_close(id_nc) |
---|
[813] | 292 | |
---|
[883] | 293 | ! now process some of the data we just read |
---|
[879] | 294 | |
---|
[813] | 295 | max_primal_deg = SIZE(primal_edge,1) |
---|
| 296 | max_dual_deg = SIZE(dual_edge,1) |
---|
| 297 | max_trisk_deg = SIZE(trisk,1) |
---|
| 298 | |
---|
[856] | 299 | DO ij=1,edge_num |
---|
| 300 | clon = COS(lon_e(ij)) |
---|
| 301 | slon = SIN(lon_e(ij)) |
---|
| 302 | clat = COS(lat_e(ij)) |
---|
| 303 | slat = SIN(lat_e(ij)) |
---|
| 304 | ! x,y,z = clat*clon, clat*slon, slat |
---|
| 305 | elon = (/ -slon, clon, 0. /) |
---|
| 306 | elat = (/ -slat*clon, -slat*slon, clat /) |
---|
[883] | 307 | ep_e(:,ij) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat |
---|
[856] | 308 | END DO |
---|
| 309 | |
---|
| 310 | DEALLOCATE(angle_e) |
---|
[879] | 311 | |
---|
| 312 | CALL swap_geometry(1) |
---|
| 313 | swap_needed = .FALSE. |
---|
[883] | 314 | |
---|
[863] | 315 | END SUBROUTINE read_local_mesh |
---|
[813] | 316 | |
---|
| 317 | END MODULE init_unstructured_mod |
---|