Use the same real kind for axis variables in restart files as real ki…

…nd used for data variables. (#697) * Write netcdf axis variables using the same real kind as data variables
NOAA-EMC · Sep 26, 2023 · 3b4423c · 3b4423c
1 parent bbc5bf8
commit 3b4423c
Show file tree

Hide file tree

Showing 7 changed files with 45 additions and 35 deletions.
diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere
diff --git a/io/fv3atm_clm_lake_io.F90 b/io/fv3atm_clm_lake_io.F90
@@ -14,7 +14,7 @@ module fv3atm_clm_lake_io
        register_restart_field, write_data, &
        register_variable_attribute, register_field, get_dimension_size
   use fv3atm_common_io,   only: create_2d_field_and_add_to_bundle, &
-       create_3d_field_and_add_to_bundle
+       create_3d_field_and_add_to_bundle, axis_type
 
   implicit none
 
@@ -179,16 +179,16 @@ subroutine clm_lake_write_axes(clm_lake, Model, Sfc_restart)
     type(GFS_control_type),      intent(in) :: Model
     type(FmsNetcdfDomainFile_t) :: Sfc_restart
     integer :: i
-    call register_field(Sfc_restart, 'levlake_clm_lake', 'double', (/'levlake_clm_lake'/))
+    call register_field(Sfc_restart, 'levlake_clm_lake', axis_type, (/'levlake_clm_lake'/))
     call register_variable_attribute(Sfc_restart, 'levlake_clm_lake', 'cartesian_axis' ,'Z', str_len=1)
 
-    call register_field(Sfc_restart, 'levsoil_clm_lake', 'double', (/'levsoil_clm_lake'/))
+    call register_field(Sfc_restart, 'levsoil_clm_lake', axis_type, (/'levsoil_clm_lake'/))
     call register_variable_attribute(Sfc_restart, 'levsoil_clm_lake', 'cartesian_axis' ,'Z', str_len=1)
 
-    call register_field(Sfc_restart, 'levsnowsoil_clm_lake', 'double', (/'levsnowsoil_clm_lake'/))
+    call register_field(Sfc_restart, 'levsnowsoil_clm_lake', axis_type, (/'levsnowsoil_clm_lake'/))
     call register_variable_attribute(Sfc_restart, 'levsnowsoil_clm_lake', 'cartesian_axis' ,'Z', str_len=1)
 
-    call register_field(Sfc_restart, 'levsnowsoil1_clm_lake', 'double', (/'levsnowsoil1_clm_lake'/))
+    call register_field(Sfc_restart, 'levsnowsoil1_clm_lake', axis_type, (/'levsnowsoil1_clm_lake'/))
     call register_variable_attribute(Sfc_restart, 'levsnowsoil1_clm_lake', 'cartesian_axis' ,'Z', str_len=1)
 
     call write_data(Sfc_restart, 'levlake_clm_lake', clm_lake%levlake_clm_lake)

diff --git a/io/fv3atm_common_io.F90 b/io/fv3atm_common_io.F90
@@ -31,6 +31,12 @@ module fv3atm_common_io
 
   public :: get_nx_ny_from_atm
 
+#ifdef CCPP_32BIT
+  character(len=5), parameter, public :: axis_type = 'float'
+#else
+  character(len=6), parameter, public :: axis_type = 'double'
+#endif
+
   !>\defgroup fv3atm_common_io FV3ATM Common I/O Utilities Module
   !> @{
 

diff --git a/io/fv3atm_restart_io.F90 b/io/fv3atm_restart_io.F90
@@ -17,7 +17,7 @@ module fv3atm_restart_io_mod
                                 get_global_io_domain_indices, get_dimension_size
   use mpp_domains_mod,    only: domain2d
   use fv3atm_common_io,   only: create_2d_field_and_add_to_bundle, &
-       create_3d_field_and_add_to_bundle, copy_from_gfs_data
+       create_3d_field_and_add_to_bundle, copy_from_gfs_data, axis_type
   use fv3atm_sfc_io
   use fv3atm_rrfs_sd_io
   use fv3atm_clm_lake_io
@@ -913,23 +913,23 @@ subroutine phys_restart_write (GFS_Restart, Atm_block, Model, fv_domain, timesta
     amiopen=open_file(Phy_restart, trim(infile), 'overwrite', domain=fv_domain, is_restart=.true., dont_add_res_to_filename=.true.)
     if( amiopen ) then
       call register_axis(Phy_restart, 'xaxis_1', 'X')
-      call register_field(Phy_restart, 'xaxis_1', 'double', (/'xaxis_1'/))
+      call register_field(Phy_restart, 'xaxis_1', axis_type, (/'xaxis_1'/))
       call register_variable_attribute(Phy_restart, 'xaxis_1', 'cartesian_axis', 'X', str_len=1)
       call get_global_io_domain_indices(Phy_restart, 'xaxis_1', is, ie, indices=buffer)
       call write_data(Phy_restart, "xaxis_1", buffer)
       deallocate(buffer)
       call get_dimension_size(Phy_restart, 'xaxis_1', xaxis_1_chunk)
 
       call register_axis(Phy_restart, 'yaxis_1', 'Y')
-      call register_field(Phy_restart, 'yaxis_1', 'double', (/'yaxis_1'/))
+      call register_field(Phy_restart, 'yaxis_1', axis_type, (/'yaxis_1'/))
       call register_variable_attribute(Phy_restart, 'yaxis_1', 'cartesian_axis', 'Y', str_len=1)
       call get_global_io_domain_indices(Phy_restart, 'yaxis_1', is, ie, indices=buffer)
       call write_data(Phy_restart, "yaxis_1", buffer)
       deallocate(buffer)
       call get_dimension_size(Phy_restart, 'yaxis_1', yaxis_1_chunk)
 
       call register_axis(Phy_restart, 'zaxis_1', phy%npz)
-      call register_field(Phy_restart, 'zaxis_1', 'double', (/'zaxis_1'/))
+      call register_field(Phy_restart, 'zaxis_1', axis_type, (/'zaxis_1'/))
       call register_variable_attribute(Phy_restart, 'zaxis_1', 'cartesian_axis', 'Z', str_len=1)
       allocate( buffer(phy%npz) )
       do i=1, phy%npz
@@ -939,7 +939,7 @@ subroutine phys_restart_write (GFS_Restart, Atm_block, Model, fv_domain, timesta
       deallocate(buffer)
 
       call register_axis(Phy_restart, 'Time', unlimited)
-      call register_field(Phy_restart, 'Time', 'double', (/'Time'/))
+      call register_field(Phy_restart, 'Time', axis_type, (/'Time'/))
       call register_variable_attribute(Phy_restart, 'Time', 'cartesian_axis', 'T', str_len=1)
       call write_data(Phy_restart, "Time", 1)
     else

diff --git a/io/fv3atm_rrfs_sd_io.F90 b/io/fv3atm_rrfs_sd_io.F90
@@ -10,7 +10,7 @@ module fv3atm_rrfs_sd_io
                                 get_dimension_size
   use GFS_typedefs,       only: GFS_sfcprop_type, GFS_control_type, kind_phys
   use fv3atm_common_io,   only: get_nx_ny_from_atm, create_2d_field_and_add_to_bundle, &
-                                create_3d_field_and_add_to_bundle
+                                create_3d_field_and_add_to_bundle, axis_type
 
   implicit none
 
@@ -114,7 +114,7 @@ subroutine rrfs_sd_state_write_axis(data,Model,Sfc_restart)
     type(FmsNetcdfDomainFile_t) :: Sfc_restart
     type(GFS_control_type),      intent(in) :: Model
 
-    call register_field(Sfc_restart, 'fire_aux_data_levels', 'double', (/'fire_aux_data_levels'/))
+    call register_field(Sfc_restart, 'fire_aux_data_levels', axis_type, (/'fire_aux_data_levels'/))
     call register_variable_attribute(Sfc_restart, 'fire_aux_data_levels', 'cartesian_axis' ,'Z', str_len=1)
     call write_data(Sfc_restart, 'fire_aux_data_levels', data%fire_aux_data_levels)
   end subroutine rrfs_sd_state_write_axis

diff --git a/io/fv3atm_sfc_io.F90 b/io/fv3atm_sfc_io.F90
@@ -11,7 +11,7 @@ module fv3atm_sfc_io
                                 register_variable_attribute, register_field, &
                                 get_global_io_domain_indices, variable_exists, &
                                 get_dimension_size
-  use fv3atm_common_io,   only: GFS_Data_transfer, &
+  use fv3atm_common_io,   only: GFS_Data_transfer, axis_type, &
        create_2d_field_and_add_to_bundle, create_3d_field_and_add_to_bundle
   use GFS_typedefs,       only: GFS_sfcprop_type, GFS_control_type, kind_phys
   use mpp_mod,            only: mpp_error,  NOTE
@@ -309,19 +309,19 @@ subroutine Sfc_io_write_axes(sfc, Model, Sfc_restart)
     integer :: i, is, ie
     logical :: mand
 
-    call register_field(Sfc_restart, 'xaxis_1', 'double', (/'xaxis_1'/))
+    call register_field(Sfc_restart, 'xaxis_1', axis_type, (/'xaxis_1'/))
     call register_variable_attribute(Sfc_restart, 'xaxis_1', 'cartesian_axis', 'X', str_len=1)
     call get_global_io_domain_indices(Sfc_restart, 'xaxis_1', is, ie, indices=buffer)
     call write_data(Sfc_restart, "xaxis_1", buffer)
     deallocate(buffer)
 
-    call register_field(Sfc_restart, 'yaxis_1', 'double', (/'yaxis_1'/))
+    call register_field(Sfc_restart, 'yaxis_1', axis_type, (/'yaxis_1'/))
     call register_variable_attribute(Sfc_restart, 'yaxis_1', 'cartesian_axis', 'Y', str_len=1)
     call get_global_io_domain_indices(Sfc_restart, 'yaxis_1', is, ie, indices=buffer)
     call write_data(Sfc_restart, "yaxis_1", buffer)
     deallocate(buffer)
 
-    call register_field(Sfc_restart, 'zaxis_1', 'double', (/'zaxis_1'/))
+    call register_field(Sfc_restart, 'zaxis_1', axis_type, (/'zaxis_1'/))
     call register_variable_attribute(Sfc_restart, 'zaxis_1', 'cartesian_axis', 'Z', str_len=1)
     allocate( buffer(Model%kice) )
     do i=1, Model%kice
@@ -331,7 +331,7 @@ subroutine Sfc_io_write_axes(sfc, Model, Sfc_restart)
     deallocate(buffer)
 
     if (Model%lsm == Model%lsm_noah .or. Model%lsm == Model%lsm_noahmp) then
-      call register_field(Sfc_restart, 'zaxis_2', 'double', (/'zaxis_2'/))
+      call register_field(Sfc_restart, 'zaxis_2', axis_type, (/'zaxis_2'/))
       call register_variable_attribute(Sfc_restart, 'zaxis_2', 'cartesian_axis', 'Z', str_len=1)
       allocate( buffer(Model%lsoil) )
       do i=1, Model%lsoil
@@ -342,7 +342,7 @@ subroutine Sfc_io_write_axes(sfc, Model, Sfc_restart)
     endif
 
     if(Model%lsm == Model%lsm_noahmp) then
-      call register_field(Sfc_restart, 'zaxis_3', 'double', (/'zaxis_3'/))
+      call register_field(Sfc_restart, 'zaxis_3', axis_type, (/'zaxis_3'/))
       call register_variable_attribute(Sfc_restart, 'zaxis_3', 'cartesian_axis', 'Z', str_len=1)
       allocate(buffer(3))
       do i=1, 3
@@ -351,7 +351,7 @@ subroutine Sfc_io_write_axes(sfc, Model, Sfc_restart)
       call write_data(Sfc_restart, 'zaxis_3', buffer)
       deallocate(buffer)
 
-      call register_field(Sfc_restart, 'zaxis_4', 'double', (/'zaxis_4'/))
+      call register_field(Sfc_restart, 'zaxis_4', axis_type, (/'zaxis_4'/))
       call register_variable_attribute(Sfc_restart, 'zaxis_4', 'cartesian_axis' ,'Z', str_len=1)
       allocate(buffer(7))
       do i=1, 7
@@ -360,7 +360,7 @@ subroutine Sfc_io_write_axes(sfc, Model, Sfc_restart)
       call write_data(Sfc_restart, 'zaxis_4', buffer)
       deallocate(buffer)
     end if
-    call register_field(Sfc_restart, 'Time', 'double', (/'Time'/))
+    call register_field(Sfc_restart, 'Time', axis_type, (/'Time'/))
     call register_variable_attribute(Sfc_restart, 'Time', 'cartesian_axis', 'T', str_len=1)
     call write_data( Sfc_restart, 'Time', 1)
   end subroutine Sfc_io_write_axes

diff --git a/io/module_write_restart_netcdf.F90 b/io/module_write_restart_netcdf.F90
@@ -74,7 +74,7 @@ subroutine write_restart_netcdf(wrtfb, filename, &
     integer :: ncerr,ierr
     integer :: ncid
     integer :: oldMode
-    integer :: dimid
+    integer :: dimid, dimtype
     integer :: im_dimid, im_p1_dimid, jm_dimid, jm_p1_dimid, time_dimid
     integer :: im_varid, im_p1_varid, jm_varid, jm_p1_varid, time_varid
     integer, dimension(:), allocatable :: dimids_2d, dimids_3d
@@ -188,6 +188,15 @@ subroutine write_restart_netcdf(wrtfb, filename, &
           deallocate(maxIndexPTile)
           deallocate(deToTileMap)
           deallocate(localDeToDeMap)
+
+          if (typekind == ESMF_TYPEKIND_R4) then
+             dimtype = NF90_FLOAT
+          else if (typekind == ESMF_TYPEKIND_R8) then
+             dimtype = NF90_DOUBLE
+          else
+             if (mype==0) write(0,*)'Unsupported typekind ', typekind
+             call ESMF_Finalize(endflag=ESMF_END_ABORT)
+          end if
        end if
 
        if (fieldDimCount > gridDimCount) then
@@ -236,29 +245,29 @@ subroutine write_restart_netcdf(wrtfb, filename, &
        if ( .not.is_restart_core ) then
 
           ncerr = nf90_def_dim(ncid, "xaxis_1", im, im_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "xaxis_1", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "xaxis_1", dimtype, im_dimid, im_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, im_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr)
 
           ncerr = nf90_def_dim(ncid, "yaxis_1", jm, jm_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "yaxis_1", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "yaxis_1", dimtype, jm_dimid, jm_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, jm_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr)
 
        else
 
           ncerr = nf90_def_dim(ncid, "xaxis_1", im, im_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "xaxis_1", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "xaxis_1", dimtype, im_dimid, im_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, im_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr)
 
           ncerr = nf90_def_dim(ncid, "xaxis_2", im+1, im_p1_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "xaxis_2", NF90_DOUBLE, im_p1_dimid, im_p1_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "xaxis_2", dimtype, im_p1_dimid, im_p1_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, im_p1_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr)
 
           ncerr = nf90_def_dim(ncid, "yaxis_1", jm+1, jm_p1_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "yaxis_1", NF90_DOUBLE, jm_p1_dimid, jm_p1_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "yaxis_1", dimtype, jm_p1_dimid, jm_p1_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, jm_p1_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr)
 
           ncerr = nf90_def_dim(ncid, "yaxis_2", jm, jm_dimid); NC_ERR_STOP(ncerr)
-          ncerr = nf90_def_var(ncid, "yaxis_2", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr)
+          ncerr = nf90_def_var(ncid, "yaxis_2", dimtype, jm_dimid, jm_varid); NC_ERR_STOP(ncerr)
           ncerr = nf90_put_att(ncid, jm_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr)
 
        end if
@@ -291,7 +300,7 @@ subroutine write_restart_netcdf(wrtfb, filename, &
 
        ncerr = nf90_def_dim(ncid, "Time", NF90_UNLIMITED, time_dimid); NC_ERR_STOP(ncerr)
        ! ncerr = nf90_def_dim(ncid, "Time", 1, time_dimid); NC_ERR_STOP(ncerr)
-       ncerr = nf90_def_var(ncid, "Time", NF90_DOUBLE, time_dimid, time_varid); NC_ERR_STOP(ncerr)
+       ncerr = nf90_def_var(ncid, "Time", dimtype, time_dimid, time_varid); NC_ERR_STOP(ncerr)
        if (par) then
           ncerr = nf90_var_par_access(ncid, time_varid, NF90_COLLECTIVE); NC_ERR_STOP(ncerr)
        end if
@@ -565,12 +574,7 @@ subroutine write_out_ungridded_dim_atts_from_field(field, dimLabel, dimid, rc)
         ncerr = nf90_def_dim(ncid, trim(dimLabel), valueCount, dimid=dimid); NC_ERR_STOP(ncerr); NC_ERR_STOP(ncerr)
       endif
       if( typekind == ESMF_TYPEKIND_R4 ) then
-        !!! FIXME Use NF90_DOUBLE as axis type, even though axis data are float
-        !!!       This is needed to make phy/sfc restart files identical to FMS
-        !!!       restart files which always defines all axis as double
-
-        ! ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_FLOAT, dimids=(/dimid/), varid=varid); NC_ERR_STOP(ncerr)
-        ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_DOUBLE, dimids=(/dimid/), varid=varid); NC_ERR_STOP(ncerr)
+        ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_FLOAT, dimids=(/dimid/), varid=varid); NC_ERR_STOP(ncerr)
         ncerr = nf90_put_att(ncid, varid, trim(axis_attr_name), "Z"); NC_ERR_STOP(ncerr)
         ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
         ncerr = nf90_put_var(ncid, varid, values=valueListr4); NC_ERR_STOP(ncerr)
+9 −3		driver/fvGFS/fv_ufs_restart_io.F90
+12 −6		tools/fv_io.F90