Changeset 883

Timestamp:

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

Author:

dubos

Message:

devel/unstructured : XIOS output

Location:

codes/icosagcm/devel/src

Files:

• base/field.f90 (modified) (8 diffs)
• base/math_const.f90 (modified) (1 diff)
• output/output_field.f90 (modified) (1 diff)
• output/write_field.f90 (modified) (1 diff)
• unstructured/init_unstructured.f90 (modified) (6 diffs)

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

@@ -1 +1 @@
 MODULE field_mod
   USE genmod
+  USE grid_param
+
   IMPLICIT NONE
   
@@ -44 +46 @@
   END INTERFACE
 
-  PRIVATE :: allocate_field_, deallocate_field_
+  ! This module is PUBLIC, we do not want to propagate symbols from modules USEd at module level
+  PRIVATE :: allocate_field_, deallocate_field_, grid_type, grid_unst, grid_ico
 
 CONTAINS
@@ -150 +153 @@
   SUBROUTINE allocate_field_XvalY(field, dom)
     USE domain_mod, ONLY : t_domain
-    USE grid_param
     TYPE(t_field) :: field
     TYPE(t_domain), INTENT(IN) :: dom
@@ -184 +186 @@
        
     CASE(grid_unst)
-       PRINT *, 'Allocating field ', field%name
        SELECT CASE(field%field_type)
        CASE(field_T)
@@ -194 +195 @@
        END SELECT
 
-       PRINT *, 'Allocating field ', field%name
-       PRINT *, '          with ij_size = ', ij_size
        
        IF (field%ndim==4) THEN
+          PRINT *, 'Allocating field ', TRIM(field%name), ' with shape = ', dim1, ij_size, dim2
           IF (data_type==type_integer) ALLOCATE(field%ival4d(dim1,ij_size,dim2))
           IF (data_type==type_real)    ALLOCATE(field%rval4d(dim1,ij_size,dim2))
           IF (data_type==type_logical) ALLOCATE(field%lval4d(dim1,ij_size,dim2))
        ELSE IF (field%ndim==3) THEN
+          PRINT *, 'Allocating field ', TRIM(field%name), ' with shape = ', dim1, ij_size
           IF (data_type==type_integer) ALLOCATE(field%ival3d(dim1,ij_size))
           IF (data_type==type_real)    ALLOCATE(field%rval3d(dim1,ij_size))
           IF (data_type==type_logical) ALLOCATE(field%lval3d(dim1,ij_size))
        ELSE IF (field%ndim==2) THEN
+          PRINT *, 'Allocating field ', TRIM(field%name), ' with shape = ', ij_size
           IF (data_type==type_integer) ALLOCATE(field%ival2d(ij_size))
           IF (data_type==type_real)    ALLOCATE(field%rval2d(ij_size))
@@ -362 +364 @@
     END IF
     field_pt=>field%rval2d
+    IF(grid_type == grid_unst) THEN
+       PRINT *, 'getval_r2d for ' // TRIM(field%name) // ' : ', SHAPE(field_pt)
+    END IF
   END SUBROUTINE  getval_r2d
 
@@ -375 +380 @@
     END IF
     field_pt=>field%rval3d
+    IF(grid_type == grid_unst) THEN
+       PRINT *, 'getval_r3d for ' // TRIM(field%name) // ' : ', SHAPE(field_pt)
+    END IF
   END SUBROUTINE  getval_r3d
 
@@ -387 +395 @@
     END IF
     field_pt=>field%rval4d
+    IF(grid_type == grid_unst) THEN
+       PRINT *, 'getval_r4d for ' // TRIM(field%name) // ' : ', SHAPE(field_pt)
+    END IF
   END SUBROUTINE  getval_r4d  
 

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

@@ -2 +2 @@
   USE PREC
   
-  REAL(rstd),PARAMETER :: Pi=acos(-1._rstd)
+  REAL(rstd), PARAMETER :: Pi=acos(-1._rstd), radian_to_degree=180._rstd/Pi
   COMPLEX(cstd), PARAMETER :: Imag=(0,1)
   

codes/icosagcm/devel/src/output/output_field.f90

@@ -36 +36 @@
     
     IF (xios_output) THEN
-       IF(grid_type == grid_ico) CALL xios_init_write_field
+       CALL xios_init_write_field
     ENDIF
   END SUBROUTINE output_field_init

codes/icosagcm/devel/src/output/write_field.f90

@@ -175 +175 @@
           DO n=1, cells(ind)%ncell
              ij=cells(ind)%ij(n)
-             PRINT *, 'write_2d :', ind, n, n_begin+n, ij 
              field_val2d(n_begin+n) = field(ind)%rval2d(ij)
           END DO

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