Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

various io fixes #192

Merged
merged 4 commits into from
May 12, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion driver/GFDL/atmosphere.F90
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Surf_diff, Grid_box)
if ( Atm(mygrid)%flagstruct%nudge ) then
call fv_ada_nudge_init( Time, Atm(mygrid)%atmos_axes, npz, zvir, Atm(mygrid)%ak, Atm(mygrid)%bk, &
Atm(mygrid)%ts, Atm(mygrid)%phis, Atm(mygrid)%gridstruct, Atm(mygrid)%ks, Atm(mygrid)%npx, &
Atm(mygrid)%neststruct, Atm(mygrid)%bd, Atm(mygrid)%domain)
Atm(mygrid)%neststruct, Atm(mygrid)%bd, Atm(mygrid)%domain_for_read)
call mpp_error(NOTE, 'ADA nudging is active')
endif
#else
Expand Down
5 changes: 3 additions & 2 deletions driver/SHiELD/atmosphere.F90
Original file line number Diff line number Diff line change
Expand Up @@ -713,15 +713,16 @@ subroutine set_atmosphere_pelist ()
end subroutine set_atmosphere_pelist


subroutine atmosphere_domain ( fv_domain, layout, regional, bounded_domain )
type(domain2d), intent(out) :: fv_domain
subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, bounded_domain )
type(domain2d), intent(out) :: fv_domain, rd_domain
integer, intent(out) :: layout(2)
logical, intent(out) :: regional
logical, intent(out) :: bounded_domain
! returns the domain2d variable associated with the coupling grid
! note: coupling is done using the mass/temperature grid with no halos

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
Expand Down
2 changes: 2 additions & 0 deletions model/fv_arrays.F90
Original file line number Diff line number Diff line change
Expand Up @@ -511,6 +511,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
Expand Down Expand Up @@ -1280,6 +1281,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
Expand Down
7 changes: 5 additions & 2 deletions model/fv_control.F90
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,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
Expand Down Expand Up @@ -550,7 +551,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)
Expand Down Expand Up @@ -728,6 +730,7 @@ subroutine set_namelist_pointers(Atm)
write_restart_with_bcs => Atm%flagstruct%write_restart_with_bcs
regional_bcs_from_gsi => Atm%flagstruct%regional_bcs_from_gsi
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
Expand Down Expand Up @@ -940,7 +943,7 @@ subroutine read_namelist_fv_core_nml(Atm)
w_limiter, write_coarse_restart_files, write_coarse_diagnostics,&
write_only_coarse_intermediate_restarts, &
write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, &
pass_full_omega_to_physics_in_non_hydrostatic_mode
pass_full_omega_to_physics_in_non_hydrostatic_mode, ignore_rst_cksum


! Read FVCORE namelist
Expand Down
57 changes: 25 additions & 32 deletions tools/external_ic.F90
Original file line number Diff line number Diff line change
Expand Up @@ -87,10 +87,9 @@ module external_ic_mod

contains

subroutine get_external_ic( Atm, fv_domain, cold_start, icdir )
subroutine get_external_ic( Atm, cold_start, icdir )

type(fv_atmos_type), intent(inout), target :: Atm
type(domain2d), intent(inout) :: fv_domain
logical, intent(IN) :: cold_start
character(len=*), intent(in), optional :: icdir
real:: alpha = 0.
Expand Down Expand Up @@ -139,14 +138,14 @@ subroutine get_external_ic( Atm, fv_domain, cold_start, icdir )
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
Expand All @@ -155,32 +154,32 @@ subroutine get_external_ic( Atm, fv_domain, cold_start, icdir )
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)
Expand Down Expand Up @@ -219,9 +218,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
Expand All @@ -244,13 +242,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
Expand All @@ -276,7 +274,7 @@ subroutine get_cubed_sphere_terrain( Atm, fv_domain )
end subroutine get_cubed_sphere_terrain


subroutine get_nggps_ic (Atm, fv_domain)
subroutine get_nggps_ic (Atm)
! read in data after it has been preprocessed with
! NCEP/EMC orography maker and global_chgres
! and has been horiztontally interpolated to the
Expand All @@ -301,7 +299,6 @@ subroutine get_nggps_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
Expand Down Expand Up @@ -424,7 +421,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)
Expand All @@ -444,7 +441,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
Expand Down Expand Up @@ -757,7 +754,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)
Expand Down Expand Up @@ -803,7 +800,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
!
Expand All @@ -824,7 +821,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
Expand Down Expand Up @@ -939,7 +935,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")
Expand All @@ -956,7 +952,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
Expand Down Expand Up @@ -999,7 +995,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)
Expand Down Expand Up @@ -1194,9 +1190,8 @@ subroutine get_hrrr_ic (Atm, fv_domain)
end subroutine get_hrrr_ic
!------------------------------------------------------------------
!------------------------------------------------------------------
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
Expand Down Expand Up @@ -1652,9 +1647,8 @@ subroutine get_ncep_ic( Atm, fv_domain, nq )
end subroutine get_ncep_ic
!------------------------------------------------------------------
!------------------------------------------------------------------
subroutine get_ecmwf_ic( Atm, fv_domain )
subroutine get_ecmwf_ic( Atm )
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, &
Expand Down Expand Up @@ -1871,7 +1865,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
Expand All @@ -1891,7 +1885,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, "lev", size(o3mr_gfs,3))
Expand Down Expand Up @@ -2397,9 +2391,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
Expand Down
Loading