module map_lnd2rof_irrig_mod !--------------------------------------------------------------------- ! ! Purpose: ! ! This module contains routines for mapping the irrigation field from the LND grid onto ! the ROF grid. ! ! These routines could go in prep_rof_mod, but are separated into their own module for ! the sake of (1) testability: this module has fewer dependencies than prep_rof_mod; ! and (2) symmetry with the lnd2glc and glc2lnd custom mapping routines, which also ! have their own modules. #include "shr_assert.h" use shr_kind_mod, only : r8 => shr_kind_r8 use mct_mod use seq_map_type_mod, only : seq_map use seq_map_mod, only : seq_map_map use shr_log_mod, only : errMsg => shr_log_errMsg implicit none private ! ------------------------------------------------------------------------ ! Public interfaces ! ------------------------------------------------------------------------ public :: map_lnd2rof_irrig ! map irrigation from lnd -> rof grid ! ------------------------------------------------------------------------ ! Private interfaces ! ------------------------------------------------------------------------ private :: map_rof2lnd_volr ! map volr from rof -> lnd grid character(len=*), parameter, private :: sourcefile = & __FILE__ contains subroutine map_lnd2rof_irrig(l2r_l, r2x_r, irrig_flux_field, & avwts_s, avwtsfld_s, mapper_Fl2r, mapper_Fr2l, l2r_r) !--------------------------------------------------------------- ! Description ! Do custom mapping for the irrigation flux, from land -> rof. ! ! The basic idea is that we want to pull irrigation out of ROF cells proportionally to ! the river volume (volr) in each cell. This is important in cases where the various ! ROF cells overlapping a CLM cell have very different volr: If we didn't do this ! volr-normalized remapping, we'd try to extract the same amount of water from each ! of the ROF cells, which would be more likely to have withdrawals exceeding ! available volr. ! ! (Both RTM and MOSART have code to handle excess withdrawals, by pulling the excess ! directly out of the ocean, but we'd like to avoid resorting to this as much as ! possible.) ! ! This mapping works by: ! ! (1) Normalizing the land's irrigation flux by volr ! ! (2) Mapping this volr-normalized flux to the rof grid ! ! (3) Converting the mapped, volr-normalized flux back to a normal ! (non-volr-normalized) flux on the rof grid. ! ! This assumes that the following fields are contained in the attribute vector ! arguments: ! ! - l2r_l: field given by irrig_flux_field (read) ! - l2r_r: field given by irrig_flux_field (set) ! - r2x_r: 'Flrr_volrmch' (read) ! ! Arguments type(mct_aVect) , intent(in) :: l2r_l ! lnd -> rof fields on the land grid type(mct_aVect) , intent(in) :: r2x_r ! rof -> cpl fields on the rof grid character(len=*) , intent(in) :: irrig_flux_field ! name of irrigation field to remap type(mct_aVect) , intent(in) :: avwts_s ! attr vect for source weighting character(len=*) , intent(in) :: avwtsfld_s ! field in avwts_s to use type(seq_map) , intent(inout) :: mapper_Fl2r ! flux mapper for mapping lnd -> rof type(seq_map) , intent(inout) :: mapper_Fr2l ! flux mapper for mapping rof -> lnd type(mct_aVect) , intent(inout) :: l2r_r ! lnd -> rof fields on the rof grid ! ! Local variables integer :: r, l integer :: lsize_l ! number of land points integer :: lsize_r ! number of rof points type(mct_avect) :: irrig_l_av ! temporary attribute vector holding irrigation fluxes on the land grid type(mct_avect) :: irrig_r_av ! temporary attribute vector holding irrigation fluxes on the rof grid ! The following need to be pointers to satisfy the MCT interface: real(r8), pointer :: volr_r(:) ! river volume on the rof grid real(r8), pointer :: volr_l(:) ! river volume on the land grid real(r8), pointer :: irrig_flux_l(:) ! irrigation flux on the land grid [kg m-2 s-1] real(r8), pointer :: irrig_flux_r(:) ! irrigation flux on the rof grid [kg m-2 s-1] real(r8), pointer :: irrig_normalized_l(:) ! irrigation normalized by volr, land grid real(r8), pointer :: irrig_normalized_r(:) ! irrigation normalized by volr, rof grid real(r8), pointer :: irrig_volr0_l(:) ! irrigation where volr <= 0, land grid real(r8), pointer :: irrig_volr0_r(:) ! irrigation where volr <= 0, rof grid character(len=*), parameter :: volr_field = 'Flrr_volrmch' character(len=*), parameter :: irrig_normalized_field = 'Flrl_irrig_normalized' character(len=*), parameter :: irrig_volr0_field = 'Flrl_irrig_volr0' character(len=*), parameter :: fields_to_remap = & irrig_normalized_field // ':' // irrig_volr0_field !--------------------------------------------------------------- ! ------------------------------------------------------------------------ ! Determine attribute vector sizes ! ------------------------------------------------------------------------ lsize_l = mct_aVect_lsize(l2r_l) lsize_r = mct_aVect_lsize(l2r_r) ! ------------------------------------------------------------------------ ! Extract the necessary fields from attribute vectors ! ------------------------------------------------------------------------ allocate(irrig_flux_l(lsize_l)) call mct_aVect_exportRattr(l2r_l, irrig_flux_field, irrig_flux_l) allocate(volr_r(lsize_r)) call mct_aVect_exportRattr(r2x_r, volr_field, volr_r) ! ------------------------------------------------------------------------ ! Adjust volr_r, and map it to the land grid ! ------------------------------------------------------------------------ ! Treat any rof point with volr < 0 as if it had volr = 0. Negative volr values can ! arise in RTM. This fix is needed to avoid mapping negative irrigation to those ! cells: while conservative, this would be unphysical (it would mean that irrigation ! actually adds water to those cells). do r = 1, lsize_r if (volr_r(r) < 0._r8) then volr_r(r) = 0._r8 end if end do allocate(volr_l(lsize_l)) call map_rof2lnd_volr(volr_r, mapper_Fr2l, volr_l) ! ------------------------------------------------------------------------ ! Determine irrigation normalized by volr ! ! In order to avoid possible divide by 0, as well as to handle non-sensical negative ! volr on the land grid, we divide the land's irrigation flux into two separate flux ! components: a component where we have positive volr on the land grid (put in ! irrig_normalized_l, which is mapped using volr-normalization) and a component where ! we have zero or negative volr on the land grid (put in irrig_volr0_l, which is ! mapped as a standard flux). We then remap both of these components to the rof grid, ! and then finally add the two components to determine the total irrigation flux on ! the rof grid. ! ------------------------------------------------------------------------ allocate(irrig_normalized_l(lsize_l)) allocate(irrig_volr0_l(lsize_l)) do l = 1, lsize_l if (volr_l(l) > 0._r8) then irrig_normalized_l(l) = irrig_flux_l(l) / volr_l(l) irrig_volr0_l(l) = 0._r8 else irrig_normalized_l(l) = 0._r8 irrig_volr0_l(l) = irrig_flux_l(l) end if end do ! ------------------------------------------------------------------------ ! Map irrigation ! ------------------------------------------------------------------------ call mct_aVect_init(irrig_l_av, rList = fields_to_remap, lsize = lsize_l) call mct_aVect_importRattr(irrig_l_av, irrig_normalized_field, irrig_normalized_l) call mct_aVect_importRattr(irrig_l_av, irrig_volr0_field, irrig_volr0_l) call mct_aVect_init(irrig_r_av, rList = fields_to_remap, lsize = lsize_r) ! This mapping uses the same options (such as avwts) as is used for mapping all other ! fields in prep_rof_calc_l2r_rx call seq_map_map(mapper = mapper_Fl2r, & av_s = irrig_l_av, & av_d = irrig_r_av, & fldlist = fields_to_remap, & norm = .true., & avwts_s = avwts_s, & avwtsfld_s = avwtsfld_s) allocate(irrig_normalized_r(lsize_r)) allocate(irrig_volr0_r(lsize_r)) call mct_aVect_exportRattr(irrig_r_av, irrig_normalized_field, irrig_normalized_r) call mct_aVect_exportRattr(irrig_r_av, irrig_volr0_field, irrig_volr0_r) ! ------------------------------------------------------------------------ ! Convert to a total irrigation flux on the ROF grid, and put this in the l2r_rx ! attribute vector ! ------------------------------------------------------------------------ allocate(irrig_flux_r(lsize_r)) do r = 1, lsize_r irrig_flux_r(r) = (irrig_normalized_r(r) * volr_r(r)) + irrig_volr0_r(r) end do call mct_aVect_importRattr(l2r_r, irrig_flux_field, irrig_flux_r) ! ------------------------------------------------------------------------ ! Clean up ! ------------------------------------------------------------------------ deallocate(volr_r) deallocate(volr_l) deallocate(irrig_flux_l) deallocate(irrig_flux_r) deallocate(irrig_normalized_l) deallocate(irrig_normalized_r) deallocate(irrig_volr0_l) deallocate(irrig_volr0_r) call mct_aVect_clean(irrig_l_av) call mct_aVect_clean(irrig_r_av) end subroutine map_lnd2rof_irrig subroutine map_rof2lnd_volr(volr_r, mapper_Fr2l, volr_l) !--------------------------------------------------------------- ! Description ! Map volr from the rof grid to the lnd grid. ! ! This is needed for the volr-normalization that is done in map_lnd2rof_irrig. ! ! Note that this mapping is also done in the course of mapping all rof -> lnd fields ! in prep_lnd_calc_r2x_lx. However, we do this mapping ourselves here for two reasons: ! ! (1) For the sake of this normalization, we change all volr < 0 values to 0; this is ! not done for the standard rof -> lnd mapping. ! ! (2) It's possible that the driver sequencing would be changed such that this rof -> ! lnd mapping happens before the lnd -> rof mapping. If that happened, then volr_l ! (i.e., volr that has been mapped to the land grid by prep_lnd_calc_r2x_lx) would ! be inconsistent with volr_r, which would be a Bad Thing for the ! volr-normalizated mapping (this mapping would no longer be conservative). So we ! do the rof -> lnd remapping here to ensure we have a volr_l that is consistent ! with volr_r. ! ! The pointer arguments to this routine should already be allocated to be the ! appropriate size. ! ! Arguments real(r8), pointer, intent(in) :: volr_r(:) ! river volume on the rof grid (input) type(seq_map) , intent(inout) :: mapper_Fr2l ! flux mapper for mapping rof -> lnd real(r8), pointer, intent(inout) :: volr_l(:) ! river volume on the lnd grid (output) (technically intent(in) since intent gives the association status of a pointer, but given as intent(inout) to avoid confusion, since its data are modified) ! ! Local variables integer :: lsize_r ! number of rof points integer :: lsize_l ! number of lnd points type(mct_avect) :: volr_r_av ! temporary attribute vector holding volr on the rof grid type(mct_avect) :: volr_l_av ! temporary attribute vector holding volr on the land grid ! This volr field name does not need to agree with the volr field name used in the ! 'real' attribute vectors character(len=*), parameter :: volr_field = 'volr' !--------------------------------------------------------------- SHR_ASSERT_FL(associated(volr_r), sourcefile, __LINE__) SHR_ASSERT_FL(associated(volr_l), sourcefile, __LINE__) lsize_r = size(volr_r) lsize_l = size(volr_l) call mct_aVect_init(volr_r_av, rList = volr_field, lsize = lsize_r) call mct_aVect_importRattr(volr_r_av, volr_field, volr_r) call mct_aVect_init(volr_l_av, rList = volr_field, lsize = lsize_l) ! This mapping uses the same options as the standard rof -> lnd mapping done in ! prep_lnd_calc_r2x_lx. If that mapping ever changed (e.g., introducing an avwts_s ! argument), then it's *possible* that we'd want this mapping to change, too. call seq_map_map(mapper = mapper_Fr2l, & av_s = volr_r_av, & av_d = volr_l_av, & fldlist = volr_field, & norm = .true.) call mct_aVect_exportRattr(volr_l_av, volr_field, volr_l) call mct_aVect_clean(volr_r_av) call mct_aVect_clean(volr_l_av) end subroutine map_rof2lnd_volr end module map_lnd2rof_irrig_mod