MODULE init_unstructured_mod USE mpipara, ONLY : is_mpi_master USE data_unstructured_mod IMPLICIT NONE SAVE PRIVATE ! these buffers are used when reading from the grid file REAL(8), ALLOCATABLE :: Ddata_read1(:),Ddata_read2(:,:),Ddata_read3(:,:,:) INTEGER, ALLOCATABLE :: Idata_read1(:),Idata_read2(:,:),Idata_read3(:,:,:) PUBLIC :: read_dump_partition CONTAINS SUBROUTINE free_data_read() IF(ALLOCATED(Idata_read1)) DEALLOCATE(Idata_read1) IF(ALLOCATED(Idata_read2)) DEALLOCATE(Idata_read2) IF(ALLOCATED(Idata_read3)) DEALLOCATE(Idata_read3) IF(ALLOCATED(Ddata_read1)) DEALLOCATE(Ddata_read1) IF(ALLOCATED(Ddata_read2)) DEALLOCATE(Ddata_read2) IF(ALLOCATED(Ddata_read3)) DEALLOCATE(Ddata_read3) END SUBROUTINE free_data_read SUBROUTINE read_from_gridfile(id_nc, data_type, name) use netcdf_mod CHARACTER(*) :: data_type, name INTEGER :: id_nc, id_var, status INTEGER :: dim_1, dim_2, dim_3 INTEGER :: numDims, dimIDs(3) !max_var_dims !-------------------Reading variable------------------------------- status = nf90_inq_varid(id_nc, name,id_var) IF(status /= 0) THEN print *, "Not able to read variable from gridfile :", name STOP "Exit" ENDIF !inquire dimension status = nf90_inquire_variable(id_nc,id_var,dimids=dimIDs,ndims=numDims) status = nf90_inquire_dimension(id_nc, dimIDs(1), len = dim_1) status = nf90_inquire_dimension(id_nc, dimIDs(2), len = dim_2) status = nf90_inquire_dimension(id_nc, dimIDs(3), len = dim_3) SELECT CASE(numDims) CASE(3) print *,"Size of array, ",name,":", dim_1,dim_2,dim_3 CASE(2) print *,"Size of array, ",name,":", dim_1,dim_2 CASE DEFAULT print *,"Size of array, ",name,":", dim_1 END SELECT CALL free_data_read SELECT CASE(data_type) CASE('integer') SELECT CASE(numDims) CASE(3) allocate(Idata_read3(dim_1,dim_2,dim_3)) status = nf90_get_var(id_nc, id_var,Idata_read3) print *,"First value of array, ",name,":", Idata_read3(1,1,1) CASE(2) allocate(Idata_read2(dim_1,dim_2)) status = nf90_get_var(id_nc, id_var,Idata_read2) print *,"First value of array, ",name,":", Idata_read2(1,1) CASE DEFAULT allocate(Idata_read1(dim_1)) status = nf90_get_var(id_nc, id_var,Idata_read1) print *,"First value of array, ",name,":", Idata_read1(1) END SELECT CASE DEFAULT SELECT CASE(numDims) CASE(3) allocate(Ddata_read3(dim_1,dim_2,dim_3)) status = nf90_get_var(id_nc, id_var,Ddata_read3) print *,"First value of array, ",name,":", Ddata_read3(1,1,1) CASE(2) allocate(Ddata_read2(dim_1,dim_2)) status = nf90_get_var(id_nc, id_var,Ddata_read2) print *,"First value of array, ",name,":", Ddata_read2(1,1) CASE DEFAULT allocate(Ddata_read1(dim_1)) status = nf90_get_var(id_nc, id_var,Ddata_read1) print *,"First value of array, ",name,":", Ddata_read1(1) END SELECT END SELECT IF(status /= nf90_NoErr) THEN print *, "Error when reading from grid file : ", name STOP "Exit" ENDIF END SUBROUTINE read_from_gridfile SUBROUTINE read_dump_partition use netcdf_mod USE ioipsl USE field_mod USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e IMPLICIT NONE !!-------------Declare local variables------------------- CHARACTER(LEN=*),PARAMETER :: f="input/mesh_information.nc" INTEGER :: id_nc, ierr, status, descriptionLength, ij CHARACTER(LEN= 80) :: description REAL(rstd), ALLOCATABLE :: angle_e(:) REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat elon(3), elat(3) ! lon/lat unit vectors print *,"------------------ READING FILE " , f, "----------------------- " !open and read the input file ierr = nf90_open(f, NF90_NOWRITE, id_nc) if (ierr /= nf90_noerr) then print *, trim(nf90_strerror(ierr)) stop "Error reading file" end if status = nf90_inquire_attribute(id_nc, nf90_global, "description", len =descriptionLength) IF(status /= 0 .or. len(description) < descriptionLength) THEN print *, "Not enough space to put NetCDF attribute values." STOP "Error reading file" ENDIF !-------------------Reading global attributes----------------------- status = nf90_get_att(id_nc, nf90_global, "description", description) print *,"Data file description :",description !status = nf90_get_att(id_nc, nf90_global, "primal_num", primal_num) !status = nf90_get_att(id_nc, nf90_global, "dual_num", dual_num) !status = nf90_get_att(id_nc, nf90_global, "edge_num", edge_num) !print *,primal_num,dual_num,edge_num CALL read_from_gridfile(id_nc, 'integer', 'primal_deg') ALLOCATE(primal_deg, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'primal_edge') ALLOCATE(primal_edge, source = Idata_read2) CALL read_from_gridfile(id_nc, 'integer', 'primal_ne') ALLOCATE(primal_ne, source = Idata_read2) CALL read_from_gridfile(id_nc, 'integer', 'dual_deg') ALLOCATE(dual_deg, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'dual_edge') ALLOCATE(dual_edge, source = Idata_read2) CALL read_from_gridfile(id_nc, 'integer', 'dual_ne') ALLOCATE(dual_ne, source = Idata_read2) CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') ALLOCATE(primal_vertex, source = Idata_read2) CALL read_from_gridfile(id_nc, 'integer', 'left') ALLOCATE(left, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'right') ALLOCATE(right, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'up') ALLOCATE(up, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'down') ALLOCATE(down, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg') ALLOCATE(trisk_deg, source = Idata_read1) CALL read_from_gridfile(id_nc, 'integer', 'trisk') ALLOCATE(trisk, source = Idata_read2) CALL read_from_gridfile(id_nc, 'float', 'Ai') ALLOCATE(Ai, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'Av') ALLOCATE(Av, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'le_de') ALLOCATE(le_de, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'Riv2') ALLOCATE(Riv2, source = Ddata_read2) CALL read_from_gridfile(id_nc, 'float', 'wee') ALLOCATE(wee, source = Ddata_read2) CALL read_from_gridfile(id_nc, 'float', 'fv') ! ALLOCATE(fv, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') ALLOCATE(dual_vertex, source = Idata_read2) CALL read_from_gridfile(id_nc, 'float', 'lon_i') ALLOCATE(lon_i, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'lat_i') ALLOCATE(lat_i, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'lon_e') ALLOCATE(lon_e, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'lat_e') ALLOCATE(lat_e, source = Ddata_read1) CALL read_from_gridfile(id_nc, 'float', 'angle_e') ALLOCATE(angle_e, source = Ddata_read1) CALL free_data_read ! free buffers after reading all data from grid file edge_num = SIZE(le_de) primal_num = SIZE(Ai) dual_num = SIZE(Av) max_primal_deg = SIZE(primal_edge,1) max_dual_deg = SIZE(dual_edge,1) max_trisk_deg = SIZE(trisk,1) ! now post-process some of the data we just read ALLOCATE(ep_e(edge_num,3)) DO ij=1,edge_num clon = COS(lon_e(ij)) slon = SIN(lon_e(ij)) clat = COS(lat_e(ij)) slat = SIN(lat_e(ij)) ! x,y,z = clat*clon, clat*slon, slat elon = (/ -slon, clon, 0. /) elat = (/ -slat*clon, -slat*slon, clat /) ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat END DO DEALLOCATE(angle_e) END SUBROUTINE read_dump_partition SUBROUTINE init_grid_type USE grid_param, ONLY : grid_type, grid_unst, grid_ico USE getin_mod, ONLY : getin CHARACTER(len=255) :: grid_type_var grid_type_var='icosahedral' CALL getin("grid_type",grid_type_var) SELECT CASE(grid_type_var) CASE('unstructured') grid_type = grid_unst ! is_omp_level_master=.TRUE. ! omp_level_size=1 CALL read_dump_partition IF (is_mpi_master) PRINT *,'Using unstructured grid type' CASE DEFAULT grid_type = grid_ico IF (is_mpi_master) PRINT *,'Using default grid type' END SELECT END SUBROUTINE init_grid_type END MODULE init_unstructured_mod