Changeset 879

Timestamp:

06/04/19 18:03:54 (5 years ago)

Author:

dubos

Message:

devel : introduced derived type to store cell bounds

Location:

codes/icosagcm/devel/src

Files:

• base/field.f90 (modified) (1 diff)
• parallel/domain.f90 (modified) (14 diffs)
• sphere/geometry.f90 (modified) (3 diffs)
• unstructured/init_unstructured.f90 (modified) (6 diffs)

codes/icosagcm/devel/src/base/field.f90

@@ -133 +133 @@
         field(ind)%ndim=3
         field(ind)%dim3=dim1
+        field(ind)%dim4=1
       ELSE
         field(ind)%ndim=2
+        field(ind)%dim3=1
+        field(ind)%dim4=1
       END IF
     

codes/icosagcm/devel/src/parallel/domain.f90

@@ -1 +1 @@
 MODULE domain_mod
   USE domain_param
+  IMPLICIT NONE
+  SAVE
+
+  TYPE t_cellset
+     ! derived type containing arrays for lat,lon of center and bounds of cells of primal or dual mesh
+     ! this information is needed to write CF-compliant NetCDF files
+     ! and by XIOS
+     INTEGER :: ncell
+     INTEGER, ALLOCATABLE :: ij(:), & ! ij(n) = local index of cell n
+          ind_glo(:), &  ! ind_glo(n) = global index of cell n
+          sgn(:)         ! sign to apply before writing to disk (velocities)
+     REAL, ALLOCATABLE :: lon(:), lat(:), bnds_lon(:,:), bnds_lat(:,:)
+  END type t_cellset
 
   TYPE t_domain
@@ -60 +73 @@
     INTEGER :: u_pos(6)
     INTEGER :: z_pos(6)
-     
+
+    TYPE(t_cellset), POINTER :: primal_own, dual_own, edge_own, &
+         primal_all, dual_all, edge_all ! halo included, used for debug by writefield
   END TYPE t_domain
 
-  INTEGER,SAVE :: ndomain
-  INTEGER,SAVE :: ndomain_glo
+  TYPE t_mesh
+     INTEGER :: ndomain, primal_num, edge_num, dual_num
+     TYPE(t_domain), ALLOCATABLE :: domain(:)
+     TYPE(t_cellset), ALLOCATABLE :: primal_own(:), dual_own(:), edge_own(:), &
+       primal_all(:), dual_all(:), edge_all(:)
+  END type t_mesh
+
+  INTEGER :: ndomain
+  INTEGER :: ndomain_glo
 
   LOGICAL :: swap_needed=.TRUE. ! .FALSE. if a thread always computes on the same domain
   !$OMP THREADPRIVATE(swap_needed)
 
-  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:)
-  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:)
-
-  INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:)  ! size : ndomain_glo : mpi rank assigned to a domain 
-  INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) ! size : ndomain_glo : corresponding local indice for a global domain indice
-  INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) ! size : ndomain : corresponding global indice for a local domain indice
-
-  LOGICAL, ALLOCATABLE, SAVE :: assigned_domain(:) ! size : ndomain : is an omp task is assigned to this domain
+  TYPE(t_mesh), TARGET :: mesh_loc, mesh_glo
+  TYPE(t_domain), POINTER :: domain(:), domain_glo(:)
+
+  INTEGER,ALLOCATABLE :: domglo_rank(:)  ! size : ndomain_glo : mpi rank assigned to a domain 
+  INTEGER,ALLOCATABLE :: domglo_loc_ind(:) ! size : ndomain_glo : corresponding local indice for a global domain indice
+  INTEGER,ALLOCATABLE :: domloc_glo_ind(:) ! size : ndomain : corresponding global indice for a local domain indice
+
+  LOGICAL, ALLOCATABLE :: assigned_domain(:) ! size : ndomain : is an omp task is assigned to this domain
 !$OMP THREADPRIVATE(assigned_domain)
     
 CONTAINS
+
+  SUBROUTINE allocate_mesh(ndom, mesh)
+    INTEGER, INTENT(IN)       :: ndom
+    TYPE(t_mesh), TARGET, INTENT(OUT) :: mesh
+    INTEGER :: ind
+    mesh%ndomain = ndom
+    ALLOCATE(mesh%domain(ndom))
+    ALLOCATE(mesh%primal_own(ndom))
+    ALLOCATE(mesh%primal_all(ndom))
+    ALLOCATE(mesh%dual_own(ndom))
+    ALLOCATE(mesh%dual_all(ndom))
+    ALLOCATE(mesh%edge_own(ndom))
+    ALLOCATE(mesh%edge_all(ndom))
+    DO ind=1, ndom
+       mesh%domain(ind)%primal_own => mesh%primal_own(ind)
+       mesh%domain(ind)%primal_all => mesh%primal_all(ind)
+       mesh%domain(ind)%dual_own => mesh%dual_own(ind)
+       mesh%domain(ind)%dual_all => mesh%dual_all(ind)
+       mesh%domain(ind)%edge_own => mesh%edge_own(ind)
+       mesh%domain(ind)%edge_all => mesh%edge_all(ind)
+    END DO
+  END SUBROUTINE allocate_mesh
 
   SUBROUTINE create_domain
@@ -85 +129 @@
   USE mpipara
   USE ioipsl
-  IMPLICIT NONE
   INTEGER :: ind,nf,ni,nj,i,j
   INTEGER :: quotient, rest
@@ -93 +136 @@
  
     ndomain_glo=nsplit_i*nsplit_j*nb_face
-    ALLOCATE(domain_glo(ndomain_glo))
+    CALL allocate_mesh(ndomain_glo, mesh_glo)
+    domain_glo => mesh_glo%domain
     ALLOCATE(domglo_rank(ndomain_glo))
     ALLOCATE(domglo_loc_ind(ndomain_glo))
@@ -178 +222 @@
   
   SUBROUTINE copy_domain(d1,d2)
-  IMPLICIT NONE
   INTEGER :: face
   TYPE(t_domain),TARGET,INTENT(IN) :: d1
@@ -245 +288 @@
   SUBROUTINE assign_cell
   USE metric
-  IMPLICIT NONE
     INTEGER :: ind_d,ind,ind2,e,v
     INTEGER :: nf,nf2
@@ -365 +407 @@
   SUBROUTINE compute_boundary
   USE spherical_geom_mod
-  IMPLICIT NONE
     INTEGER :: ind_d
     INTEGER :: i,j,k
@@ -388 +429 @@
   SUBROUTINE set_neighbour_indice
   USE metric
-  IMPLICIT NONE
     INTEGER :: ind_d
     TYPE(t_domain),POINTER :: d 
@@ -444 +484 @@
   USE mpipara
   USE grid_param
-  IMPLICIT NONE
     INTEGER :: nb_domain(0:mpi_size-1)
     INTEGER :: rank, ind,ind_glo
@@ -456 +495 @@
     
     ndomain=nb_domain(mpi_rank)
-    ALLOCATE(domain(ndomain))
+    
+    ! ALLOCATE(domain(ndomain))
+    CALL allocate_mesh(ndomain, mesh_loc)
+    domain => mesh_loc%domain
+
     ALLOCATE(domloc_glo_ind(ndomain))
     
@@ -588 +631 @@
   USE omp_para
   USE mpipara
-  IMPLICIT NONE
     INTEGER :: nb_domain
     INTEGER :: rank, ind, i
@@ -629 +671 @@
   SUBROUTINE compute_domain
     USE grid_param, ONLY : grid_type, grid_unst, grid_ico
-    IMPLICIT NONE
+    INTEGER :: ind
     SELECT CASE(grid_type)
     CASE(grid_unst)
        ! FIXME temporary, sequential hack
        ndomain_glo=1
-       ALLOCATE(domain_glo(ndomain_glo))
+       CALL allocate_mesh(ndomain_glo, mesh_glo)
+       domain_glo => mesh_glo%domain
        ALLOCATE(domglo_rank(ndomain_glo))
        ALLOCATE(domglo_loc_ind(ndomain_glo))
@@ -641 +684 @@
 
        ndomain=1
-       ALLOCATE(domain(ndomain))
+       CALL allocate_mesh(ndomain, mesh_loc)
+       domain => mesh_loc%domain
        ALLOCATE(domloc_glo_ind(ndomain))
        ALLOCATE(assigned_domain(ndomain))
@@ -658 +702 @@
        CALL assign_domain
     END SELECT
+
+!    ALLOCATE(primal_cells(ndomain))
+!    ALLOCATE(dual_cells(ndomain))
+!    ALLOCATE(primal_cells_single(ndomain))
+!    ALLOCATE(dual_cells_single(ndomain))
+!    DO ind=1, ndomain
+!       domain(ind)%primal_cells => primal_cells(ind)
+!       domain(ind)%dual_cells => dual_cells(ind)
+!       domain(ind)%primal_cells_single => primal_cells_single(ind)
+!       domain(ind)%dual_cells_single => dual_cells_single(ind)
+!    END DO
+
   END SUBROUTINE compute_domain
           

codes/icosagcm/devel/src/sphere/geometry.f90

@@ -104 +104 @@
   USE field_mod
   IMPLICIT NONE
-  
+  INTEGER, PARAMETER :: nvertex=6    ! FIXME unstructured
+
     CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
 
     CALL allocate_field(geom%xyz_i,field_t,type_real,3)
-    CALL allocate_field(geom%lon_i,field_t,type_real)
-    CALL allocate_field(geom%lat_i,field_t,type_real)
+    CALL allocate_field(geom%lon_i,field_t,type_real, name='lon_i')
+    CALL allocate_field(geom%lat_i,field_t,type_real, name='lat_i')
     CALL allocate_field(geom%elon_i,field_t,type_real,3)
     CALL allocate_field(geom%elat_i,field_t,type_real,3)
@@ -115 +116 @@
 
     CALL allocate_field(geom%xyz_e,field_u,type_real,3)
-    CALL allocate_field(geom%lon_e,field_u,type_real)
-    CALL allocate_field(geom%lat_e,field_u,type_real)
+    CALL allocate_field(geom%lon_e,field_u,type_real, name='lon_e')
+    CALL allocate_field(geom%lat_e,field_u,type_real, name='lat_e')
     CALL allocate_field(geom%elon_e,field_u,type_real,3)
     CALL allocate_field(geom%elat_e,field_u,type_real,3)
@@ -123 +124 @@
 
     CALL allocate_field(geom%xyz_v,field_z,type_real,3)
-    CALL allocate_field(geom%de,field_u,type_real)
-    CALL allocate_field(geom%le,field_u,type_real)
-    CALL allocate_field(geom%le_de,field_u,type_real)
+    CALL allocate_field(geom%de,field_u,type_real, name='de')
+    CALL allocate_field(geom%le,field_u,type_real, name='le')
+    CALL allocate_field(geom%le_de,field_u,type_real, name='le_de')
     CALL allocate_field(geom%bi,field_t,type_real)
-    CALL allocate_field(geom%Av,field_z,type_real)
-    CALL allocate_field(geom%S1,field_t,type_real,6)
-    CALL allocate_field(geom%S2,field_t,type_real,6)
-    CALL allocate_field(geom%Riv,field_t,type_real,6)
-    CALL allocate_field(geom%Riv2,field_t,type_real,6)
-    CALL allocate_field(geom%ne,field_t,type_integer,6)
-    CALL allocate_field(geom%Wee,field_u,type_real,5,2)
+    CALL allocate_field(geom%Av,field_z,type_real, name='Av')
+    CALL allocate_field(geom%S1,field_t,type_real,nvertex)
+    CALL allocate_field(geom%S2,field_t,type_real,nvertex)
+    CALL allocate_field(geom%Riv,field_t,type_real,nvertex)
+    CALL allocate_field(geom%Riv2,field_t,type_real,nvertex)
+    CALL allocate_field(geom%ne,field_t,type_integer,nvertex)
+    CALL allocate_field(geom%Wee,field_u,type_real,5,2) ! FIXME unstructured
     CALL allocate_field(geom%bi,field_t,type_real)
-    CALL allocate_field(geom%fv,field_z,type_real)
+    CALL allocate_field(geom%fv,field_z,type_real, name='fv')
 
   END SUBROUTINE allocate_geometry

codes/icosagcm/devel/src/unstructured/init_unstructured.f90

@@ -2 +2 @@
   USE mpipara, ONLY : is_mpi_master
   USE data_unstructured_mod
-  USE geometry, ONLY : de
+  USE geometry, ONLY : swap_geometry
   IMPLICIT NONE
   SAVE
@@ -128 +128 @@
 
 
+  SUBROUTINE read_field(id_nc, field)
+    USE field_mod, ONLY : t_field, type_integer, type_real
+    INTEGER :: id_nc
+    TYPE(t_field), POINTER :: field(:), fld
+    fld=>field(1)
+    SELECT CASE(fld%data_type)
+    CASE(type_integer)
+       CALL read_from_gridfile(id_nc, 'integer', TRIM(fld%name))
+       SELECT CASE(field(1)%ndim)
+       CASE(2)
+          fld%ival2d = Idata_read1
+       CASE(3)
+          fld%ival3d = Idata_read2
+       END SELECT
+    CASE(type_real)
+       CALL read_from_gridfile(id_nc, 'float', TRIM(fld%name))
+       SELECT CASE(field(1)%ndim)
+       CASE(2)
+          fld%rval2d = Ddata_read1
+       CASE(3)
+          fld%rval3d = Ddata_read2
+       END SELECT
+    END SELECT
+  END SUBROUTINE read_field
+
   SUBROUTINE read_local_mesh
     USE field_mod
-    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e
+    USE geometry, ONLY : geom, lon_i, lat_i, lon_e, lat_e, ep_e
+    USE domain_mod, ONLY : domain, domain_glo, t_domain, swap_needed
     IMPLICIT NONE
     INTEGER :: ij
+    INTEGER, ALLOCATABLE :: cell_ij(:)
     REAL(rstd), ALLOCATABLE :: angle_e(:)
     REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat
          elon(3), elat(3) ! lon/lat unit vectors
     
-    !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
-
+    ! if possible data is read into fields allocated in geom(1) using read_field
+    CALL read_field(id_nc, geom%Ai)
+    CALL read_field(id_nc, geom%Av)
+    CALL read_field(id_nc, geom%fv)
+    CALL read_field(id_nc, geom%le_de)
+    CALL read_field(id_nc, geom%le)
+    CALL read_field(id_nc, geom%de)
+    CALL read_field(id_nc, geom%lon_i)
+    CALL read_field(id_nc, geom%lat_i)
+    CALL read_field(id_nc, geom%lon_e)
+    CALL read_field(id_nc, geom%lat_e)
+
+    swap_needed = .TRUE.
+    CALL swap_geometry(1)
+
+    ! cell centers and bounds are read into domain(1)%primal_own
+    domain(1)%primal_own%ncell = primal_num
+    domain_glo(1)%primal_own%ncell = primal_num
+    ALLOCATE(domain(1)%primal_own%lon,     source=lon_i)
+    ALLOCATE(domain_glo(1)%primal_own%lon, source=lon_i)
+    ALLOCATE(domain(1)%primal_own%lat,     source=lat_i)
+    ALLOCATE(domain_glo(1)%primal_own%lat, source=lat_i)
+
+    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lon')
+    ALLOCATE(domain(1)%primal_own%bnds_lon,     source=Ddata_read2)
+    ALLOCATE(domain_glo(1)%primal_own%bnds_lon, source=Ddata_read2)
+    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lat')
+    ALLOCATE(domain(1)%primal_own%bnds_lat,     source=Ddata_read2)
+    ALLOCATE(domain_glo(1)%primal_own%bnds_lat, source=Ddata_read2)
+    PRINT *, 'read_local_mesh : primal_num =', primal_num, domain_glo(1)%primal_own%ncell
+
+    ! other fields are defined only in the case of an unstructured mesh
+    ! those fields are read into simple arrays
     CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
     ALLOCATE(primal_deg, source = Idata_read1)
@@ -148 +203 @@
     CALL read_from_gridfile(id_nc, 'integer', 'primal_ne')
     ALLOCATE(primal_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', 'dual_deg')
     ALLOCATE(dual_deg, source = Idata_read1)
@@ -154 +212 @@
     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', 'dual_vertex') 
+    ALLOCATE(dual_vertex, source = Idata_read2)
+
     CALL read_from_gridfile(id_nc, 'integer', 'left')
     ALLOCATE(left, source = Idata_read1)
@@ -164 +223 @@
     CALL read_from_gridfile(id_nc, 'integer', 'down')
     ALLOCATE(down, source = Idata_read1)
+    CALL read_from_gridfile(id_nc, 'float', 'angle_e')
+    ALLOCATE(angle_e, source = Ddata_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', 'le')
-    ALLOCATE(le, source = Ddata_read1)
-    CALL read_from_gridfile(id_nc, 'float', 'de')
-    ALLOCATE(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)
+
+!    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)
@@ -221 +258 @@
     
     DEALLOCATE(angle_e)
+
+    ALLOCATE(domain(1)%primal_own%ij(primal_num))
+    domain(1)%primal_own%ij(:) = (/(ij,ij=1,primal_num)/)
+    PRINT *, 'read_local_mesh :', SHAPE(domain), primal_num, &
+         SHAPE(domain(1)%primal_own%ij)
+
+    ALLOCATE(cell_ij(dual_num))
+    cell_ij(:) = (/(ij,ij=1,dual_num)/)
+    ALLOCATE(domain(1)%dual_own%ij, source=cell_ij)
+    DEALLOCATE(cell_ij)
+
+    CALL swap_geometry(1)
+    swap_needed = .FALSE.
   END SUBROUTINE read_local_mesh