Merge pull request NOAA-EMC#197 from bensonr/emc_io_fixes

Emc io fixes

driver/fvGFS/atmosphere.F90

@@ -324,10 +324,12 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
    logical :: dycore_only  = .false.
    logical :: debug        = .false.
    logical :: sync         = .false.
+   logical :: ignore_rst_cksum = .false.
    integer, parameter     :: maxhr = 4096
    real, dimension(maxhr) :: fdiag = 0.
    real                   :: fhmax=384.0, fhmaxhf=120.0, fhout=3.0, fhouthf=1.0,avg_max_length=3600.
-   namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, avg_max_length
+   namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, &
+                            & avg_max_length, ignore_rst_cksum
    ! *DH 20210326
 
    !For regional
@@ -357,7 +359,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 #ifdef MOVING_NEST
    call fv_tracker_init(size(Atm))
    if (mygrid .eq. 2) call allocate_tracker(mygrid, Atm(mygrid)%bd%isc, Atm(mygrid)%bd%iec, Atm(mygrid)%bd%jsc, Atm(mygrid)%bd%jec)
-#endif   
+#endif
 
    Atm(mygrid)%Time_init = Time_init
 
@@ -450,6 +452,15 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 !--- allocate pref
    allocate(pref(npz+1,2), dum1d(npz+1))
 
+   ! DH* 20210326
+   ! First, read atmos_model_nml namelist section - this is a workaround to avoid
+   ! unnecessary additional changes to the input namelists, in anticipation of the
+   ! implementation of a generic interface for GFDL and CCPP fast physics soon
+   read(input_nml_file, nml=atmos_model_nml, iostat=io)
+   ierr = check_nml_error(io, 'atmos_model_nml')
+   !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'"
+   ! *DH 20210326
+
    call fv_restart(Atm(mygrid)%domain, Atm, seconds, days, cold_start, Atm(mygrid)%gridstruct%grid_type, mygrid)
 
    fv_time = Time
@@ -491,15 +502,6 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 
    ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics)
 
-   ! DH* 20210326
-   ! First, read atmos_model_nml namelist section - this is a workaround to avoid
-   ! unnecessary additional changes to the input namelists, in anticipation of the
-   ! implementation of a generic interface for GFDL and CCPP fast physics soon
-   read(input_nml_file, nml=atmos_model_nml, iostat=io)
-   ierr = check_nml_error(io, 'atmos_model_nml')
-   !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'"
-   ! *DH 20210326
-
    ! For fast physics running over the entire domain, block
    ! and thread number are not used; set to safe values
    cdata%blk_no = 1
@@ -954,10 +956,10 @@ end subroutine get_nth_domain_info
 !! the "domain2d" variable associated with the coupling grid and the
 !! decomposition for the current cubed-sphere tile.
 !>@detail Coupling is done using the mass/temperature grid with no halos.
- subroutine atmosphere_domain ( fv_domain, layout, regional, nested, &
+ subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, nested, &
                                 moving_nest_parent, is_moving_nest, &
                                 ngrids_atmos, mygrid_atmos, pelist )
-   type(domain2d), intent(out) :: fv_domain
+   type(domain2d), intent(out) :: fv_domain, rd_domain
    integer, intent(out) :: layout(2)
    logical, intent(out) :: regional
    logical, intent(out) :: nested
@@ -970,6 +972,7 @@ subroutine atmosphere_domain ( fv_domain, layout, regional, nested, &
    integer :: n
 
    fv_domain = Atm(mygrid)%domain_for_coupler
+   rd_domain = Atm(mygrid)%domain_for_read
    layout(1:2) =  Atm(mygrid)%layout(1:2)
    regional = Atm(mygrid)%flagstruct%regional
    nested = ngrids > 1

model/fv_arrays.F90

@@ -516,6 +516,7 @@ module fv_arrays_mod
 !-----------------------------------------------------------------------------------------------
 
    logical :: reset_eta = .false.
+   logical :: ignore_rst_cksum = .false. !< enfore (.false.) or override (.true.) data integrity restart checksums
    real    :: p_fac = 0.05  !< Safety factor for minimum nonhydrostatic pressures, which
                             !< will be limited so the full pressure is no less than p_fac
                             !< times the hydrostatic pressure. This is only of concern in mid-top
@@ -1303,6 +1304,7 @@ module fv_arrays_mod
 #if defined(SPMD)
 
      type(domain2D) :: domain_for_coupler !< domain used in coupled model with halo = 1.
+     type(domain2D) :: domain_for_read    !< domain used for reads to increase performance when io_layout=(1,1)
 
      !global tile and tile_of_mosaic only have a meaning for the CURRENT pe
      integer :: num_contact, npes_per_tile, global_tile, tile_of_mosaic, npes_this_grid

model/fv_control.F90

@@ -198,7 +198,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
      real,                intent(in)    :: dt_atmos
      integer,             intent(OUT)   :: this_grid
      logical, allocatable, intent(OUT) :: grids_on_this_pe(:)
-     character(len=32), optional,      intent(in)    :: nml_filename_in ! alternate nml 
+     character(len=32), optional,      intent(in)    :: nml_filename_in ! alternate nml
      logical, optional,                intent(in)    :: skip_nml_read_in ! use previously loaded nml
 
      integer, intent(INOUT) :: p_split
@@ -291,6 +291,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
      real(kind=R_GRID) , pointer :: target_lon
 
      logical , pointer :: reset_eta
+     logical , pointer :: ignore_rst_cksum
      real    , pointer :: p_fac
      real    , pointer :: a_imp
      integer , pointer :: n_split
@@ -676,7 +677,8 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
           Atm(this_grid)%flagstruct%grid_type,Atm(this_grid)%neststruct%nested, &
           Atm(this_grid)%layout,Atm(this_grid)%io_layout,Atm(this_grid)%bd,Atm(this_grid)%tile_of_mosaic, &
           Atm(this_grid)%gridstruct%square_domain,Atm(this_grid)%npes_per_tile,Atm(this_grid)%domain, &
-          Atm(this_grid)%domain_for_coupler,Atm(this_grid)%num_contact,Atm(this_grid)%pelist)
+          Atm(this_grid)%domain_for_coupler,Atm(this_grid)%domain_for_read,Atm(this_grid)%num_contact, &
+          Atm(this_grid)%pelist)
      call broadcast_domains(Atm,Atm(this_grid)%pelist,size(Atm(this_grid)%pelist))
      do n=1,ngrids
         tile_id = mpp_get_tile_id(Atm(n)%domain)
@@ -857,6 +859,7 @@ subroutine set_namelist_pointers(Atm)
        regional_bcs_from_gsi         => Atm%flagstruct%regional_bcs_from_gsi
        write_restart_with_bcs        => Atm%flagstruct%write_restart_with_bcs
        reset_eta                     => Atm%flagstruct%reset_eta
+       ignore_rst_cksum              => Atm%flagstruct%ignore_rst_cksum
        p_fac                         => Atm%flagstruct%p_fac
        a_imp                         => Atm%flagstruct%a_imp
        n_split                       => Atm%flagstruct%n_split
@@ -1075,7 +1078,7 @@ subroutine read_namelist_fv_core_nml(Atm)
             write_coarse_restart_files,&
             write_coarse_diagnostics,&
             write_only_coarse_intermediate_restarts, &
-            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst
+            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum
 
 
        ! Read FVCORE namelist

tools/external_ic.F90

@@ -206,10 +206,9 @@ module external_ic_mod
 
 contains
 
-   subroutine get_external_ic( Atm, fv_domain, cold_start )
+   subroutine get_external_ic( Atm, cold_start )
 
       type(fv_atmos_type), intent(inout), target :: Atm
-      type(domain2d),      intent(inout) :: fv_domain
       logical, intent(IN) :: cold_start
       real:: alpha = 0.
       real rdg
@@ -260,14 +259,14 @@ subroutine get_external_ic( Atm, fv_domain, cold_start )
          enddo
       enddo
 
-      call mpp_update_domains( f0, fv_domain )
+      call mpp_update_domains( f0, Atm%domain )
       if ( Atm%gridstruct%cubed_sphere .and. (.not. Atm%gridstruct%bounded_domain))then
          call fill_corners(f0, Atm%npx, Atm%npy, YDir)
       endif
 
 ! Read in cubed_sphere terrain
       if ( Atm%flagstruct%mountain ) then
-           call get_cubed_sphere_terrain(Atm, fv_domain)
+           call get_cubed_sphere_terrain(Atm)
       else
          if (.not. Atm%neststruct%nested) Atm%phis = 0. !TODO: Not sure about this line --- lmh 30 may 18
       endif
@@ -276,32 +275,32 @@ subroutine get_external_ic( Atm, fv_domain, cold_start )
       if ( Atm%flagstruct%ncep_ic ) then
            nq = 1
                              call timing_on('NCEP_IC')
-           call get_ncep_ic( Atm, fv_domain, nq )
+           call get_ncep_ic( Atm, nq )
                              call timing_off('NCEP_IC')
 #ifdef FV_TRACERS
            if (.not. cold_start) then
-              call fv_io_read_tracers( fv_domain, Atm )
+              call fv_io_read_tracers( Atm )
               if(is_master()) write(*,*) 'All tracers except sphum replaced by FV IC'
            endif
 #endif
       elseif ( Atm%flagstruct%nggps_ic ) then
                              call timing_on('NGGPS_IC')
-           call get_nggps_ic( Atm, fv_domain )
+           call get_nggps_ic( Atm )
                              call timing_off('NGGPS_IC')
       elseif ( Atm%flagstruct%hrrrv3_ic ) then
                              call timing_on('HRRR_IC')
-           call get_hrrr_ic( Atm, fv_domain )
+           call get_hrrr_ic( Atm )
                              call timing_off('HRRR_IC')
       elseif ( Atm%flagstruct%ecmwf_ic ) then
            if( is_master() ) write(*,*) 'Calling get_ecmwf_ic'
                              call timing_on('ECMWF_IC')
-           call get_ecmwf_ic( Atm, fv_domain )
+           call get_ecmwf_ic( Atm )
                              call timing_off('ECMWF_IC')
       else
 ! The following is to read in legacy lat-lon FV core restart file
 !  is Atm%q defined in all cases?
            nq = size(Atm%q,4)
-           call get_fv_ic( Atm, fv_domain, nq )
+           call get_fv_ic( Atm, nq )
       endif
 
       call prt_maxmin('PS', Atm%ps, is, ie, js, je, ng, 1, 0.01)
@@ -368,9 +367,8 @@ end subroutine get_external_ic
 
 
 !------------------------------------------------------------------
-  subroutine get_cubed_sphere_terrain( Atm, fv_domain )
+  subroutine get_cubed_sphere_terrain( Atm )
     type(fv_atmos_type), intent(inout), target :: Atm
-    type(domain2d),      intent(inout) :: fv_domain
     type(FmsNetcdfDomainFile_t) :: Fv_core
     integer :: tile_id(1)
     character(len=64)    :: fname
@@ -393,13 +391,13 @@ subroutine get_cubed_sphere_terrain( Atm, fv_domain )
     jed = Atm%bd%jed
     ng  = Atm%bd%ng
 
-    tile_id = mpp_get_tile_id( fv_domain )
+    tile_id = mpp_get_tile_id( Atm%domain )
 
     fname = 'INPUT/fv_core.res.nc'
     call mpp_error(NOTE, 'external_ic: looking for '//fname)
 
 
-       if( open_file(Fv_core, fname, "read", fv_domain, is_restart=.true.) ) then
+       if( open_file(Fv_core, fname, "read", Atm%domain_for_read, is_restart=.true.) ) then
          call read_data(Fv_core, 'phis', Atm%phis(is:ie,js:je))
          call close_file(Fv_core)
        else
@@ -428,7 +426,7 @@ end subroutine get_cubed_sphere_terrain
 !>@brief The subroutine 'get_nggps_ic' reads in data after it has been preprocessed with
 !!    NCEP/EMC orography maker and 'global_chgres', and has been horiztontally
 !! interpolated to the current cubed-sphere grid
-  subroutine get_nggps_ic (Atm, fv_domain)
+  subroutine get_nggps_ic (Atm)
 
 !>variables read in from 'gfs_ctrl.nc'
 !>       VCOORD  -  level information
@@ -456,7 +454,6 @@ subroutine get_nggps_ic (Atm, fv_domain)
 #endif
 
     type(fv_atmos_type), intent(inout) :: Atm
-    type(domain2d),      intent(inout) :: fv_domain
 ! local:
     real, dimension(:), allocatable:: ak, bk
     real, dimension(:,:), allocatable:: wk2, ps, oro_g
@@ -585,7 +582,7 @@ subroutine get_nggps_ic (Atm, fv_domain)
 
     !--- read in surface temperature (k) and land-frac
     ! surface skin temperature
-   if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+   if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
       naxis_dims = get_variable_num_dimensions(SFC_restart, 'tsea')
       allocate (dim_names_alloc(naxis_dims))
       call get_variable_dimension_names(SFC_restart, 'tsea', dim_names_alloc)
@@ -605,7 +602,7 @@ subroutine get_nggps_ic (Atm, fv_domain)
     dim_names_2d(2) = "lon"
 
     ! terrain surface height -- (needs to be transformed into phis = zs*grav)
-    if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+    if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
       call register_axis(ORO_restart, "lat", "y")
       call register_axis(ORO_restart, "lon", "x")
       if (filtered_terrain) then
@@ -915,7 +912,7 @@ subroutine read_gfs_ic()
         dim_names_3d4(1) = "levp"
 
         ! surface pressure (Pa)
-        if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+        if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
           call register_axis(GFS_restart, "lat", "y")
           call register_axis(GFS_restart, "lon", "x")
           call register_axis(GFS_restart, "lonp", "x", domain_position=east)
@@ -945,7 +942,7 @@ subroutine read_gfs_ic()
           do nt = 1, ntracers
 
               q(:,:,:,nt) = -999.99
-             
+
             call get_tracer_names(MODEL_ATMOS, nt, tracer_name)
             call register_restart_field(GFS_restart, trim(tracer_name), q(:,:,:,nt), dim_names_3d3, is_optional=.true.)
           enddo
@@ -964,7 +961,7 @@ subroutine read_gfs_ic()
   end subroutine get_nggps_ic
 !------------------------------------------------------------------
 !------------------------------------------------------------------
-  subroutine get_hrrr_ic (Atm, fv_domain)
+  subroutine get_hrrr_ic (Atm)
 !    read in data after it has been preprocessed with
 !    NCEP/EMC orography maker
 !
@@ -990,7 +987,6 @@ subroutine get_hrrr_ic (Atm, fv_domain)
 
 
       type(fv_atmos_type), intent(inout) :: Atm
-      type(domain2d),      intent(inout) :: fv_domain
 ! local:
       real, dimension(:), allocatable:: ak, bk
       real, dimension(:,:), allocatable:: wk2, ps, oro_g
@@ -1104,7 +1100,7 @@ subroutine get_hrrr_ic (Atm, fv_domain)
 
 !--- read in surface temperature (k) and land-frac
         ! surface skin temperature
-       if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+       if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
           call get_variable_dimension_names(SFC_restart, 'tsea', dim_names_2d)
           call register_axis(SFC_restart, dim_names_2d(2), "y")
           call register_axis(SFC_restart, dim_names_2d(1), "x")
@@ -1121,7 +1117,7 @@ subroutine get_hrrr_ic (Atm, fv_domain)
         dim_names_2d(2) = "lon"
 
         ! terrain surface height -- (needs to be transformed into phis = zs*grav)
-        if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+        if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
           call register_axis(ORO_restart, "lat", "y")
           call register_axis(ORO_restart, "lon", "x")
           if (filtered_terrain) then
@@ -1164,7 +1160,7 @@ subroutine get_hrrr_ic (Atm, fv_domain)
         dim_names_3d4(1) = "levp"
 
         ! edge pressure (Pa)
-        if( open_file(HRRR_restart, fn_hrr_ics, "read", Atm%domain,is_restart=.true., dont_add_res_to_filename=.true.) ) then
+        if( open_file(HRRR_restart, fn_hrr_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
           call register_axis(HRRR_restart, "lat", "y")
           call register_axis(HRRR_restart, "lon", "x")
           call register_axis(HRRR_restart, "lonp", "x", domain_position=east)
@@ -1360,9 +1356,8 @@ end subroutine get_hrrr_ic
 !------------------------------------------------------------------
 !------------------------------------------------------------------
 !>@brief The subroutine 'get_ncep_ic' reads in the specified NCEP analysis or reanalysis dataset
-  subroutine get_ncep_ic( Atm, fv_domain, nq )
+  subroutine get_ncep_ic( Atm, nq )
       type(fv_atmos_type), intent(inout) :: Atm
-      type(domain2d),      intent(inout) :: fv_domain
       integer, intent(in):: nq
 ! local:
 #ifdef HIWPP_ETA
@@ -1818,7 +1813,7 @@ end subroutine get_ncep_ic
 !>@brief The subroutine 'get_ecmwf_ic' reads in initial conditions from ECMWF analyses
 !! (EXPERIMENTAL: contact Jan-Huey Chen jan-huey.chen@noaa.gov for support)
 !>@authors Jan-Huey Chen, Xi Chen, Shian-Jiann Lin
-  subroutine get_ecmwf_ic( Atm, fv_domain )
+  subroutine get_ecmwf_ic( Atm )
 
 #ifdef __PGI
       use GFS_restart, only : GFS_restart_type
@@ -1827,7 +1822,6 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
 #endif
 
       type(fv_atmos_type), intent(inout) :: Atm
-      type(domain2d),      intent(inout) :: fv_domain
 ! local:
       real :: ak_ec(138), bk_ec(138)
       data ak_ec/ 0.000000,     2.000365,     3.102241,     4.666084,     6.827977,     9.746966, &
@@ -2062,7 +2056,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
       dim_names_3d4(1) = "levp"
 
 !! Read in model terrain from oro_data.tile?.nc
-      if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+      if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
         call register_axis(ORO_restart, "lat", "y")
         call register_axis(ORO_restart, "lon", "x")
         if (filtered_terrain) then
@@ -2088,7 +2082,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
       allocate (ps_gfs(is:ie,js:je))
       allocate (zh_gfs(is:ie,js:je,levp_gfs+1))
 
-      if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then
+      if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then
         call register_axis(GFS_restart, "lat", "y")
         call register_axis(GFS_restart, "lon", "x")
         call register_axis(GFS_restart, "levp", size(zh_gfs,3))
@@ -2622,9 +2616,8 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
   end subroutine get_ecmwf_ic
 !------------------------------------------------------------------
 !------------------------------------------------------------------
-  subroutine get_fv_ic( Atm, fv_domain, nq )
+  subroutine get_fv_ic( Atm, nq )
       type(fv_atmos_type), intent(inout) :: Atm
-      type(domain2d),      intent(inout) :: fv_domain
       integer, intent(in):: nq
 
       type(FmsNetcdfFile_t) :: Latlon_dyn, Latlon_tra