prep_glc_mod.F90
module prep_glc_mod
#include "shr_assert.h"
use shr_kind_mod , only: r8 => SHR_KIND_R8
use shr_kind_mod , only: cl => SHR_KIND_CL
use shr_sys_mod , only: shr_sys_abort, shr_sys_flush
use seq_comm_mct , only: num_inst_glc, num_inst_lnd, num_inst_frc, &
num_inst_ocn
use seq_comm_mct , only: CPLID, GLCID, logunit
use seq_comm_mct , only: seq_comm_getData=>seq_comm_setptrs
use seq_infodata_mod, only: seq_infodata_type, seq_infodata_getdata
use seq_map_type_mod
use seq_map_mod
use seq_flds_mod
use t_drv_timers_mod
use mct_mod
use perf_mod
use component_type_mod, only: component_get_x2c_cx, component_get_c2x_cx
use component_type_mod, only: component_get_dom_cx
use component_type_mod, only: glc, lnd, ocn
use glc_elevclass_mod, only : glc_get_num_elevation_classes, glc_elevclass_as_string
use glc_elevclass_mod, only : glc_all_elevclass_strings, GLC_ELEVCLASS_STRLEN
implicit none
save
private
!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------
public :: prep_glc_init
public :: prep_glc_mrg_lnd
public :: prep_glc_accum_lnd
public :: prep_glc_accum_ocn
public :: prep_glc_accum_avg
public :: prep_glc_calc_l2x_gx
public :: prep_glc_calc_o2x_gx
public :: prep_glc_zero_fields
public :: prep_glc_get_l2x_gx
public :: prep_glc_get_l2gacc_lx
public :: prep_glc_get_l2gacc_lx_one_instance
public :: prep_glc_get_l2gacc_lx_cnt
public :: prep_glc_get_o2x_gx
public :: prep_glc_get_x2gacc_gx
public :: prep_glc_get_x2gacc_gx_cnt
public :: prep_glc_get_mapper_Sl2g
public :: prep_glc_get_mapper_Fl2g
public :: prep_glc_get_mapper_So2g
public :: prep_glc_get_mapper_Fo2g
public :: prep_glc_calculate_subshelf_boundary_fluxes
!--------------------------------------------------------------------------
! Private interfaces
!--------------------------------------------------------------------------
private :: prep_glc_do_renormalize_smb
private :: prep_glc_set_g2x_lx_fields
private :: prep_glc_merge_lnd_forcing
private :: prep_glc_map_one_state_field_lnd2glc
private :: prep_glc_map_qice_conservative_lnd2glc
private :: prep_glc_renormalize_smb
!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------
! mappers
type(seq_map), pointer :: mapper_Sl2g
type(seq_map), pointer :: mapper_Fl2g
type(seq_map), pointer :: mapper_So2g
type(seq_map), pointer :: mapper_Fo2g
type(seq_map), pointer :: mapper_Fg2l
! attribute vectors
type(mct_aVect), pointer :: l2x_gx(:) ! Lnd export, glc grid, cpl pes - allocated in driver
type(mct_aVect), pointer :: o2x_gx(:) ! Ocn export, glc grid, cpl pes - allocated in driver
! accumulation variables
type(mct_aVect), pointer :: x2gacc_gx(:) ! Glc export, glc grid, cpl pes - allocated in driver
integer , target :: x2gacc_gx_cnt ! x2gacc_gx: number of time samples accumulated
type(mct_aVect), pointer :: l2gacc_lx(:) ! Lnd export, lnd grid, cpl pes - allocated in driver
integer , target :: l2gacc_lx_cnt ! l2gacc_lx: number of time samples accumulated
! other module variables
integer :: mpicom_CPLID ! MPI cpl communicator
! Whether to renormalize the SMB for conservation.
! Should be set to true for 2-way coupled runs with evolving ice sheets.
! Does not need to be true for 1-way coupling.
logical :: smb_renormalize
logical :: glc_present ! .true. => glc is present
! Name of flux field giving surface mass balance
character(len=*), parameter :: qice_fieldname = 'Flgl_qice'
! Names of some other fields
character(len=*), parameter :: Sg_frac_field = 'Sg_ice_covered'
character(len=*), parameter :: Sg_topo_field = 'Sg_topo'
character(len=*), parameter :: Sg_icemask_field = 'Sg_icemask'
! Fields needed in the g2x_lx attribute vector used as part of mapping qice from lnd to glc
character(len=:), allocatable :: g2x_lx_fields
type(mct_aVect), pointer :: o2gacc_ox(:) ! Ocn export, lnd grid, cpl pes - allocated in driver
integer , target :: o2gacc_ox_cnt ! number of time samples accumulated
real(r8), allocatable :: oceanTemperature(:)
real(r8), allocatable :: oceanSalinity(:)
real(r8), allocatable :: oceanHeatTransferVelocity(:)
real(r8), allocatable :: oceanSaltTransferVelocity(:)
real(r8), allocatable :: interfacePressure(:)
real(r8), allocatable :: iceTemperature(:)
real(r8), allocatable :: iceTemperatureDistance(:)
integer, allocatable :: iceFloatingMask(:)
real(r8), allocatable :: outInterfaceSalinity(:)
real(r8), allocatable :: outInterfaceTemperature(:)
real(r8), allocatable :: outFreshwaterFlux(:)
real(r8), allocatable :: outOceanHeatFlux(:)
real(r8), allocatable :: outIceHeatFlux(:)
!================================================================================================
contains
!================================================================================================
subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf)
!---------------------------------------------------------------
! Description
! Initialize module attribute vectors and mapping variables
!
! Arguments
type (seq_infodata_type) , intent(inout) :: infodata
logical , intent(in) :: lnd_c2_glc ! .true. => lnd to glc coupling on
logical , intent(in) :: ocn_c2_glcshelf ! .true. => ocn to glc coupling on
!
! Local Variables
integer :: eli, egi, eoi
integer :: lsize_l
integer :: lsize_g
integer :: lsize_o
logical :: samegrid_lg ! samegrid land and glc
logical :: samegrid_go ! .true. => samegrid ocean and glc
logical :: esmf_map_flag ! .true. => use esmf for mapping
logical :: iamroot_CPLID ! .true. => CPLID masterproc
logical :: do_hist_l2x1yrg ! .true. => create aux files: l2x 1yr glc forcings
character(CL) :: lnd_gnam ! lnd grid
character(CL) :: glc_gnam ! glc grid
character(CL) :: ocn_gnam ! ocn grid
type(mct_avect), pointer :: l2x_lx
type(mct_avect), pointer :: x2g_gx
type(mct_avect), pointer :: o2x_ox
character(*), parameter :: subname = '(prep_glc_init)'
character(*), parameter :: F00 = "('"//subname//" : ', 4A )"
!---------------------------------------------------------------
call seq_infodata_getData(infodata , &
esmf_map_flag=esmf_map_flag , &
glc_present=glc_present , &
histaux_l2x1yrg=do_hist_l2x1yrg, &
lnd_gnam=lnd_gnam , &
glc_gnam=glc_gnam , &
ocn_gnam=ocn_gnam)
allocate(mapper_Sl2g)
allocate(mapper_Fl2g)
allocate(mapper_So2g)
allocate(mapper_Fo2g)
allocate(mapper_Fg2l)
smb_renormalize = prep_glc_do_renormalize_smb(infodata)
if ((glc_present .and. lnd_c2_glc) .or. do_hist_l2x1yrg) then
l2x_lx => component_get_c2x_cx(lnd(1))
lsize_l = mct_aVect_lsize(l2x_lx)
allocate(l2gacc_lx(num_inst_lnd))
do eli = 1,num_inst_lnd
call mct_aVect_init(l2gacc_lx(eli), rList=seq_flds_l2x_fields_to_glc, lsize=lsize_l)
call mct_aVect_zero(l2gacc_lx(eli))
end do
l2gacc_lx_cnt = 0
end if
if (glc_present .and. lnd_c2_glc) then
call seq_comm_getData(CPLID, &
mpicom=mpicom_CPLID, iamroot=iamroot_CPLID)
x2g_gx => component_get_x2c_cx(glc(1))
lsize_g = mct_aVect_lsize(x2g_gx)
allocate(l2x_gx(num_inst_lnd))
do eli = 1,num_inst_lnd
call mct_aVect_init(l2x_gx(eli), rList=seq_flds_x2g_fields, lsize=lsize_g)
call mct_aVect_zero(l2x_gx(eli))
enddo
if (lnd_c2_glc) then
samegrid_lg = .true.
if (trim(lnd_gnam) /= trim(glc_gnam)) samegrid_lg = .false.
if (iamroot_CPLID) then
write(logunit,*) ' '
write(logunit,F00) 'Initializing mapper_Sl2g'
end if
call seq_map_init_rcfile(mapper_Sl2g, lnd(1), glc(1), &
'seq_maps.rc', 'lnd2glc_smapname:', 'lnd2glc_smaptype:', samegrid_lg, &
'mapper_Sl2g initialization', esmf_map_flag)
if (iamroot_CPLID) then
write(logunit,*) ' '
write(logunit,F00) 'Initializing mapper_Fl2g'
end if
call seq_map_init_rcfile(mapper_Fl2g, lnd(1), glc(1), &
'seq_maps.rc', 'lnd2glc_fmapname:', 'lnd2glc_fmaptype:', samegrid_lg, &
'mapper_Fl2g initialization', esmf_map_flag)
! We need to initialize our own Fg2l mapper because in some cases (particularly
! TG compsets - dlnd forcing CISM) the system doesn't otherwise create a Fg2l
! mapper.
if (iamroot_CPLID) then
write(logunit,*) ' '
write(logunit,F00) 'Initializing mapper_Fg2l'
end if
call seq_map_init_rcfile(mapper_Fg2l, glc(1), lnd(1), &
'seq_maps.rc', 'glc2lnd_fmapname:', 'glc2lnd_fmaptype:', samegrid_lg, &
'mapper_Fg2l initialization', esmf_map_flag)
call prep_glc_set_g2x_lx_fields()
end if
call shr_sys_flush(logunit)
end if
if (glc_present .and. ocn_c2_glcshelf) then
call seq_comm_getData(CPLID, &
mpicom=mpicom_CPLID, iamroot=iamroot_CPLID)
o2x_ox => component_get_c2x_cx(ocn(1))
lsize_o = mct_aVect_lsize(o2x_ox)
x2g_gx => component_get_x2c_cx(glc(1))
lsize_g = mct_aVect_lsize(x2g_gx)
allocate(o2x_gx(num_inst_ocn))
do eoi = 1,num_inst_ocn
call mct_aVect_init(o2x_gx(eoi), rList=seq_flds_o2x_fields, lsize=lsize_g)
call mct_aVect_zero(o2x_gx(eoi))
enddo
allocate(x2gacc_gx(num_inst_glc))
do egi = 1,num_inst_glc
call mct_aVect_init(x2gacc_gx(egi), x2g_gx, lsize_g)
call mct_aVect_zero(x2gacc_gx(egi))
end do
x2gacc_gx_cnt = 0
samegrid_go = .true.
if (trim(ocn_gnam) /= trim(glc_gnam)) samegrid_go = .false.
if (iamroot_CPLID) then
write(logunit,*) ' '
write(logunit,F00) 'Initializing mapper_So2g'
end if
call seq_map_init_rcfile(mapper_So2g, ocn(1), glc(1), &
'seq_maps.rc','ocn2glc_smapname:','ocn2glc_smaptype:',samegrid_go, &
'mapper_So2g initialization',esmf_map_flag)
if (iamroot_CPLID) then
write(logunit,*) ' '
write(logunit,F00) 'Initializing mapper_Fo2g'
end if
call seq_map_init_rcfile(mapper_Fo2g, ocn(1), glc(1), &
'seq_maps.rc','ocn2glc_fmapname:','ocn2glc_fmaptype:',samegrid_go, &
'mapper_Fo2g initialization',esmf_map_flag)
!Initialize module-level arrays associated with compute_melt_fluxes
allocate(oceanTemperature(lsize_g))
allocate(oceanSalinity(lsize_g))
allocate(oceanHeatTransferVelocity(lsize_g))
allocate(oceanSaltTransferVelocity(lsize_g))
allocate(interfacePressure(lsize_g))
allocate(iceTemperature(lsize_g))
allocate(iceTemperatureDistance(lsize_g))
allocate(iceFloatingMask(lsize_g))
allocate(outInterfaceSalinity(lsize_g))
allocate(outInterfaceTemperature(lsize_g))
allocate(outFreshwaterFlux(lsize_g))
allocate(outOceanHeatFlux(lsize_g))
allocate(outIceHeatFlux(lsize_g))
! TODO: Can we allocate these only while used or are we worried about performance hit?
! TODO: add deallocates!
call shr_sys_flush(logunit)
end if
end subroutine prep_glc_init
!================================================================================================
function prep_glc_do_renormalize_smb(infodata) result(do_renormalize_smb)
! Returns a logical saying whether we should do the smb renormalization
logical :: do_renormalize_smb ! function return value
!
! Arguments
type (seq_infodata_type) , intent(in) :: infodata
! Local variables
character(len=cl) :: glc_renormalize_smb ! namelist option saying whether to do smb renormalization
logical :: glc_coupled_fluxes ! does glc send fluxes to other components?
logical :: lnd_prognostic ! is lnd a prognostic component?
character(len=*), parameter :: subname = '(prep_glc_do_renormalize_smb)'
!---------------------------------------------------------------
call seq_infodata_getdata(infodata, &
glc_renormalize_smb = glc_renormalize_smb, &
glc_coupled_fluxes = glc_coupled_fluxes, &
lnd_prognostic = lnd_prognostic)
select case (glc_renormalize_smb)
case ('on')
do_renormalize_smb = .true.
case ('off')
do_renormalize_smb = .false.
case ('on_if_glc_coupled_fluxes')
if (.not. lnd_prognostic) then
! Do not renormalize if we're running glc with dlnd (T compsets): In this case
! there is no feedback from glc to lnd, and conservation is not important
do_renormalize_smb = .false.
else if (.not. glc_coupled_fluxes) then
! Do not renormalize if glc isn't sending fluxes to other components: In this
! case conservation is not important
do_renormalize_smb = .false.
else
! lnd_prognostic is true and glc_coupled_fluxes is true
do_renormalize_smb = .true.
end if
case default
write(logunit,*) subname,' ERROR: unknown value for glc_renormalize_smb: ', &
trim(glc_renormalize_smb)
call shr_sys_abort(subname//' ERROR: unknown value for glc_renormalize_smb')
end select
end function prep_glc_do_renormalize_smb
!================================================================================================
subroutine prep_glc_set_g2x_lx_fields()
!---------------------------------------------------------------
! Description
! Sets the module-level g2x_lx_fields variable.
!
! This gives the fields needed in the g2x_lx attribute vector used as part of mapping
! qice from lnd to glc.
!
! Local Variables
character(len=GLC_ELEVCLASS_STRLEN), allocatable :: all_elevclass_strings(:)
character(len=:), allocatable :: frac_fields
character(len=:), allocatable :: topo_fields
integer :: strlen
! 1 is probably enough, but use 10 to be safe, in case the length of the delimiter
! changes
integer, parameter :: extra_len_for_list_merge = 10
character(len=*), parameter :: subname = '(prep_glc_set_g2x_lx_fields)'
!---------------------------------------------------------------
allocate(all_elevclass_strings(0:glc_get_num_elevation_classes()))
all_elevclass_strings = glc_all_elevclass_strings(include_zero = .true.)
frac_fields = shr_string_listFromSuffixes( &
suffixes = all_elevclass_strings, &
strBase = Sg_frac_field)
! Sg_topo is not actually needed on the land grid in
! prep_glc_map_qice_conservative_lnd2glc, but it is required by the current interface
! for map_glc2lnd_ec.
topo_fields = shr_string_listFromSuffixes( &
suffixes = all_elevclass_strings, &
strBase = Sg_topo_field)
strlen = len_trim(frac_fields) + len_trim(topo_fields) + extra_len_for_list_merge
allocate(character(len=strlen) :: g2x_lx_fields)
call shr_string_listMerge(frac_fields, topo_fields, g2x_lx_fields)
end subroutine prep_glc_set_g2x_lx_fields
!================================================================================================
subroutine prep_glc_accum_lnd(timer)
!---------------------------------------------------------------
! Description
! Accumulate glc inputs from lnd
!
! Arguments
character(len=*), intent(in) :: timer
!
! Local Variables
integer :: eli
type(mct_avect), pointer :: l2x_lx
character(*), parameter :: subname = '(prep_glc_accum_lnd)'
!---------------------------------------------------------------
call t_drvstartf (trim(timer),barrier=mpicom_CPLID)
do eli = 1,num_inst_lnd
l2x_lx => component_get_c2x_cx(lnd(eli))
if (l2gacc_lx_cnt == 0) then
call mct_avect_copy(l2x_lx, l2gacc_lx(eli))
else
call mct_avect_accum(l2x_lx, l2gacc_lx(eli))
endif
end do
l2gacc_lx_cnt = l2gacc_lx_cnt + 1
call t_drvstopf (trim(timer))
end subroutine prep_glc_accum_lnd
!================================================================================================
subroutine prep_glc_accum_ocn(timer)
!---------------------------------------------------------------
! Description
! Accumulate glc inputs from ocn
!
! Arguments
character(len=*), intent(in) :: timer
!
! Local Variables
integer :: egi
type(mct_avect), pointer :: x2g_gx
character(*), parameter :: subname = '(prep_glc_accum_ocn)'
!---------------------------------------------------------------
call t_drvstartf (trim(timer),barrier=mpicom_CPLID)
do egi = 1,num_inst_glc
x2g_gx => component_get_x2c_cx(glc(egi))
if (x2gacc_gx_cnt == 0) then
call mct_avect_copy(x2g_gx, x2gacc_gx(egi))
else
call mct_avect_accum(x2g_gx, x2gacc_gx(egi))
endif
end do
x2gacc_gx_cnt = x2gacc_gx_cnt + 1
call t_drvstopf (trim(timer))
end subroutine prep_glc_accum_ocn
!================================================================================================
subroutine prep_glc_accum_avg(timer, lnd2glc_averaged_now)
!---------------------------------------------------------------
! Description
! Finalize accumulation of glc inputs
! Note: There could be separate accum_avg routines for forcing coming
! from each component (LND and OCN), but they can be combined here
! by taking advantage of l2gacc_lx_cnt and x2gacc_gx_cnt variables
! that will only be greater than 0 if corresponding coupling is enabled.
!
! Arguments
character(len=*), intent(in) :: timer
logical, intent(inout) :: lnd2glc_averaged_now ! Set to .true. if lnd2glc averages were taken this timestep (otherwise left unchanged)
!
! Local Variables
integer :: eli, egi
type(mct_avect), pointer :: x2g_gx
character(*), parameter :: subname = '(prep_glc_accum_avg)'
!---------------------------------------------------------------
! Accumulation for LND
call t_drvstartf (trim(timer),barrier=mpicom_CPLID)
if (l2gacc_lx_cnt > 0) then
lnd2glc_averaged_now = .true.
end if
if (l2gacc_lx_cnt > 1) then
do eli = 1,num_inst_lnd
call mct_avect_avg(l2gacc_lx(eli), l2gacc_lx_cnt)
end do
end if
l2gacc_lx_cnt = 0
! Accumulation for OCN
if (x2gacc_gx_cnt > 1) then
do egi = 1,num_inst_glc
! temporary formation of average
call mct_avect_avg(x2gacc_gx(egi), x2gacc_gx_cnt)
! ***NOTE***THE FOLLOWING ACTUALLY MODIFIES x2g_gx
x2g_gx => component_get_x2c_cx(glc(egi))
call mct_avect_copy(x2gacc_gx(egi), x2g_gx)
enddo
end if
x2gacc_gx_cnt = 0
call t_drvstopf (trim(timer))
end subroutine prep_glc_accum_avg
!================================================================================================
subroutine prep_glc_mrg_lnd(infodata, fractions_gx, timer_mrg)
!---------------------------------------------------------------
! Description
! Merge glc inputs
!
! Arguments
type(seq_infodata_type) , intent(in) :: infodata
type(mct_aVect) , intent(in) :: fractions_gx(:)
character(len=*) , intent(in) :: timer_mrg
!
! Local Variables
integer :: egi, eli, efi
type(mct_avect), pointer :: x2g_gx
character(*), parameter :: subname = '(prep_glc_mrg_lnd)'
!---------------------------------------------------------------
call t_drvstartf (trim(timer_mrg),barrier=mpicom_CPLID)
do egi = 1,num_inst_glc
! Use fortran mod to address ensembles in merge
eli = mod((egi-1),num_inst_lnd) + 1
efi = mod((egi-1),num_inst_frc) + 1
x2g_gx => component_get_x2c_cx(glc(egi))
call prep_glc_merge_lnd_forcing(l2x_gx(eli), fractions_gx(efi), x2g_gx)
enddo
call t_drvstopf (trim(timer_mrg))
end subroutine prep_glc_mrg_lnd
!================================================================================================
subroutine prep_glc_merge_lnd_forcing( l2x_g, fractions_g, x2g_g )
!-----------------------------------------------------------------------
! Description
! "Merge" land forcing for glc input.
!
! State fields are copied directly, meaning that averages are taken just over the
! land-covered portion of the glc domain.
!
! Flux fields are downweighted by landfrac, which effectively sends a 0 flux from the
! non-land-covered portion of the glc domain.
!
! Arguments
type(mct_aVect), intent(inout) :: l2x_g ! input
type(mct_aVect), intent(in) :: fractions_g
type(mct_aVect), intent(inout) :: x2g_g ! output
!-----------------------------------------------------------------------
integer :: num_flux_fields
integer :: num_state_fields
integer :: nflds
integer :: i,n
integer :: mrgstr_index
integer :: index_l2x
integer :: index_x2g
integer :: index_lfrac
integer :: lsize
logical :: iamroot
logical, save :: first_time = .true.
character(CL),allocatable :: mrgstr(:) ! temporary string
character(CL) :: field ! string converted to char
character(*), parameter :: subname = '(prep_glc_merge_lnd_forcing) '
!-----------------------------------------------------------------------
call seq_comm_getdata(CPLID, iamroot=iamroot)
lsize = mct_aVect_lsize(x2g_g)
num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_fluxes_from_lnd))
num_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_states_from_lnd))
if (first_time) then
nflds = num_flux_fields + num_state_fields
allocate(mrgstr(nflds))
end if
mrgstr_index = 1
do i = 1, num_state_fields
call seq_flds_getField(field, i, seq_flds_x2g_states)
index_l2x = mct_aVect_indexRA(l2x_g, trim(field))
index_x2g = mct_aVect_indexRA(x2g_g, trim(field))
if (first_time) then
mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // &
' = l2x%'//trim(field)
end if
do n = 1, lsize
x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n)
end do
mrgstr_index = mrgstr_index + 1
enddo
index_lfrac = mct_aVect_indexRA(fractions_g,"lfrac")
do i = 1, num_flux_fields
call seq_flds_getField(field, i, seq_flds_x2g_fluxes_from_lnd)
index_l2x = mct_aVect_indexRA(l2x_g, trim(field))
index_x2g = mct_aVect_indexRA(x2g_g, trim(field))
if (trim(field) == qice_fieldname) then
if (first_time) then
mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // &
' = l2x%'//trim(field)
end if
! treat qice as if it were a state variable, with a simple copy.
do n = 1, lsize
x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n)
end do
else
write(logunit,*) subname,' ERROR: Flux fields other than ', &
qice_fieldname, ' currently are not handled in lnd2glc remapping.'
write(logunit,*) '(Attempt to handle flux field <', trim(field), '>.)'
write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes'
write(logunit,*) 'in a smooth, conservative manner.'
call shr_sys_abort(subname//&
' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.')
endif ! qice_fieldname
mrgstr_index = mrgstr_index + 1
end do
if (first_time) then
if (iamroot) then
write(logunit,'(A)') subname//' Summary:'
do i = 1,nflds
write(logunit,'(A)') trim(mrgstr(i))
enddo
endif
deallocate(mrgstr)
endif
first_time = .false.
end subroutine prep_glc_merge_lnd_forcing
subroutine prep_glc_calc_o2x_gx(timer)
!---------------------------------------------------------------
! Description
! Create o2x_gx
! Arguments
character(len=*), intent(in) :: timer
character(*), parameter :: subname = '(prep_glc_calc_o2x_gx)'
! Local Variables
integer eoi
type(mct_avect), pointer :: o2x_ox
call t_drvstartf (trim(timer),barrier=mpicom_CPLID)
do eoi = 1,num_inst_ocn
o2x_ox => component_get_c2x_cx(ocn(eoi))
call seq_map_map(mapper_So2g, o2x_ox, o2x_gx(eoi), &
fldlist=seq_flds_x2g_states_from_ocn,norm=.true.)
enddo
call t_drvstopf (trim(timer))
end subroutine prep_glc_calc_o2x_gx
!================================================================================================
!================================================================================================
subroutine prep_glc_calc_l2x_gx(fractions_lx, timer)
!---------------------------------------------------------------
! Description
! Create l2x_gx (note that l2x_gx is a local module variable)
! Also l2x_gx is really the accumulated l2xacc_lx mapped to l2x_gx
!
use shr_string_mod, only : shr_string_listGetNum
! Arguments
type(mct_aVect) , intent(in) :: fractions_lx(:)
character(len=*), intent(in) :: timer
!
! Local Variables
integer :: egi, eli, efi
integer :: num_flux_fields
integer :: num_state_fields
integer :: field_num
character(len=cl) :: fieldname
character(*), parameter :: subname = '(prep_glc_calc_l2x_gx)'
!---------------------------------------------------------------
call t_drvstartf (trim(timer),barrier=mpicom_CPLID)
num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_fluxes_from_lnd))
num_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_states_from_lnd))
do egi = 1,num_inst_glc
! Use fortran mod to address ensembles in merge
eli = mod((egi-1),num_inst_lnd) + 1
efi = mod((egi-1),num_inst_frc) + 1
do field_num = 1, num_flux_fields
call seq_flds_getField(fieldname, field_num, seq_flds_x2g_fluxes_from_lnd)
if (trim(fieldname) == qice_fieldname) then
! Use a bilinear (Sl2g) mapper, as for states.
! The Fg2l mapper is needed to map some glc fields to the land grid
! for purposes of conservation.
call prep_glc_map_qice_conservative_lnd2glc(egi=egi, eli=eli, &
fractions_lx = fractions_lx(efi), &
mapper_Sl2g = mapper_Sl2g, &
mapper_Fg2l = mapper_Fg2l)
else
write(logunit,*) subname,' ERROR: Flux fields other than ', &
qice_fieldname, ' currently are not handled in lnd2glc remapping.'
write(logunit,*) '(Attempt to handle flux field <', trim(fieldname), '>.)'
write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes'
write(logunit,*) 'in a smooth, conservative manner.'
call shr_sys_abort(subname//&
' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.')
endif ! qice_fieldname
end do
do field_num = 1, num_state_fields
call seq_flds_getField(fieldname, field_num, seq_flds_x2g_states_from_lnd)
call prep_glc_map_one_state_field_lnd2glc(egi=egi, eli=eli, &
fieldname = fieldname, &
fractions_lx = fractions_lx(efi), &
mapper = mapper_Sl2g)
end do
enddo ! egi
call t_drvstopf (trim(timer))
end subroutine prep_glc_calc_l2x_gx
!================================================================================================
subroutine prep_glc_map_one_state_field_lnd2glc(egi, eli, fieldname, fractions_lx, mapper)
! Maps a single field from the land grid to the glc grid.
!
! This mapping is not conservative, so should only be used for state fields.
!
! NOTE(wjs, 2017-05-10) We used to map each field separately because each field needed
! its own vertical gradient calculator. Now that we don't need vertical gradient
! calculators, we may be able to change this to map multiple fields at once, at least
! for part of map_lnd2glc.
use map_lnd2glc_mod, only : map_lnd2glc
! Arguments
integer, intent(in) :: egi ! glc instance index
integer, intent(in) :: eli ! lnd instance index
character(len=*), intent(in) :: fieldname ! base name of field to map (without elevation class suffix)
type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid, for this frac instance
type(seq_map), intent(inout) :: mapper
!
! Local Variables
type(mct_avect), pointer :: g2x_gx ! glc export, glc grid, cpl pes - allocated in driver
!---------------------------------------------------------------
g2x_gx => component_get_c2x_cx(glc(egi))
call map_lnd2glc(l2x_l = l2gacc_lx(eli), &
landfrac_l = fractions_lx, &
g2x_g = g2x_gx, &
fieldname = fieldname, &
mapper = mapper, &
l2x_g = l2x_gx(eli))
end subroutine prep_glc_map_one_state_field_lnd2glc
!================================================================================================
subroutine prep_glc_calculate_subshelf_boundary_fluxes
!---------------------------------------------------------------
! Description
! On the ice sheet grid, calculate shelf boundary fluxes
use shr_const_mod , only: SHR_CONST_KAPPA_LAND_ICE
! Local Variables
integer :: gsize, n, egi
type(mct_aVect), pointer :: o2x_ox ! Ocn export, ocn grid, cpl pes
type(mct_aVect), pointer :: x2g_gx ! Glc import, glc grid, cpl pes
type(mct_aVect), pointer :: g2x_gx ! Glc import, glc grid, cpl pes
integer :: index_x2g_So_blt
integer :: index_x2g_So_bls
integer :: index_x2g_So_htv
integer :: index_x2g_So_stv
integer :: index_x2g_So_rhoeff
integer :: index_g2x_Sg_tbot
integer :: index_g2x_Sg_dztbot
integer :: index_g2x_Sg_lithop
integer :: index_g2x_Sg_icemask_floating
integer :: index_g2x_Sg_blis
integer :: index_g2x_Sg_blit
integer :: index_g2x_Fogx_qiceho
integer :: index_g2x_Fogx_qicelo
integer :: index_x2g_Fogx_qiceli
integer :: index_x2g_Fogx_qicehi
character(*), parameter :: subname = '(prep_glc_calculate_subshelf_boundary_fluxes)'
!---------------------------------------------------------------
if (.not.(glc_present)) return
do egi = 1,num_inst_glc
o2x_ox => component_get_c2x_cx(ocn(egi))
g2x_gx => component_get_c2x_cx(glc(egi))
x2g_gx => component_get_x2c_cx(glc(egi))
!Remap relevant ocean variables to ice sheet grid.
!Done here instead of in glc-frequency mapping so it happens within ocean coupling interval.
! Also could map o2x_ox->o2x_gx(1) but using x2g_gx as destination allows us to see
! these fields on the GLC grid of the coupler history file, which helps with debugging.
call seq_map_map(mapper_So2g, o2x_ox, x2g_gx, &
fldlist=seq_flds_x2g_states_from_ocn,norm=.true.)
! inputs to melt flux calculation
index_x2g_So_blt = mct_avect_indexra(x2g_gx,'So_blt',perrwith='quiet')
index_x2g_So_bls = mct_avect_indexra(x2g_gx,'So_bls',perrwith='quiet')
index_x2g_So_htv = mct_avect_indexra(x2g_gx,'So_htv',perrwith='quiet')
index_x2g_So_stv = mct_avect_indexra(x2g_gx,'So_stv',perrwith='quiet')
index_x2g_So_rhoeff = mct_avect_indexra(x2g_gx,'So_rhoeff',perrwith='quiet')
index_g2x_Sg_tbot = mct_avect_indexra(g2x_gx,'Sg_tbot',perrwith='quiet')
index_g2x_Sg_dztbot = mct_avect_indexra(g2x_gx,'Sg_dztbot',perrwith='quiet')
index_g2x_Sg_lithop = mct_avect_indexra(g2x_gx,'Sg_lithop',perrwith='quiet')
index_g2x_Sg_icemask_floating = mct_avect_indexra(g2x_gx,'Sg_icemask_floating',perrwith='quiet')
! outputs to melt flux calculation
index_g2x_Sg_blis = mct_avect_indexra(g2x_gx,'Sg_blis',perrwith='quiet')
index_g2x_Sg_blit = mct_avect_indexra(g2x_gx,'Sg_blit',perrwith='quiet')
index_g2x_Fogx_qiceho = mct_avect_indexra(g2x_gx,'Fogx_qiceho',perrwith='quiet')
index_g2x_Fogx_qicelo = mct_avect_indexra(g2x_gx,'Fogx_qicelo',perrwith='quiet')
index_x2g_Fogx_qiceli = mct_avect_indexra(x2g_gx,'Fogx_qiceli',perrwith='quiet')
index_x2g_Fogx_qicehi = mct_avect_indexra(x2g_gx,'Fogx_qicehi',perrwith='quiet')
gsize = mct_aVect_lsize(g2x_gx)
do n=1,gsize
!Extract glc and ocn-sourced coupler fields used as input to compute_melt_fluxes to local arrays...
! Fields from the ocean, now on the GLC grid
oceanTemperature(n) = x2g_gx%rAttr(index_x2g_So_blt,n)
oceanSalinity(n) = x2g_gx%rAttr(index_x2g_So_bls,n)
oceanHeatTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_htv,n)
oceanSaltTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_stv,n)
! Fields from the ice sheet model (still on the GLC grid)
iceTemperature(n) = g2x_gx%rAttr(index_g2x_Sg_tbot,n)
iceTemperatureDistance(n) = g2x_gx%rAttr(index_g2x_Sg_dztbot,n)
interfacePressure(n) = g2x_gx%rAttr(index_g2x_Sg_lithop,n)
iceFloatingMask(n) = g2x_gx%rAttr(index_g2x_Sg_icemask_floating,n)
!... initialize local compute_melt_fluxes output arrays...
outInterfaceSalinity(n) = 0.0_r8
outInterfaceTemperature(n) = 0.0_r8
outFreshwaterFlux(n) = 0.0_r8
outOceanHeatFlux(n) = 0.0_r8
outIceHeatFlux(n) = 0.0_r8
end do
!...calculate fluxes...
call compute_melt_fluxes(oceanTemperature=oceanTemperature,&
oceanSalinity=oceanSalinity,&
oceanHeatTransferVelocity=oceanHeatTransferVelocity,&
oceanSaltTransferVelocity=oceanSaltTransferVelocity,&
interfacePressure=interfacePressure,&
iceTemperature=iceTemperature,&
iceTemperatureDistance=iceTemperatureDistance, &
iceFloatingMask=iceFloatingMask, &
outInterfaceSalinity=outInterfaceSalinity,&
outInterfaceTemperature=outInterfaceTemperature,&
outFreshwaterFlux=outFreshwaterFlux,&
outOceanHeatFlux=outOceanHeatFlux,&
outIceHeatFlux=outIceHeatFlux,&
gsize=gsize)
!...and assign fluxes to glc and ocn-directed coupler fields
do n=1,gsize
!Assign outputs from compute_melt_fluxes back into coupler attributes
g2x_gx%rAttr(index_g2x_Sg_blis,n) = outInterfaceSalinity(n) !to ocean
g2x_gx%rAttr(index_g2x_Sg_blit,n) = outInterfaceTemperature(n) !to ocean
g2x_gx%rAttr(index_g2x_Fogx_qiceho,n) = outOceanHeatFlux(n) !to ocean
g2x_gx%rAttr(index_g2x_Fogx_qicelo,n)= outFreshwaterFlux(n) !to ocean
x2g_gx%rAttr(index_x2g_Fogx_qicehi,n) = outIceHeatFlux(n) !to ice sheet
x2g_gx%rAttr(index_x2g_Fogx_qiceli,n) = -1.0_r8 * outFreshwaterFlux(n) !to ice sheet
end do
!Note: remap ocean-side outputs back onto ocean grid done in call to prep_ocn_shelf_calc_g2x_ox
end do ! loop over GLC instances
end subroutine prep_glc_calculate_subshelf_boundary_fluxes
!================================================================================================
subroutine prep_glc_zero_fields()
!---------------------------------------------------------------
! Description
! Set glc inputs to zero
!
! This is appropriate during time intervals when we're not sending valid data to glc.
! In principle we shouldn't need to zero the fields at these times (instead, glc
! should just ignore the fields at these times). However, some tests (like an ERS or
! ERI test that stops the final run segment mid-year) can fail if we don't explicitly
! zero the fields, because these x2g fields can then differ upon restart.
! Local Variables
integer :: egi
type(mct_avect), pointer :: x2g_gx
!---------------------------------------------------------------
do egi = 1,num_inst_glc
x2g_gx => component_get_x2c_cx(glc(egi))
call mct_aVect_zero(x2g_gx)
end do
end subroutine prep_glc_zero_fields
!================================================================================================
subroutine prep_glc_map_qice_conservative_lnd2glc(egi, eli, fractions_lx, &
mapper_Sl2g, mapper_Fg2l)
! Maps the surface mass balance field (qice) from the land grid to the glc grid.
!
! Use a smooth, non-conservative (bilinear) mapping, followed by a correction for
! conservation.
!
! For high-level design, see:
! https://docs.google.com/document/d/1H_SuK6SfCv1x6dK91q80dFInPbLYcOkUj_iAa6WRnqQ/edit
use map_lnd2glc_mod, only : map_lnd2glc
! Arguments
integer, intent(in) :: egi ! glc instance index
integer, intent(in) :: eli ! lnd instance index
type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid, for this frac instance
type(seq_map), intent(inout) :: mapper_Sl2g ! state mapper from land to glc grid; non-conservative
type(seq_map), intent(inout) :: mapper_Fg2l ! flux mapper from glc to land grid; conservative
!
! Local Variables
type(mct_aVect), pointer :: g2x_gx ! glc export, glc grid
logical :: iamroot
!Note: The sums in this subroutine use the coupler areas aream_l and aream_g.
! The coupler areas can differ from the native areas area_l and area_g.
! (For CISM with a polar stereographic projection, area_g can differ from aream_g
! by up to ~10%.)
! If so, then the calls to subroutine mct_avect_vecmult in component_mod.F90
! (just before and after the call to comp_run) should adjust the SMB fluxes
! such that in each grid cell, the native value of area*flux is equal to the
! coupler value of aream*flux. This assumes that the SMB field is contained in
! seq_fields l2x_fluxes and seq_fields_x2g_fluxes.
real(r8), dimension(:), allocatable :: aream_g ! cell areas on glc grid, for mapping
real(r8), dimension(:), allocatable :: area_g ! cell areas on glc grid, according to glc model
type(mct_ggrid), pointer :: dom_g ! glc grid info
integer :: lsize_g ! number of points on glc grid
integer :: n
integer :: km, ka
real(r8), pointer :: qice_g(:) ! qice data on glc grid
!---------------------------------------------------------------
call seq_comm_getdata(CPLID, iamroot=iamroot)
if (iamroot) then
write(logunit,*) ' '
write(logunit,*) 'In prep_glc_map_qice_conservative_lnd2glc'
write(logunit,*) 'smb_renormalize = ', smb_renormalize
endif
! Get attribute vector needed for mapping and conservation
g2x_gx => component_get_c2x_cx(glc(egi))
! get grid size
lsize_g = mct_aVect_lsize(l2x_gx(eli))
! allocate and fill area arrays on the glc grid
! (Note that we get domain information from instance 1, following what's done in
! other parts of the coupler.)
dom_g => component_get_dom_cx(glc(1))
allocate(aream_g(lsize_g))
km = mct_aVect_indexRa(dom_g%data, "aream" )
aream_g(:) = dom_g%data%rAttr(km,:)
allocate(area_g(lsize_g))
ka = mct_aVect_indexRa(dom_g%data, "area" )
area_g(:) = dom_g%data%rAttr(ka,:)
! Map the SMB from the land grid to the glc grid, using a non-conservative state mapper.
call map_lnd2glc(l2x_l = l2gacc_lx(eli), &
landfrac_l = fractions_lx, &
g2x_g = g2x_gx, &
fieldname = qice_fieldname, &
mapper = mapper_Sl2g, &
l2x_g = l2x_gx(eli))
! Export the remapped SMB to a local array
allocate(qice_g(lsize_g))
call mct_aVect_exportRattr(l2x_gx(eli), trim(qice_fieldname), qice_g)
! Make a preemptive adjustment to qice_g to account for area differences between CISM and the coupler.
! In component_mod.F90, there is a call to mct_avect_vecmult, which multiplies the fluxes
! by aream_g/area_g for conservation purposes. Where CISM areas are larger (area_g > aream_g),
! the fluxes are reduced, and where CISM areas are smaller, the fluxes are increased.
! As a result, an SMB of 1 m/yr in CLM would be converted to an SMB ranging from
! ~0.9 to 1.05 m/yr in CISM (with smaller values where CISM areas are larger, and larger
! values where CISM areas are smaller).
! Here, to keep CISM values close to the CLM values in the corresponding locations,
! we anticipate the later correction and multiply qice_g by area_g/aream_g.
! Then the later call to mct_avect_vecmult will bring qice back to the original values
! obtained from bilinear remapping.
! If Flgl_qice were changed to a state (and not included in seq_flds_x2g_fluxes),
! then we could skip this adjustment.
!
! Note that we are free to do this or any other adjustments we want to qice at this
! point in the remapping, because the conservation correction will ensure that we
! still conserve globally despite these adjustments (and smb_renormalize = .false.
! should only be used in cases where conservation doesn't matter anyway).
do n = 1, lsize_g
if (aream_g(n) > 0.0_r8) then
qice_g(n) = qice_g(n) * area_g(n)/aream_g(n)
else
qice_g(n) = 0.0_r8
endif
enddo
if (smb_renormalize) then
call prep_glc_renormalize_smb( &
eli = eli, &
fractions_lx = fractions_lx, &
g2x_gx = g2x_gx, &
mapper_Fg2l = mapper_Fg2l, &
aream_g = aream_g, &
qice_g = qice_g)
end if
! Put the adjusted SMB back into l2x_gx.
!
! If we are doing renormalization, then this is the renormalized SMB. Whether or not
! we are doing renormalization, this captures the preemptive adjustment to qice_g to
! account for area differences between CISM and the coupler.
call mct_aVect_importRattr(l2x_gx(eli), qice_fieldname, qice_g)
! clean up
deallocate(aream_g)
deallocate(area_g)
deallocate(qice_g)
end subroutine prep_glc_map_qice_conservative_lnd2glc
!================================================================================================
subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, aream_g, qice_g)
! Renormalizes surface mass balance (smb, here named qice_g) so that the global
! integral on the glc grid is equal to the global integral on the land grid.
!
! This is required for conservation - although conservation is only necessary if we
! are running with a fully-interactive, two-way-coupled glc.
!
! For high-level design, see:
! https://docs.google.com/document/d/1H_SuK6SfCv1x6dK91q80dFInPbLYcOkUj_iAa6WRnqQ/edit
use map_glc2lnd_mod, only : map_glc2lnd_ec
! Arguments
integer , intent(in) :: eli ! lnd instance index
type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid, for this frac instance
type(mct_aVect) , intent(in) :: g2x_gx ! glc export, glc grid
type(seq_map) , intent(inout) :: mapper_Fg2l ! flux mapper from glc to land grid; conservative
real(r8) , intent(in) :: aream_g(:) ! cell areas on glc grid, for mapping
real(r8) , intent(inout) :: qice_g(:) ! qice data on glc grid
!
! Local Variables
integer :: mpicom
logical :: iamroot
type(mct_ggrid), pointer :: dom_l ! land grid info
integer :: lsize_l ! number of points on land grid
integer :: lsize_g ! number of points on glc grid
real(r8), dimension(:), allocatable :: aream_l ! cell areas on land grid, for mapping
real(r8), pointer :: qice_l(:,:) ! SMB (Flgl_qice) on land grid
real(r8), pointer :: frac_l(:,:) ! EC fractions (Sg_ice_covered) on land grid
real(r8), pointer :: tmp_field_l(:) ! temporary field on land grid
! The following need to be pointers to satisfy the MCT interface
! Note: Sg_icemask defines where the ice sheet model can receive a nonzero SMB from the land model.
real(r8), pointer :: Sg_icemask_g(:) ! icemask on glc grid
real(r8), pointer :: Sg_icemask_l(:) ! icemask on land grid
real(r8), pointer :: lfrac(:) ! land fraction on land grid
type(mct_aVect) :: g2x_lx ! glc export, lnd grid (not a pointer: created locally)
type(mct_avect) :: Sg_icemask_l_av ! temporary attribute vector holding Sg_icemask on the land grid
integer :: nEC ! number of elevation classes
integer :: n
integer :: ec
integer :: km
! various strings for building field names
character(len=:), allocatable :: elevclass_as_string
character(len=:), allocatable :: qice_field
character(len=:), allocatable :: frac_field
! local and global sums of accumulation and ablation; used to compute renormalization factors
real(r8) :: local_accum_on_land_grid
real(r8) :: global_accum_on_land_grid
real(r8) :: local_accum_on_glc_grid
real(r8) :: global_accum_on_glc_grid
real(r8) :: local_ablat_on_land_grid
real(r8) :: global_ablat_on_land_grid
real(r8) :: local_ablat_on_glc_grid
real(r8) :: global_ablat_on_glc_grid
! renormalization factors (should be close to 1, e.g. in range 0.95 to 1.05)
real(r8) :: accum_renorm_factor ! ratio between global accumulation on the two grids
real(r8) :: ablat_renorm_factor ! ratio between global ablation on the two grids
real(r8) :: effective_area ! grid cell area multiplied by min(lfrac,Sg_icemask_l).
! This is the area that can contribute SMB to the ice sheet model.
!---------------------------------------------------------------
lsize_g = size(qice_g)
SHR_ASSERT_FL((size(aream_g) == lsize_g), __FILE__, __LINE__)
call seq_comm_setptrs(CPLID, mpicom=mpicom)
call seq_comm_getdata(CPLID, iamroot=iamroot)
lsize_l = mct_aVect_lsize(l2gacc_lx(eli))
! allocate and fill area arrays on the land grid
! (Note that we get domain information from instance 1, following what's done in
! other parts of the coupler.)
dom_l => component_get_dom_cx(lnd(1))
allocate(aream_l(lsize_l))
km = mct_aVect_indexRa(dom_l%data, "aream" )
aream_l(:) = dom_l%data%rAttr(km,:)
! Export land fractions from fractions_lx to a local array
allocate(lfrac(lsize_l))
call mct_aVect_exportRattr(fractions_lx, "lfrac", lfrac)
! Map Sg_icemask from the glc grid to the land grid.
! This may not be necessary, if Sg_icemask_l has already been mapped from Sg_icemask_g.
! It is done here for two reasons:
! (1) The mapping will *not* have been done if we are running with dlnd (e.g., a TG case).
! (2) Because of coupler lags, the current Sg_icemask_l might not be up to date with
! Sg_icemask_g. This probably isn't a problem in practice, but doing the mapping
! here ensures the mask is up to date.
!
! This mapping uses the same options as the standard glc -> lnd mapping done in
! prep_lnd_calc_g2x_lx. If that mapping ever changed (e.g., changing norm to
! .false.), then we should change this mapping, too.
!
! BUG(wjs, 2017-05-11, #1516) I think we actually want norm = .false. here, but this
! requires some more thought
call mct_aVect_init(Sg_icemask_l_av, rList = Sg_icemask_field, lsize = lsize_l)
call seq_map_map(mapper = mapper_Fg2l, &
av_s = g2x_gx, &
av_d = Sg_icemask_l_av, &
fldlist = Sg_icemask_field, &
norm = .true.)
! Export Sg_icemask_l from the temporary attribute vector to a local array
allocate(Sg_icemask_l(lsize_l))
call mct_aVect_exportRattr(Sg_icemask_l_av, Sg_icemask_field, Sg_icemask_l)
! Clean the temporary attribute vector
call mct_aVect_clean(Sg_icemask_l_av)
! Map Sg_ice_covered from the glc grid to the land grid.
! This gives the fields Sg_ice_covered00, Sg_ice_covered01, etc. on the land grid.
! These fields are needed to integrate the total SMB on the land grid, for conservation purposes.
! As above, the mapping may not be necessary, because Sg_ice_covered might already have been mapped.
! However, the mapping will not have been done in a TG case with dlnd, and it might not
! be up to date because of coupler lags (though the latter probably isn't a problem
! in practice).
!
! Note that, for a case with full two-way coupling, we will only conserve if the
! actual land cover used over the course of the year matches these currently-remapped
! values. This should generally be the case with the current coupling setup.
!
! One could argue that it would be safer (for conservation purposes) if LND sent its
! grid cell average SMB values, or if it sent its own notion of the area in each
! elevation class for the purpose of creating grid cell average SMB values here. But
! these options cause problems if we're not doing full two-way coupling (e.g., in a TG
! case with dlnd, or in the common case where GLC is a diagnostic component that
! doesn't cause updates in the glacier areas in LND). In these cases without full
! two-way coupling, if we use the LND's notion of the area in each elevation class,
! then the conservation corrections would end up correcting for discrepancies in
! elevation class areas between LND and GLC, rather than just correcting for
! discrepancies arising from the remapping of SMB. (And before you get worried: It
! doesn't matter that we are not conserving in these cases without full two-way
! coupling, because GLC isn't connected with the rest of the system in terms of energy
! and mass in these cases. So in these cases, it's okay that the LND integral computed
! here differs from the integral that LND itself would compute.)
! Create an attribute vector g2x_lx to hold the mapped fields
call mct_aVect_init(g2x_lx, rList=g2x_lx_fields, lsize=lsize_l)
! Map Sg_ice_covered and Sg_topo from glc to land
call map_glc2lnd_ec( &
g2x_g = g2x_gx, &
frac_field = Sg_frac_field, &
topo_field = Sg_topo_field, &
icemask_field = Sg_icemask_field, &
extra_fields = ' ', & ! no extra fields
mapper = mapper_Fg2l, &
g2x_l = g2x_lx)
! Export qice and Sg_ice_covered in each elevation class to local arrays.
! Note: qice comes from l2gacc_lx; frac comes from g2x_lx.
nEC = glc_get_num_elevation_classes()
allocate(qice_l(lsize_l,0:nEC))
allocate(frac_l(lsize_l,0:nEC))
allocate(tmp_field_l(lsize_l))
do ec = 0, nEC
elevclass_as_string = glc_elevclass_as_string(ec)
frac_field = Sg_frac_field // elevclass_as_string ! Sg_ice_covered01, etc.
call mct_aVect_exportRattr(g2x_lx, trim(frac_field), tmp_field_l)
frac_l(:,ec) = tmp_field_l(:)
qice_field = qice_fieldname // elevclass_as_string ! Flgl_qice01, etc.
call mct_aVect_exportRattr(l2gacc_lx(eli), trim(qice_field), tmp_field_l)
qice_l(:,ec) = tmp_field_l(:)
enddo
! clean the temporary attribute vector g2x_lx
call mct_aVect_clean(g2x_lx)
! Sum qice over local land grid cells
! initialize qice sum
local_accum_on_land_grid = 0.0_r8
local_ablat_on_land_grid = 0.0_r8
do n = 1, lsize_l
effective_area = min(lfrac(n),Sg_icemask_l(n)) * aream_l(n)
do ec = 0, nEC
if (qice_l(n,ec) >= 0.0_r8) then
local_accum_on_land_grid = local_accum_on_land_grid &
+ effective_area * frac_l(n,ec) * qice_l(n,ec)
else
local_ablat_on_land_grid = local_ablat_on_land_grid &
+ effective_area * frac_l(n,ec) * qice_l(n,ec)
endif
enddo ! ec
enddo ! n
call shr_mpi_sum(local_accum_on_land_grid, &
global_accum_on_land_grid, &
mpicom, 'accum_l')
call shr_mpi_sum(local_ablat_on_land_grid, &
global_ablat_on_land_grid, &
mpicom, 'ablat_l')
call shr_mpi_bcast(global_accum_on_land_grid, mpicom)
call shr_mpi_bcast(global_ablat_on_land_grid, mpicom)
! Sum qice_g over local glc grid cells.
! Note: This sum uses the coupler areas (aream_g), which differ from the native CISM areas.
! But since the original qice_g (from bilinear remapping) has been multiplied by
! area_g/aream_g above, this calculation is equivalent to multiplying the original qice_g
! by the native CISM areas (area_g).
! If Flgl_qice were changed to a state (and not included in seq_flds_x2g_fluxes),
! then it would be appropriate to use the native CISM areas in this sum.
! Export Sg_icemask from g2x_gx to a local array
allocate(Sg_icemask_g(lsize_g))
call mct_aVect_exportRattr(g2x_gx, Sg_icemask_field, Sg_icemask_g)
local_accum_on_glc_grid = 0.0_r8
local_ablat_on_glc_grid = 0.0_r8
do n = 1, lsize_g
if (qice_g(n) >= 0.0_r8) then
local_accum_on_glc_grid = local_accum_on_glc_grid &
+ Sg_icemask_g(n) * aream_g(n) * qice_g(n)
else
local_ablat_on_glc_grid = local_ablat_on_glc_grid &
+ Sg_icemask_g(n) * aream_g(n) * qice_g(n)
endif
enddo ! n
call shr_mpi_sum(local_accum_on_glc_grid, &
global_accum_on_glc_grid, &
mpicom, 'accum_g')
call shr_mpi_sum(local_ablat_on_glc_grid, &
global_ablat_on_glc_grid, &
mpicom, 'ablat_g')
call shr_mpi_bcast(global_accum_on_glc_grid, mpicom)
call shr_mpi_bcast(global_ablat_on_glc_grid, mpicom)
! Renormalize
if (global_accum_on_glc_grid > 0.0_r8) then
accum_renorm_factor = global_accum_on_land_grid / global_accum_on_glc_grid
else
accum_renorm_factor = 0.0_r8
endif
if (global_ablat_on_glc_grid < 0.0_r8) then ! negative by definition
ablat_renorm_factor = global_ablat_on_land_grid / global_ablat_on_glc_grid
else
ablat_renorm_factor = 0.0_r8
endif
if (iamroot) then
write(logunit,*) 'accum_renorm_factor = ', accum_renorm_factor
write(logunit,*) 'ablat_renorm_factor = ', ablat_renorm_factor
endif
do n = 1, lsize_g
if (qice_g(n) >= 0.0_r8) then
qice_g(n) = qice_g(n) * accum_renorm_factor
else
qice_g(n) = qice_g(n) * ablat_renorm_factor
endif
enddo
deallocate(aream_l)
deallocate(lfrac)
deallocate(Sg_icemask_l)
deallocate(Sg_icemask_g)
deallocate(tmp_field_l)
deallocate(qice_l)
deallocate(frac_l)
end subroutine prep_glc_renormalize_smb
!================================================================================================
function prep_glc_get_l2x_gx()
type(mct_aVect), pointer :: prep_glc_get_l2x_gx(:)
prep_glc_get_l2x_gx => l2x_gx(:)
end function prep_glc_get_l2x_gx
function prep_glc_get_l2gacc_lx()
type(mct_aVect), pointer :: prep_glc_get_l2gacc_lx(:)
prep_glc_get_l2gacc_lx => l2gacc_lx(:)
end function prep_glc_get_l2gacc_lx
function prep_glc_get_l2gacc_lx_one_instance(lnd_inst)
integer, intent(in) :: lnd_inst
type(mct_aVect), pointer :: prep_glc_get_l2gacc_lx_one_instance
prep_glc_get_l2gacc_lx_one_instance => l2gacc_lx(lnd_inst)
end function prep_glc_get_l2gacc_lx_one_instance
function prep_glc_get_l2gacc_lx_cnt()
integer, pointer :: prep_glc_get_l2gacc_lx_cnt
prep_glc_get_l2gacc_lx_cnt => l2gacc_lx_cnt
end function prep_glc_get_l2gacc_lx_cnt
function prep_glc_get_o2x_gx()
type(mct_aVect), pointer :: prep_glc_get_o2x_gx(:)
prep_glc_get_o2x_gx => o2x_gx(:)
end function prep_glc_get_o2x_gx
function prep_glc_get_x2gacc_gx()
type(mct_aVect), pointer :: prep_glc_get_x2gacc_gx(:)
prep_glc_get_x2gacc_gx => x2gacc_gx(:)
end function prep_glc_get_x2gacc_gx
function prep_glc_get_x2gacc_gx_cnt()
integer, pointer :: prep_glc_get_x2gacc_gx_cnt
prep_glc_get_x2gacc_gx_cnt => x2gacc_gx_cnt
end function prep_glc_get_x2gacc_gx_cnt
function prep_glc_get_mapper_Sl2g()
type(seq_map), pointer :: prep_glc_get_mapper_Sl2g
prep_glc_get_mapper_Sl2g => mapper_Sl2g
end function prep_glc_get_mapper_Sl2g
function prep_glc_get_mapper_Fl2g()
type(seq_map), pointer :: prep_glc_get_mapper_Fl2g
prep_glc_get_mapper_Fl2g => mapper_Fl2g
end function prep_glc_get_mapper_Fl2g
function prep_glc_get_mapper_So2g()
type(seq_map), pointer :: prep_glc_get_mapper_So2g
prep_glc_get_mapper_So2g=> mapper_So2g
end function prep_glc_get_mapper_So2g
function prep_glc_get_mapper_Fo2g()
type(seq_map), pointer :: prep_glc_get_mapper_Fo2g
prep_glc_get_mapper_Fo2g=> mapper_Fo2g
end function prep_glc_get_mapper_Fo2g
!***********************************************************************
!
! routine compute_melt_fluxes
!
!> \brief Computes ocean and ice melt fluxes, etc.
!> \author Xylar Asay-Davis
!> \date 3/27/2015
!> This routine computes melt fluxes (melt rate, temperature fluxes
!> into the ice and the ocean, and salt flux) as well as the interface
!> temperature and salinity. This routine expects an ice temperature
!> in the bottom layer of ice and ocean temperature and salinity in
!> the top ocean layer as well as the pressure at the ice/ocean interface.
!>
!> The ocean heat and salt transfer velocities are determined based on
!> observations of turbulent mixing rates in the under-ice boundary layer.
!> They should be the product of the friction velocity and a (possibly
!> spatially variable) non-dimenional transfer coefficient.
!>
!> The iceTemperatureDistance is the distance between the location
!> where the iceTemperature is supplied and the ice-ocean interface,
!> used to compute a temperature gradient. The ice thermal conductivity,
!> SHR_CONST_KAPPA_LAND_ICE, is zero for the freezing solution from Holland and Jenkins
!> (1999) in which the ice is purely insulating.
!
!-----------------------------------------------------------------------
subroutine compute_melt_fluxes( &
oceanTemperature, &
oceanSalinity, &
oceanHeatTransferVelocity, &
oceanSaltTransferVelocity, &
interfacePressure, &
iceTemperature, &
iceTemperatureDistance, &
iceFloatingMask, &
outInterfaceSalinity, &
outInterfaceTemperature, &
outFreshwaterFlux, &
outOceanHeatFlux, &
outIceHeatFlux, &
gsize)
use shr_const_mod, only: SHR_CONST_CPICE, &
SHR_CONST_CPSW, &
SHR_CONST_LATICE, &
SHR_CONST_RHOICE, &
SHR_CONST_RHOSW, &
SHR_CONST_DTF_DP, &
SHR_CONST_DTF_DS, &
SHR_CONST_DTF_DPDS, &
SHR_CONST_TF0, &
SHR_CONST_KAPPA_LAND_ICE
!-----------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------
real (kind=r8), dimension(:), intent(in) :: &
oceanTemperature, & !< Input: ocean temperature in top layer
oceanSalinity, & !< Input: ocean salinity in top layer
oceanHeatTransferVelocity, & !< Input: ocean heat transfer velocity
oceanSaltTransferVelocity, & !< Input: ocean salt transfer velocity
interfacePressure, & !< Input: pressure at the ice-ocean interface
iceTemperature, & !< Input: ice temperature in bottom layer
iceTemperatureDistance !< Input: distance to ice temperature from ice-ocean interface
integer, dimension(:), intent(in) :: &
iceFloatingMask !< Input: mask of cells that contain floating ice
integer, intent(in) :: gsize !< Input: number of values in each array
!-----------------------------------------------------------------
!
! output variables
!
!-----------------------------------------------------------------
real (kind=r8), dimension(:), intent(out) :: &
outInterfaceSalinity, & !< Output: ocean salinity at the interface
outInterfaceTemperature, & !< Output: ice/ocean temperature at the interface
outFreshwaterFlux, & !< Output: ocean thickness flux (melt rate)
outOceanHeatFlux, & !< Output: the temperature flux into the ocean
outIceHeatFlux !< Output: the temperature flux into the ice
!-----------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------
real (kind=r8) :: T0, transferVelocityRatio, Tlatent, nu, a, b, c, eta, &
iceHeatFluxCoeff, iceDeltaT, dTf_dS
integer :: n
character(*), parameter :: subname = '(compute_melt_fluxes)'
real (kind=r8), parameter :: minInterfaceSalinity = 0.001_r8
real (kind=r8), parameter :: referencePressure = 0.0_r8 ! Using reference pressure of 0
real (kind=r8) :: pressureOffset
Tlatent = SHR_CONST_LATICE/SHR_CONST_CPSW
do n = 1, gsize
if (iceFloatingMask(n) == 0) cycle ! Only calculate on floating cells
if (oceanHeatTransferVelocity(n) == 0.0_r8) then
write(logunit,*) 'compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0 at index ', n
call shr_sys_abort('compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0')
end if
iceHeatFluxCoeff = SHR_CONST_RHOICE*SHR_CONST_CPICE*SHR_CONST_KAPPA_LAND_ICE/iceTemperatureDistance(n)
nu = iceHeatFluxCoeff/(SHR_CONST_RHOSW*SHR_CONST_CPSW*oceanHeatTransferVelocity(n))
pressureOffset = max(interfacePressure(n) - referencePressure, 0.0_r8)
T0 = SHR_CONST_TF0 + SHR_CONST_DTF_DP * pressureOffset
!Note: These two terms for T0 are not needed because we are evaluating at salinity=0:
!+ SHR_CONST_DTF_DS * oceanSalinity(n) + SHR_CONST_DTF_DPDS * pressureOffset * oceanSalinity(n)
iceDeltaT = T0 - iceTemperature(n)
dTf_dS = SHR_CONST_DTF_DS + SHR_CONST_DTF_DPDS * pressureOffset
transferVelocityRatio = oceanSaltTransferVelocity(n)/oceanHeatTransferVelocity(n)
a = -1.0_r8 * dTf_dS * (1.0_r8 + nu)
b = transferVelocityRatio*Tlatent - nu*iceDeltaT + oceanTemperature(n) - T0
c = -transferVelocityRatio*Tlatent*max(oceanSalinity(n), 0.0_r8)
! a is non-negative; c is strictly non-positive so we never get imaginary roots.
! Since a can be zero, we need a solution of the quadratic equation for 1/Si instead of Si.
! Following: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
! Since a and -c are are non-negative, the term in the square root is also always >= |b|.
! In all reasonable cases, b will be strictly positive, since transferVelocityRatio*Tlatent ~ 2 C,
! T0 ~ -1.8 C and oceanTemperature should never be able to get below about -3 C
! As long as either b or both a and c are greater than zero, the strictly non-negative root is
outInterfaceSalinity(n) = max(-(2.0_r8*c)/(b + sqrt(b**2 - 4.0_r8*a*c)), minInterfaceSalinity)
outInterfaceTemperature(n) = dTf_dS*outInterfaceSalinity(n)+T0
outFreshwaterFlux(n) = SHR_CONST_RHOSW*oceanSaltTransferVelocity(n) &
* (oceanSalinity(n)/outInterfaceSalinity(n) - 1.0_r8)
! According to Jenkins et al. (2001), the temperature fluxes into the ocean are:
! 1. the advection of meltwater into the top layer (or removal for freezing)
! 2. the turbulent transfer of heat across the boundary layer, based on the termal driving
outOceanHeatFlux(n) = SHR_CONST_CPSW*(outFreshwaterFlux(n)*outInterfaceTemperature(n) &
- SHR_CONST_RHOSW*oceanHeatTransferVelocity(n)*(oceanTemperature(n)-outInterfaceTemperature(n)))
! the temperature fluxes into the ice are:
! 1. the advection of ice at the interface temperature out of the domain due to melting
! (or in due to freezing)
! 2. the diffusion (if any) of heat into the ice, based on temperature difference between
! the reference point in the ice (either the surface or the middle of the bottom layer)
! and the interface
outIceHeatFlux(n) = -SHR_CONST_CPICE*outFreshwaterFlux(n)*outInterfaceTemperature(n)
outIceHeatFlux(n) = outIceHeatFlux(n) &
- iceHeatFluxCoeff*(iceTemperature(n) - outInterfaceTemperature(n))
end do
!--------------------------------------------------------------------
end subroutine compute_melt_fluxes
end module prep_glc_mod
