Changeset 883 for codes/icosagcm/devel/src/unstructured/init_unstructured.f90

Timestamp:

06/11/19 17:38:17 (5 years ago)

Author:

dubos

Message:

devel/unstructured : XIOS output

File:

• codes/icosagcm/devel/src/unstructured/init_unstructured.f90 (modified) (6 diffs)

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

@@ -1 +1 @@
 MODULE init_unstructured_mod
+  USE prec, ONLY : rstd
+  USE math_const, ONLY : radian_to_degree
   USE mpipara, ONLY : is_mpi_master
   USE data_unstructured_mod
@@ -95 +97 @@
   END SUBROUTINE read_from_gridfile
 
-
   SUBROUTINE open_local_mesh_file
     USE netcdf_mod
@@ -153 +154 @@
   END SUBROUTINE read_field
 
+
+  SUBROUTINE setup_cellset(set, lon, lat)
+    USE domain_mod, ONLY : t_cellset
+    TYPE(t_cellset) :: set
+    REAL(rstd) :: lon(:), lat(:)
+    INTEGER :: ij
+    set%ncell = SIZE(lon)
+    ALLOCATE(set%lon, SOURCE=lon*radian_to_degree)
+    ALLOCATE(set%lat, SOURCE=lat*radian_to_degree)
+    ALLOCATE(set%ij, SOURCE = (/(ij,ij=1,set%ncell)/))
+    ALLOCATE(set%ind_glo, SOURCE = set%ij(:)-1)
+  END SUBROUTINE setup_cellset
+
+  SUBROUTINE read_local_mesh_bounds(dom)
+    USE domain_mod, ONLY : t_domain
+    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e
+    TYPE(t_domain), INTENT(IN) :: dom
+    REAL(rstd) :: bounds_e(2,edge_num)
+    REAL(rstd), ALLOCATABLE :: lon_v(:), lat_v(:)
+    INTEGER :: ij, vup, vdown
+
+    ! read dual cell positions, not stored in geom
+    CALL read_from_gridfile(id_nc, 'float', 'lon_v')
+    ALLOCATE(lon_v, SOURCE=Ddata_read1)
+    CALL read_from_gridfile(id_nc, 'float', 'lat_v')
+    ALLOCATE(lat_v, SOURCE=Ddata_read1)
+
+    CALL setup_cellset(dom%primal_own, lon_i, lat_i)
+    CALL setup_cellset(dom%dual_own, lon_v, lat_v)
+    CALL setup_cellset(dom%edge_own, lon_e, lat_e)
+
+    PRINT *, 'read_local_mesh :', primal_num, SHAPE(dom%primal_own%ij)
+
+    ! read primal cell bounds
+    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lon')
+    ALLOCATE(dom%primal_own%bnds_lon,     source=Ddata_read2*radian_to_degree)
+    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lat')
+    ALLOCATE(dom%primal_own%bnds_lat,     source=Ddata_read2*radian_to_degree)
+    ! read dual cell bounds
+    CALL read_from_gridfile(id_nc, 'float', 'dual_bounds_lon')
+    ALLOCATE(dom%dual_own%bnds_lon,     source=Ddata_read2*radian_to_degree)
+    CALL read_from_gridfile(id_nc, 'float', 'dual_bounds_lat')
+    ALLOCATE(dom%dual_own%bnds_lat,     source=Ddata_read2*radian_to_degree)
+
+    CALL collect_lonlat_e(lat_v)
+    ALLOCATE(dom%edge_own%bnds_lat, SOURCE=bounds_e)
+    CALL collect_lonlat_e(lon_v)
+    ALLOCATE(dom%edge_own%bnds_lon, SOURCE=bounds_e)
+
+    CONTAINS
+      SUBROUTINE collect_lonlat_e(lonlat)
+        REAL(rstd) :: lonlat(dual_num)
+        DO ij=1, edge_num
+           vdown = down(ij)
+           vup = up(ij)
+           bounds_e(1,ij) = lonlat(vdown)*radian_to_degree
+           bounds_e(2,ij) = lonlat(vup)*radian_to_degree
+        END DO
+      END SUBROUTINE collect_lonlat_e
+
+  END SUBROUTINE read_local_mesh_bounds
+
   SUBROUTINE read_local_mesh
     USE field_mod
+    USE domain_mod, ONLY : swap_needed, domain, domain_glo
     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
+    USE netcdf_mod, ONLY : nf90_close
     IMPLICIT NONE
-    INTEGER :: ij
+    INTEGER :: ij, ierr
     INTEGER, ALLOCATABLE :: cell_ij(:)
     REAL(rstd), ALLOCATABLE :: angle_e(:)
@@ -178 +242 @@
     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
 
@@ -234 +284 @@
     ALLOCATE(wee, source = Ddata_read2)
 
+    ! read cell centers and bounds
+    CALL read_local_mesh_bounds(domain(1))
+    CALL read_local_mesh_bounds(domain_glo(1))
+
     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)
+    ierr = nf90_close(id_nc)
+
+    ! now process some of the data we just read
+
     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))
@@ -254 +305 @@
        elon = (/ -slon, clon, 0. /)
        elat = (/ -slat*clon, -slat*slon, clat /)
-       ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
+       ep_e(:,ij) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
     END DO
     
     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