Skip to content

Commit

Permalink
Merge pull request #71 from NOAA-GFDL/main_to_dev
Browse files Browse the repository at this point in the history 
Main to dev
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA authored Feb 21, 2022
2 parents 149073f + 712ff9e commit 2e72b88
Show file tree
Hide file tree
Showing 2 changed files with 99 additions and 29 deletions.
107 changes: 89 additions & 18 deletions src/framework/MOM_io.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -892,9 +892,17 @@ end subroutine read_variable_1d_int

!> Read a 2d array from a netCDF input file and save to a variable.
!!
!! Start and nread ranks may exceed var, but must match the rank of the
!! variable in the netCDF file. This allows for reading slices of larger
!! arrays.
!! Start and nread lenths may exceed var rank. This allows for reading slices
!! of larger arrays.
!!
!! Previous versions of the model required a time axis on IO fields. This
!! constraint was dropped in later versions. As a result, versions both with
!! and without a time axis now exist. In order to support all such versions,
!! we use a reshaped version of start and nread in order to read the variable
!! as it exists in the file.
!!
!! Certain constraints are still applied to start and nread in order to ensure
!! that varname is a valid 2d array, or contains valid 2d slices.
!!
!! I/O occurs only on the root PE, and data is broadcast to other ranks.
!! Due to potentially large memory communication and storage, this subroutine
Expand All @@ -908,11 +916,40 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in)
integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file.
!! If absent, the file is opened and closed within this routine.

integer :: ncid, varid, ndims, rc
character(len=*), parameter :: hdr = "read_variable_2d"
integer :: ncid, varid
integer :: field_ndims, dim_len
integer, allocatable :: field_dimids(:), field_shape(:)
integer, allocatable :: field_start(:), field_nread(:)
integer :: i, rc
character(len=*), parameter :: hdr = "read_variable_2d: "
character(len=128) :: msg
logical :: size_mismatch

! Validate shape of start and nread
if (present(start)) then
if (size(start) < 2) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
// " start must have at least two dimensions.")
endif

if (present(nread)) then
if (size(nread) < 2) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
// " nread must have at least two dimensions.")

if (any(nread(3:) > 1)) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
// " nread may only read a single level in higher dimensions.")
endif

! Since start and nread may be reshaped, we cannot rely on netCDF to ensure
! that their lengths are equivalent, and must do it here.
if (present(start) .and. present(nread)) then
if (size(start) /= size(nread)) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
// " start and nread must have the same length.")
endif

! Open and read `varname` from `filename`
if (is_root_pe()) then
if (present(ncid_in)) then
ncid = ncid_in
Expand All @@ -923,23 +960,57 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in)
call get_varid(varname, ncid, filename, varid, match_case=.false.)
if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
" in "//trim(filename))
! Verify that start(:) and nread(:) ranks match variable's dimension count
rc = nf90_inquire_variable(ncid, varid, ndims=ndims)

! Query for the dimensionality of the input field
rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims)
if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
" Difficulties reading "//trim(varname)//" from "//trim(filename))
": Difficulties reading "//trim(varname)//" from "//trim(filename))

! Confirm that field is at least 2d
if (field_ndims < 2) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) // " " // &
trim(varname) // " from " // trim(filename) // " is not a 2D field.")

! If start and nread are present, then reshape them to match field dims
if (present(start) .or. present(nread)) then
allocate(field_shape(field_ndims))
allocate(field_dimids(field_ndims))

rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids)
if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
": Difficulties reading "//trim(varname)//" from "//trim(filename))

do i = 1, field_ndims
rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len)
if (rc /= NF90_NOERR) &
call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
// ": Difficulties reading dimensions from " // trim(filename))
field_shape(i) = dim_len
enddo

size_mismatch = .false.
if (present(start)) size_mismatch = size_mismatch .or. size(start) /= ndims
if (present(nread)) size_mismatch = size_mismatch .or. size(nread) /= ndims
! Reshape start(:) and nreads(:) in case ranks differ
allocate(field_start(field_ndims))
field_start(:) = 1
if (present(start)) then
dim_len = min(size(start), size(field_start))
field_start(:dim_len) = start(:dim_len)
endif

allocate(field_nread(field_ndims))
field_nread(:2) = field_shape(:2)
field_nread(3:) = 1
if (present(nread)) field_shape(:2) = nread(:2)

if (size_mismatch) then
write (msg, '("'// hdr //': size(start) ", i0, " and/or size(nread) ", &
i0, " do not match ndims ", i0)') size(start), size(nread), ndims
call MOM_error(FATAL, trim(msg))
rc = nf90_get_var(ncid, varid, var, field_start, field_nread)

deallocate(field_start)
deallocate(field_nread)
deallocate(field_shape)
deallocate(field_dimids)
else
rc = nf90_get_var(ncid, varid, var)
endif
! NOTE: We could check additional information here (type, size, ...)

rc = nf90_get_var(ncid, varid, var, start, nread)
if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
" Difficulties reading "//trim(varname)//" from "//trim(filename))

Expand Down
21 changes: 10 additions & 11 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ module MOM_diabatic_driver
type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Control structure for a child module
type(bulkmixedlayer_CS) :: bulkmixedlayer !< Bulk mixed layer control struct
type(CVMix_conv_CS) :: CVMix_conv !< CVMix convection control struct
type(energetic_PBL_CS) :: energetic_PBL !< Energetic PBL control struct
type(energetic_PBL_CS) :: ePBL !< Energetic PBL control struct
type(entrain_diffusive_CS) :: entrain_diffusive !< Diffusive entrainment control struct
type(geothermal_CS) :: geothermal !< Geothermal control struct
type(int_tide_CS) :: int_tide !< Internal tide control struct
Expand Down Expand Up @@ -838,15 +838,15 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim

call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
CS%energetic_PBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)

if (associated(Hml)) then
call energetic_PBL_get_MLD(CS%energetic_PBL, Hml(:,:), G, US)
call energetic_PBL_get_MLD(CS%ePBL, Hml(:,:), G, US)
call pass_var(Hml, G%domain, halo=1)
! If visc%MLD exists, copy ePBL's MLD into it
if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
elseif (associated(visc%MLD)) then
call energetic_PBL_get_MLD(CS%energetic_PBL, visc%MLD, G, US)
call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US)
call pass_var(visc%MLD, G%domain, halo=1)
endif

Expand Down Expand Up @@ -1375,15 +1375,15 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,

call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
CS%energetic_PBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)

if (associated(Hml)) then
call energetic_PBL_get_MLD(CS%energetic_PBL, Hml(:,:), G, US)
call energetic_PBL_get_MLD(CS%ePBL, Hml(:,:), G, US)
call pass_var(Hml, G%domain, halo=1)
! If visc%MLD exists, copy ePBL's MLD into it
if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
elseif (associated(visc%MLD)) then
call energetic_PBL_get_MLD(CS%energetic_PBL, visc%MLD, G, US)
call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US)
call pass_var(visc%MLD, G%domain, halo=1)
endif

Expand Down Expand Up @@ -2589,14 +2589,13 @@ subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit,
if (present(opacity_CSp)) opacity_CSp => CS%opacity
if (present(optics_CSp)) optics_CSp => CS%optics
if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp
if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL
if (present(energetic_PBL_CSp) .and. CS%use_energetic_PBL) energetic_PBL_CSp => CS%ePBL
if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp

! Constants within diabatic_CS
if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit
if (present(minimum_forcing_depth)) minimum_forcing_depth = CS%minimum_forcing_depth
if (present(diabatic_halo)) diabatic_halo = CS%halo_TS_diff

end subroutine extract_diabatic_member

!> Routine called for adiabatic physics
Expand Down Expand Up @@ -3487,7 +3486,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
if (CS%use_bulkmixedlayer) &
call bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS%bulkmixedlayer)
if (CS%use_energetic_PBL) &
call energetic_PBL_init(Time, G, GV, US, param_file, diag, CS%energetic_PBL)
call energetic_PBL_init(Time, G, GV, US, param_file, diag, CS%ePBL)

call regularize_layers_init(Time, G, GV, param_file, diag, CS%regularize_layers)

Expand Down Expand Up @@ -3522,7 +3521,7 @@ subroutine diabatic_driver_end(CS)
call diapyc_energy_req_end(CS%diapyc_en_rec_CSp)

if (CS%use_energetic_PBL) &
call energetic_PBL_end(CS%energetic_PBL)
call energetic_PBL_end(CS%ePBL)

call diabatic_aux_end(CS%diabatic_aux_CSp)

Expand Down

0 comments on commit 2e72b88

Please sign in to comment.