!===============================================================================
! SVN $Id: seq_diag_mct.F90 19586 2009-11-23 21:50:15Z kauff $
! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/drv/seq_mct/branch_tags/cesm1_0_rel_tags/cesm1_0_rel01_drvseq3_1_32/driver/seq_diag_mct.F90 $
!===============================================================================
!BOP ===========================================================================
!
! !MODULE: seq_diag_mod -- computes spatial \& time averages of fluxed quatities
!
! !DESCRIPTION:
! The coupler is required to do certain diagnostics, those calculations are
! located in this module.
!
! !REMARKS:
! Sign convention:
! positive value <=> the model is gaining water, heat, momentum, etc.
! Unit convention:
! heat flux ~ W/m^2
! momentum flux ~ N/m^2
! water flux ~ (kg/s)/m^2
! salt flux ~ (kg/s)/m^2
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - updated budget implementation
! 2007-may-07 - B. Kauffman - initial port to cpl7.
! 2002-nov-21 - R. Jacob - initial port to cpl6.
! 199x-mmm-dd - B. Kauffman - original version in cpl4.
!
! !INTERFACE: ------------------------------------------------------------------
module seq_diag_mct 3,10
! !USES:
use shr_kind_mod
, only: r8 => shr_kind_r8, in=>shr_kind_in
use shr_kind_mod, only: i8 => shr_kind_i8, cl=>shr_kind_cl
use shr_sys_mod ! system calls
use shr_mpi_mod ! mpi wrappers
use shr_const_mod ! shared constants
use mct_mod ! mct wrappers
use esmf_mod
use seq_flds_indices ! shared field indicies
use seq_comm_mct ! mpi comm groups & related
use seq_cdata_mod
use seq_timemgr_mod
implicit none
private
! !PUBLIC TYPES:
! none
!PUBLIC MEMBER FUNCTIONS:
public seq_diag_zero_mct
public seq_diag_atm_mct
public seq_diag_lnd_mct
public seq_diag_rtm_mct
public seq_diag_ocn_mct
public seq_diag_ice_mct
public seq_diag_accum_mct
public seq_diag_sum0_mct
public seq_diag_print_mct
public seq_diag_avect_mct
public seq_diag_avdiff_mct
!EOP
!----------------------------------------------------------------------------
! Local data
!----------------------------------------------------------------------------
!----- local constants -----
real(r8),parameter :: HFLXtoWFLX = & ! water flux implied by latent heat of fusion
& - (shr_const_ocn_ref_sal-shr_const_ice_ref_sal) / &
& (shr_const_ocn_ref_sal*shr_const_latice)
!--- C for component ---
!--- "r" is recieve in the coupler, "s" is send from the coupler
integer(in),parameter :: c_size = 18
integer(in),parameter :: c_atm_as = 1 ! model index: atm
integer(in),parameter :: c_atm_ar = 2 ! model index: atm
integer(in),parameter :: c_inh_is = 3 ! model index: ice, northern
integer(in),parameter :: c_inh_ir = 4 ! model index: ice, northern
integer(in),parameter :: c_ish_is = 5 ! model index: ice, southern
integer(in),parameter :: c_ish_ir = 6 ! model index: ice, southern
integer(in),parameter :: c_lnd_ls = 7 ! model index: lnd
integer(in),parameter :: c_lnd_lr = 8 ! model index: lnd
integer(in),parameter :: c_ocn_os = 9 ! model index: ocn
integer(in),parameter :: c_ocn_or =10 ! model index: ocn
! --- on atm grid ---
integer(in),parameter :: c_inh_as =11 ! model index: ice, northern
integer(in),parameter :: c_inh_ar =12 ! model index: ice, northern
integer(in),parameter :: c_ish_as =13 ! model index: ice, southern
integer(in),parameter :: c_ish_ar =14 ! model index: ice, southern
integer(in),parameter :: c_lnd_as =15 ! model index: lnd
integer(in),parameter :: c_lnd_ar =16 ! model index: lnd
integer(in),parameter :: c_ocn_as =17 ! model index: ocn
integer(in),parameter :: c_ocn_ar =18 ! model index: ocn
character(len=8),parameter :: cname(c_size) = &
(/' c2a_atm',' a2c_atm',' c2i_inh',' i2c_inh',' c2i_ish',' i2c_ish', &
' c2l_lnd',' l2c_lnd',' c2o_ocn',' o2c_ocn', &
' c2a_inh',' a2c_inh',' c2a_ish',' a2c_ish', &
' c2a_lnd',' a2c_lnd',' c2a_ocn',' a2c_ocn' /)
!--- F for field ---
integer(in),parameter :: f_size = 17
integer(in),parameter :: f_a = 1 ! index for area
integer(in),parameter :: f_h = 2 ! 1st index for heat
integer(in),parameter :: f_w = 11 ! 1st index for water
integer(in),parameter :: f_area = 1 ! area (wrt to unit sphere)
integer(in),parameter :: f_hfrz = 2 ! heat : latent, freezing
integer(in),parameter :: f_hmelt = 3 ! heat : latent, melting
integer(in),parameter :: f_hswnet = 4 ! heat : short wave, net
integer(in),parameter :: f_hlwdn = 5 ! heat : longwave down
integer(in),parameter :: f_hlwup = 6 ! heat : longwave up
integer(in),parameter :: f_hlatv = 7 ! heat : latent, vaporization
integer(in),parameter :: f_hlatf = 8 ! heat : latent, fusion, snow
integer(in),parameter :: f_hioff = 9 ! heat : latent, fusion, frozen runoff
integer(in),parameter :: f_hsen =10 ! heat : sensible
integer(in),parameter :: f_wfrz =11 ! water: freezing
integer(in),parameter :: f_wmelt =12 ! water: melting
integer(in),parameter :: f_wrain =13 ! water: precip, liquid
integer(in),parameter :: f_wsnow =14 ! water: precip, frozen
integer(in),parameter :: f_wevap =15 ! water: evaporation
integer(in),parameter :: f_wroff =16 ! water: runoff
integer(in),parameter :: f_wioff =17 ! water: frozen runoff
character(len=8),parameter :: fname(f_size) = &
(/' area',' hfreeze',' hmelt',' hnetsw',' hlwdn', &
' hlwup',' hlatvap',' hlatfus',' hiroff',' hsen', &
' wfreeze',' wmelt',' wrain',' wsnow', &
' wevap',' wrunoff',' wfrzrof' /)
!--- P for period ---
integer(in),parameter :: p_size = 5
integer(in),parameter :: p_inst = 1
integer(in),parameter :: p_day = 2
integer(in),parameter :: p_mon = 3
integer(in),parameter :: p_ann = 4
integer(in),parameter :: p_inf = 5
character(len=8),parameter :: pname(p_size) = &
(/' inst',' daily',' monthly',' annual','all_time' /)
! !PUBLIC DATA MEMBERS
!--- time-averaged (annual?) global budge diagnostics ---
!--- note: call sum0 then save budg_dataG and budg_ns on restart from/to root pe ---
real(r8),public :: budg_dataL(f_size,c_size,p_size) ! local sum, valid on all pes
real(r8),public :: budg_dataG(f_size,c_size,p_size) ! global sum, valid only on root pe
real(r8),public :: budg_ns (f_size,c_size,p_size) ! counter, valid only on root pe
character(len=*),parameter :: afldname = 'aream'
character(len=*),parameter :: latname = 'lat'
character(len=*),parameter :: afracname = 'afrac'
character(len=*),parameter :: lfracname = 'lfrac'
character(len=*),parameter :: ofracname = 'ofrac'
character(len=*),parameter :: ifracname = 'ifrac'
character(len=*),parameter :: ascalname = 'ascale'
character(len=*),parameter :: lfrinname = 'lfrin'
character(*),parameter :: modName = "(seq_diag_mct) "
integer(in),parameter :: debug = 0 ! internal debug level
!===============================================================================
contains
!===============================================================================
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_zero_mct - zero out global budget diagnostic data.
!
! !DESCRIPTION:
! Zero out global budget diagnostic data.
!
! !REVISION HISTORY:
! 2008-jul-11 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_zero_mct(EClock,mode) 2,3
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock), intent(in),optional :: EClock
character(len=*), intent(in),optional :: mode
!EOP
integer(IN) :: ip,yr,mon,day,sec
!----- formats -----
character(*),parameter :: subName = '(seq_diag_zero_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
if (.not. present(EClock) .and. .not. present(mode)) then
call shr_sys_abort(subName//' ERROR EClock or mode should be present')
endif
if (present(EClock)) then
call seq_timemgr_EClockGetData(EClock,curr_yr=yr, &
curr_mon=mon,curr_day=day,curr_tod=sec)
do ip = 1,p_size
if (ip == p_inst) then
budg_dataL(:,:,ip) = 0.0_r8
budg_dataG(:,:,ip) = 0.0_r8
budg_ns(:,:,ip) = 0.0_r8
endif
if (ip==p_day .and. sec==0) then
budg_dataL(:,:,ip) = 0.0_r8
budg_dataG(:,:,ip) = 0.0_r8
budg_ns(:,:,ip) = 0.0_r8
endif
if (ip==p_mon .and. day==1 .and. sec==0) then
budg_dataL(:,:,ip) = 0.0_r8
budg_dataG(:,:,ip) = 0.0_r8
budg_ns(:,:,ip) = 0.0_r8
endif
if (ip==p_ann .and. mon==1 .and. day==1 .and. sec==0) then
budg_dataL(:,:,ip) = 0.0_r8
budg_dataG(:,:,ip) = 0.0_r8
budg_ns(:,:,ip) = 0.0_r8
endif
enddo
endif
if (present(mode)) then
if (trim(mode) == 'inst') then
budg_dataL(:,:,p_inst) = 0.0_r8
budg_dataG(:,:,p_inst) = 0.0_r8
budg_ns(:,:,p_inst) = 0.0_r8
elseif (trim(mode) == 'day') then
budg_dataL(:,:,p_day) = 0.0_r8
budg_dataG(:,:,p_day) = 0.0_r8
budg_ns(:,:,p_day) = 0.0_r8
elseif (trim(mode) == 'mon') then
budg_dataL(:,:,p_mon) = 0.0_r8
budg_dataG(:,:,p_mon) = 0.0_r8
budg_ns(:,:,p_mon) = 0.0_r8
elseif (trim(mode) == 'ann') then
budg_dataL(:,:,p_ann) = 0.0_r8
budg_dataG(:,:,p_ann) = 0.0_r8
budg_ns(:,:,p_ann) = 0.0_r8
elseif (trim(mode) == 'inf') then
budg_dataL(:,:,p_inf) = 0.0_r8
budg_dataG(:,:,p_inf) = 0.0_r8
budg_ns(:,:,p_inf) = 0.0_r8
elseif (trim(mode) == 'all') then
budg_dataL(:,:,:) = 0.0_r8
budg_dataG(:,:,:) = 0.0_r8
budg_ns(:,:,:) = 0.0_r8
else
call shr_sys_abort(subname//' ERROR in mode '//trim(mode))
endif
endif
end subroutine seq_diag_zero_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_accum_mct - accum out global budget diagnostic data.
!
! !DESCRIPTION:
! Accum out global budget diagnostic data.
!
! !REVISION HISTORY:
! 2008-jul-11 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_accum_mct() 1
! !INPUT/OUTPUT PARAMETERS:
!EOP
integer(in) :: ip
!----- formats -----
character(*),parameter :: subName = '(seq_diag_accum_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
do ip = p_inst+1,p_size
budg_dataL(:,:,ip) = budg_dataL(:,:,ip) + budg_dataL(:,:,p_inst)
enddo
budg_ns(:,:,:) = budg_ns(:,:,:) + 1.0_r8
end subroutine seq_diag_accum_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_sum0_mct - sum local to global on root
!
! !DESCRIPTION:
! Sum local values to global on root
!
! !REVISION HISTORY:
! 2008-jul-19 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_sum0_mct() 1,1
! !INPUT/OUTPUT PARAMETERS:
!EOP
real(r8) :: budg_dataGtmp(f_size,c_size,p_size) ! temporary sum
integer(in) :: mpicom ! mpi comm
!----- formats -----
character(*),parameter :: subName = '(seq_diag_sum0_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
call seq_comm_setptrs(CPLID,mpicom=mpicom)
budg_dataGtmp = 0.0_r8
call shr_mpi_sum(budg_dataL,budg_dataGtmp,mpicom,subName)
budg_dataG = budg_dataG + budg_dataGtmp
budg_dataL = 0.0_r8
end subroutine seq_diag_sum0_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_atm_mct - compute global atm input/output flux diagnostics
!
! !DESCRIPTION:
! Compute global atm input/output flux diagnostics
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_atm_mct( dom_a, frac_a, a2x_a, x2a_a ) 1,1
! !INPUT/OUTPUT PARAMETERS:
type(mct_gGrid),intent(in) :: dom_a ! model domain
type(mct_aVect),intent(in) :: frac_a ! domain fractions
type(mct_aVect),intent(in),optional :: a2x_a ! model to drv bundle
type(mct_aVect),intent(in),optional :: x2a_a ! drv to model bundle
!EOP
!----- local -----
integer(in) :: k,n,ic,if,ip ! generic index
integer(in) :: kArea ! index of area field in aVect
integer(in) :: kLat ! index of lat field in aVect
integer(in) :: kl,ka,ko,ki ! fraction indices
integer(in) :: lSize ! size of aVect
real(r8) :: da,di,do,dl ! area of a grid cell
!----- formats -----
character(*),parameter :: subName = '(seq_diag_atm_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
if (.not. present(a2x_a) .and. .not. present(x2a_a)) then
call shr_sys_abort(subName//"ERROR: must input a bundle")
end if
kArea = mct_aVect_indexRA(dom_a%data,afldname)
kLat = mct_aVect_indexRA(dom_a%data,latname)
ka = mct_aVect_indexRA(frac_a,afracname)
kl = mct_aVect_indexRA(frac_a,lfracname)
ko = mct_aVect_indexRA(frac_a,ofracname)
ki = mct_aVect_indexRA(frac_a,ifracname)
!---------------------------------------------------------------------------
! add values found in this bundle to the budget table
!---------------------------------------------------------------------------
ip = p_inst
if (present(a2x_a)) then
lSize = mct_avect_lSize(a2x_a)
do n=1,lSize
do k=1,4
if (k == 1) then
ic = c_atm_ar
da = -dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ka,n)
elseif (k == 2) then
ic = c_lnd_ar
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(kl,n)
elseif (k == 3) then
ic = c_ocn_ar
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ko,n)
elseif (k == 4) then
if (dom_a%data%rAttr(kLat,n) > 0.0_r8) then
ic = c_inh_ar
else
ic = c_ish_ar
endif
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ki,n)
endif
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da
if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_swnet,n)
if = f_hlwdn ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_lwdn,n)
if = f_wrain ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_rainc,n) &
+ da*a2x_a%rAttr(index_a2x_Faxa_rainl,n)
if = f_wsnow ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_snowc,n) &
+ da*a2x_a%rAttr(index_a2x_Faxa_snowl,n)
enddo
enddo
! --- heat implied by snow flux ---
ic = c_atm_ar; budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
ic = c_lnd_ar; budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
ic = c_ocn_ar; budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
ic = c_inh_ar; budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
ic = c_ish_ar; budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
end if
if (present(x2a_a)) then
lSize = mct_avect_lSize(x2a_a)
do n=1,lSize
do k=1,4
if (k == 1) then
ic = c_atm_as
da = -dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ka,n)
elseif (k == 2) then
ic = c_lnd_as
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(kl,n)
elseif (k == 3) then
ic = c_ocn_as
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ko,n)
elseif (k == 4) then
if (dom_a%data%rAttr(kLat,n) > 0.0_r8) then
ic = c_inh_as
else
ic = c_ish_as
endif
da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ki,n)
endif
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da
if = f_hlwup; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_lwup,n)
if = f_hlatv; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_lat,n)
if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_sen,n)
if = f_wevap; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_evap,n)
enddo
enddo
end if
end subroutine seq_diag_atm_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_lnd_mct - compute global lnd input/output flux diagnostics
!
! !DESCRIPTION:
! Compute global lnd input/output flux diagnostics
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_lnd_mct( dom_l, frac_l, l2x_l, x2l_l) 1,1
type(mct_gGrid),intent(in) :: dom_l ! model domain
type(mct_aVect),intent(in) :: frac_l ! frac bundle
type(mct_aVect),intent(in),optional :: l2x_l ! model to drv bundle
type(mct_aVect),intent(in),optional :: x2l_l ! drv to model bundle
!EOP
!----- local -----
integer(in) :: k,n,ic,if,ip ! generic index
integer(in) :: kArea ! index of area field in aVect
integer(in) :: kLat ! index of lat field in aVect
integer(in) :: kl,ka,ko,ki ! fraction indices
integer(in) :: lSize ! size of aVect
real(r8) :: da,di,do,dl ! area of a grid cell
!----- formats -----
character(*),parameter :: subName = '(seq_diag_lnd_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
if (.not. present(l2x_l) .and. .not. present(x2l_l)) then
call shr_sys_abort(subName//"ERROR: must input a bundle")
end if
!---------------------------------------------------------------------------
! add values found in this bundle to the budget table
!---------------------------------------------------------------------------
ip = p_inst
kArea = mct_aVect_indexRA(dom_l%data,afldname)
kl = mct_aVect_indexRA(frac_l,lfrinname)
ka = mct_aVect_indexRA(dom_l%data,ascalname)
if (present(l2x_l)) then
lSize = mct_avect_lSize(l2x_l)
ic = c_lnd_lr
do n=1,lSize
dl = dom_l%data%rAttr(kArea,n) * frac_l%rAttr(kl,n) * dom_l%data%rAttr(ka,n)
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl
if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_swnet,n)
if = f_hlwup ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_lwup,n)
if = f_hlatv ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_lat,n)
if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_sen,n)
if = f_wevap ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_evap,n)
end do
end if
if (present(x2l_l)) then
lSize = mct_avect_lSize(x2l_l)
ic = c_lnd_ls
do n=1,lSize
dl = dom_l%data%rAttr(kArea,n) * frac_l%rAttr(kl,n) * dom_l%data%rAttr(ka,n)
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl
if = f_hlwdn; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_lwdn,n)
if = f_wrain; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_rainc,n) &
+ dl*x2l_l%rAttr(index_x2l_Faxa_rainl,n)
if = f_wsnow; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_snowc,n) &
+ dl*x2l_l%rAttr(index_x2l_Faxa_snowl,n)
end do
budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
end if
end subroutine seq_diag_lnd_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_rtm_mct - compute global rtm input/output flux diagnostics
!
! !DESCRIPTION:
! Compute global rtm input/output flux diagnostics
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_rtm_mct( dom_r, r2x_r) 1
type(mct_gGrid),intent(in) :: dom_r ! model domain
type(mct_aVect),intent(in) :: r2x_r ! model to drv bundle
!EOP
!----- local -----
integer(in) :: k,n,ic,if,ip ! generic index
integer(in) :: kArea ! index of area field in aVect
integer(in) :: kLat ! index of lat field in aVect
integer(in) :: kl,ka,ko,ki ! fraction indices
integer(in) :: lSize ! size of aVect
real(r8) :: da,di,do,dl ! area of a grid cell
!----- formats -----
character(*),parameter :: subName = '(seq_diag_rtm_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!---------------------------------------------------------------------------
! add values found in this bundle to the budget table
!---------------------------------------------------------------------------
ip = p_inst
ic = c_lnd_lr
kArea = mct_aVect_indexRA(dom_r%data,afldname)
lSize = mct_avect_lSize(r2x_r)
do n=1,lSize
dl = dom_r%data%rAttr(kArea,n)
if = f_wroff; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - dl*r2x_r%rAttr(index_r2x_Forr_roff,n)
if = f_wioff; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - dl*r2x_r%rAttr(index_r2x_Forr_ioff,n)
end do
budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice
end subroutine seq_diag_rtm_mct
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_ocn_mct - compute global ocn input/output flux diagnostics
!
! !DESCRIPTION:
! Compute global ocn input/output flux diagnostics
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_ocn_mct( dom_o, frac_o, o2x_o, x2o_o, xao_o, r2x_o) 1,1
type(mct_gGrid),intent(in) :: dom_o ! model domain
type(mct_aVect),intent(in) :: frac_o ! frac bundle
type(mct_aVect),intent(in),optional :: o2x_o ! model to drv bundle
type(mct_aVect),intent(in),optional :: x2o_o ! drv to model bundle
type(mct_aVect),intent(in),optional :: xao_o ! drv to model bundle
type(mct_aVect),intent(in),optional :: r2x_o ! roff to drv bundle
!EOP
!----- local -----
integer(in) :: k,n,if,ic,ip ! generic index
integer(in) :: kArea ! index of area field in aVect
integer(in) :: kLat ! index of lat field in aVect
integer(in) :: kl,ka,ko,ki ! fraction indices
integer(in) :: lSize ! size of aVect
real(r8) :: da,di,do,dl ! area of a grid cell
!----- formats -----
character(*),parameter :: subName = '(seq_diag_ocn_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
if (.not. present(o2x_o) .and. .not. present(x2o_o) .and. .not. present(xao_o)) then
call shr_sys_abort(subName//"ERROR: must input a bundle")
end if
!---------------------------------------------------------------------------
! add values found in this bundle to the budget table
!---------------------------------------------------------------------------
ip = p_inst
kArea = mct_aVect_indexRA(dom_o%data,afldname)
ko = mct_aVect_indexRA(frac_o,ofracname)
ki = mct_aVect_indexRA(frac_o,ifracname)
if (present(o2x_o)) then
lSize = mct_avect_lSize(o2x_o)
ic = c_ocn_or
do n=1,lSize
do = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
di = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
if = f_area; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do
if = f_hfrz; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*max(0.0_r8,o2x_o%rAttr(index_o2x_Fioo_q,n))
end do
budg_dataL(f_wfrz,ic,ip) = budg_dataL(f_hfrz,ic,ip) * HFLXtoWFLX
end if
if (present(xao_o)) then
lSize = mct_avect_lSize(xao_o)
ic = c_ocn_or
do n=1,lSize
do = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
if = f_hlwup; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_lwup,n)
if = f_hlatv; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_lat,n)
if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_sen,n)
if = f_wevap; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_evap,n)
end do
end if
if (present(x2o_o)) then
lSize = mct_avect_lSize(x2o_o)
ic = c_ocn_os
do n=1,lSize
do = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
di = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do
if = f_hmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_melth,n)
if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_swnet,n)
if = f_hlwdn ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_lwdn,n)
if = f_wmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_meltw,n)
if = f_wrain ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_rain,n)
if = f_wsnow ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_snow,n)
! if = f_wroff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Forr_roff,n)
! if = f_wioff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Forr_ioff,n)
end do
budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
end if
if (present(r2x_o)) then
lSize = mct_avect_lSize(r2x_o)
ic = c_ocn_os
do n=1,lSize
do = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
di = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
if = f_wroff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*r2x_o%rAttr(index_r2x_Forr_roff,n)
if = f_wioff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*r2x_o%rAttr(index_r2x_Forr_ioff,n)
end do
budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice
end if
end subroutine seq_diag_ocn_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_ice_mct - compute global ice input/output flux diagnostics
!
! !DESCRIPTION:
! Compute global ice input/output flux diagnostics
!
! !REVISION HISTORY:
! 2008-jul-10 - T. Craig - update
!
! !INTERFACE: ------------------------------------------------------------------
subroutine seq_diag_ice_mct( dom_i, frac_i, i2x_i, x2i_i) 2,1
type(mct_gGrid),intent(in) :: dom_i ! model domain
type(mct_aVect),intent(in) :: frac_i ! frac bundle
type(mct_aVect),intent(in),optional :: i2x_i ! model to drv bundle
type(mct_aVect),intent(in),optional :: x2i_i ! drv to model bundle
!EOP
!----- local -----
integer(in) :: k,n,ic,if,ip ! generic index
integer(in) :: kArea ! index of area field in aVect
integer(in) :: kLat ! index of lat field in aVect
integer(in) :: kl,ka,ko,ki ! fraction indices
integer(in) :: lSize ! size of aVect
real(r8) :: da,di,do,dl ! area of a grid cell
!----- formats -----
character(*),parameter :: subName = '(seq_diag_ice_mct) '
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
if (.not. present(i2x_i) .and. .not. present(x2i_i)) then
call shr_sys_abort(subName//"ERROR: must input a bundle")
end if
!---------------------------------------------------------------------------
! add values found in this bundle to the budget table
!---------------------------------------------------------------------------
ip = p_inst
kArea = mct_aVect_indexRA(dom_i%data,afldname)
kLat = mct_aVect_indexRA(dom_i%data,latname)
ki = mct_aVect_indexRA(frac_i,ifracname)
ko = mct_aVect_indexRA(frac_i,ofracname)
if (present(i2x_i)) then
lSize = mct_avect_lSize(i2x_i)
do n=1,lSize
if (dom_i%data%rAttr(kLat,n) > 0.0_r8) then
ic = c_inh_ir
else
ic = c_ish_ir
endif
di = dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ki,n)
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di
if = f_hmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - di*i2x_i%rAttr(index_i2x_Fioi_melth,n)
if = f_wmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - di*i2x_i%rAttr(index_i2x_Fioi_meltw,n)
if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_swnet,n) &
- di*i2x_i%rAttr(index_i2x_Fioi_swpen,n)
if = f_hlwup ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_lwup,n)
if = f_hlatv ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_lat,n)
if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_sen,n)
if = f_wevap ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_evap,n)
end do
end if
if (present(x2i_i)) then
lSize = mct_avect_lSize(x2i_i)
do n=1,lSize
if (dom_i%data%rAttr(kLat,n) > 0.0_r8) then
ic = c_inh_is
else
ic = c_ish_is
endif
do = dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ko,n)
di = dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ki,n)
if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di
if = f_hlwdn; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_lwdn,n)
if = f_wrain; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_rain,n)
if = f_wsnow; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_snow,n)
if = f_hfrz ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - (do+di)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_q,n))
end do
ic = c_inh_is
budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
budg_dataL(f_wfrz ,ic,ip) = budg_dataL(f_hfrz ,ic,ip)*HFLXtoWFLX
ic = c_ish_is
budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
budg_dataL(f_wfrz ,ic,ip) = budg_dataL(f_hfrz ,ic,ip)*HFLXtoWFLX
end if
end subroutine seq_diag_ice_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_print_mct - print global budget diagnostics
!
! !DESCRIPTION:
! Print global budget diagnostics.
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
SUBROUTINE seq_diag_print_mct(EClock,stop_alarm, & 1,5
budg_print_inst, budg_print_daily, budg_print_month, &
budg_print_ann, budg_print_ltann, budg_print_ltend)
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock), intent(in) :: EClock
logical,intent(in) :: stop_alarm
integer,intent(in) :: budg_print_inst
integer,intent(in) :: budg_print_daily
integer,intent(in) :: budg_print_month
integer,intent(in) :: budg_print_ann
integer,intent(in) :: budg_print_ltann
integer,intent(in) :: budg_print_ltend
!EOP
!--- local ---
integer(in) :: ic,if,ip ! data array indicies
integer(in) :: ica,icl,icn,ics,ico
integer(in) :: icar,icxs,icxr,icas
integer(in) :: n ! loop counter
integer(in) :: nday ! number of days in time avg
integer(in) :: cdate,sec ! coded date, seconds
integer(in) :: yr,mon,day ! date
integer(in) :: iam ! pe number
integer(in) :: plev ! print level
logical :: sumdone ! has a sum been computed yet
character(len=40):: str ! string
real(r8) :: dataGpr (f_size,c_size,p_size) ! values to print, scaled and such
!----- formats -----
character(*),parameter :: subName = '(seq_diag_print_mct) '
character(*),parameter :: F00 = "('(seq_diag_print_mct) ',4a)"
!----- formats -----
character(*),parameter :: FAH="(4a,i9,i6)"
character(*),parameter :: FA0="(' ',8x,6(6x,a8,1x))"
character(*),parameter :: FA1="(' ',a8,6f15.8)"
!-------------------------------------------------------------------------------
! print instantaneous budget data
!-------------------------------------------------------------------------------
sumdone = .false.
call seq_comm_setptrs(CPLID,iam=iam)
call seq_timemgr_EClockGetData(EClock,curr_yr=yr, &
curr_mon=mon,curr_day=day,curr_tod=sec)
cdate = yr*10000+mon*100+day
do ip = 1,p_size
plev = 0
if (ip == p_inst) then
plev = max(plev,budg_print_inst)
endif
if (ip==p_day .and. sec==0) then
plev = max(plev,budg_print_daily)
endif
if (ip==p_mon .and. day==1 .and. sec==0) then
plev = max(plev,budg_print_month)
endif
if (ip==p_ann .and. mon==1 .and. day==1 .and. sec==0) then
plev = max(plev,budg_print_ann)
endif
if (ip==p_inf .and. mon==1 .and. day==1 .and. sec==0) then
plev = max(plev,budg_print_ltann)
endif
if (ip==p_inf .and. stop_alarm) then
plev = max(plev,budg_print_ltend)
endif
if (plev > 0) then
! ---- doprint ---- doprint ---- doprint ----
if (.not.sumdone) then
call seq_diag_sum0_mct()
dataGpr = budg_dataG
sumdone = .true.
! old budget normalizations (global area and 1e6 for water)
dataGpr = dataGpr/(4.0_r8*shr_const_pi)
dataGpr(f_w:f_size,:,:) = dataGpr(f_w:f_size,:,:) * 1.0e6_r8
dataGpr = dataGpr/budg_ns
if (iam /= 0) return
endif
! ---------------------------------------------------------
! ---- detail atm budgets and breakdown into components ---
! ---------------------------------------------------------
if (plev >= 3) then
do ic = 1,2
if (ic == 1) then
ica = c_atm_ar
icl = c_lnd_ar
icn = c_inh_ar
ics = c_ish_ar
ico = c_ocn_ar
str = "ATM_to_CPL"
elseif (ic == 2) then
ica = c_atm_as
icl = c_lnd_as
icn = c_inh_as
ics = c_ish_as
ico = c_ocn_as
str = "CPL_TO_ATM"
else
call shr_sys_abort(subname//' ERROR in ic index code 411')
endif
write(logunit,*) ' '
write(logunit,FAH) subname,trim(str)//' AREA BUDGET (m2/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM* '
do if = f_a, f_h-1
write(logunit,FA1) fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip)
enddo
write(logunit,*) ' '
write(logunit,FAH) subname,trim(str)//' HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM* '
do if = f_h, f_w-1
write(logunit,FA1) fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip)
enddo
write(logunit,FA1) ' *SUM*',sum(dataGpr(f_h:f_w-1,ica,ip)),sum(dataGpr(f_h:f_w-1,icl,ip)), &
sum(dataGpr(f_h:f_w-1,icn,ip)),sum(dataGpr(f_h:f_w-1,ics,ip)),sum(dataGpr(f_h:f_w-1,ico,ip)), &
sum(dataGpr(f_h:f_w-1,ica,ip))+sum(dataGpr(f_h:f_w-1,icl,ip))+ &
sum(dataGpr(f_h:f_w-1,icn,ip))+sum(dataGpr(f_h:f_w-1,ics,ip))+sum(dataGpr(f_h:f_w-1,ico,ip))
write(logunit,*) ' '
write(logunit,FAH) subname,trim(str)//' WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM* '
do if = f_w, f_size
write(logunit,FA1) fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip)
enddo
write(logunit,FA1) ' *SUM*',sum(dataGpr(f_w:f_size,ica,ip)),sum(dataGpr(f_w:f_size,icl,ip)), &
sum(dataGpr(f_w:f_size,icn,ip)),sum(dataGpr(f_w:f_size,ics,ip)),sum(dataGpr(f_w:f_size,ico,ip)), &
sum(dataGpr(f_w:f_size,ica,ip))+sum(dataGpr(f_w:f_size,icl,ip))+ &
sum(dataGpr(f_w:f_size,icn,ip))+sum(dataGpr(f_w:f_size,ics,ip))+sum(dataGpr(f_w:f_size,ico,ip))
enddo
endif ! plev
! ---------------------------------------------------------
! ---- detail lnd/ocn/ice component budgets ----
! ---------------------------------------------------------
if (plev >= 2) then
do ic = 1,4
if (ic == 1) then
icar = c_lnd_ar
icxs = c_lnd_ls
icxr = c_lnd_lr
icas = c_lnd_as
str = "LND"
elseif (ic == 2) then
icar = c_ocn_ar
icxs = c_ocn_os
icxr = c_ocn_or
icas = c_ocn_as
str = "OCN"
elseif (ic == 3) then
icar = c_inh_ar
icxs = c_inh_is
icxr = c_inh_ir
icas = c_inh_as
str = "ICE_NH"
elseif (ic == 4) then
icar = c_ish_ar
icxs = c_ish_is
icxr = c_ish_ir
icas = c_ish_as
str = "ICE_SH"
else
call shr_sys_abort(subname//' ERROR in ic index code 412')
endif
write(logunit,*) ' '
write(logunit,FAH) subname,trim(str)//' HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) cname(icar),cname(icxs),cname(icxr),cname(icas),' *SUM* '
do if = f_h, f_w-1
write(logunit,FA1) fname(if),-dataGpr(if,icar,ip),dataGpr(if,icxs,ip), &
dataGpr(if,icxr,ip),-dataGpr(if,icas,ip), &
-dataGpr(if,icar,ip)+dataGpr(if,icxs,ip)+ &
dataGpr(if,icxr,ip)-dataGpr(if,icas,ip)
enddo
write(logunit,FA1) ' *SUM*',-sum(dataGpr(f_h:f_w-1,icar,ip)),sum(dataGpr(f_h:f_w-1,icxs,ip)), &
sum(dataGpr(f_h:f_w-1,icxr,ip)),-sum(dataGpr(f_h:f_w-1,icas,ip)), &
-sum(dataGpr(f_h:f_w-1,icar,ip))+sum(dataGpr(f_h:f_w-1,icxs,ip))+ &
sum(dataGpr(f_h:f_w-1,icxr,ip))-sum(dataGpr(f_h:f_w-1,icas,ip))
write(logunit,*) ' '
write(logunit,FAH) subname,trim(str)//' WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) cname(icar),cname(icxs),cname(icxr),cname(icas),' *SUM* '
do if = f_w, f_size
write(logunit,FA1) fname(if),-dataGpr(if,icar,ip),dataGpr(if,icxs,ip), &
dataGpr(if,icxr,ip),-dataGpr(if,icas,ip), &
-dataGpr(if,icar,ip)+dataGpr(if,icxs,ip)+ &
dataGpr(if,icxr,ip)-dataGpr(if,icas,ip)
enddo
write(logunit,FA1) ' *SUM*',-sum(dataGpr(f_w:f_size,icar,ip)),sum(dataGpr(f_w:f_size,icxs,ip)), &
sum(dataGpr(f_w:f_size,icxr,ip)),-sum(dataGpr(f_w:f_size,icas,ip)), &
-sum(dataGpr(f_w:f_size,icar,ip))+sum(dataGpr(f_w:f_size,icxs,ip))+ &
sum(dataGpr(f_w:f_size,icxr,ip))-sum(dataGpr(f_w:f_size,icas,ip))
enddo
endif ! plev
! ---------------------------------------------------------
! ---- net summary budgets ----
! ---------------------------------------------------------
if (plev >= 1) then
write(logunit,*) ' '
write(logunit,FAH) subname,'NET AREA BUDGET (m2/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) ' atm',' lnd',' ocn',' ice nh',' ice sh',' *SUM* '
do if = 1,f_h-1
write(logunit,FA1) fname(if),dataGpr(if,c_atm_ar,ip), &
dataGpr(if,c_lnd_lr,ip), &
dataGpr(if,c_ocn_or,ip), &
dataGpr(if,c_inh_ir,ip), &
dataGpr(if,c_ish_ir,ip), &
dataGpr(if,c_atm_ar,ip)+ &
dataGpr(if,c_lnd_lr,ip)+ &
dataGpr(if,c_ocn_or,ip)+ &
dataGpr(if,c_inh_ir,ip)+ &
dataGpr(if,c_ish_ir,ip)
enddo
write(logunit,*) ' '
write(logunit,FAH) subname,'NET HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) ' atm',' lnd',' ocn',' ice nh',' ice sh',' *SUM* '
do if = f_h, f_w-1
write(logunit,FA1) fname(if),dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip), &
dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip), &
dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip), &
dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip), &
dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip), &
dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip)+ &
dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip)+ &
dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip)+ &
dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip)+ &
dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip)
enddo
write(logunit,FA1) ' *SUM*',sum(dataGpr(f_h:f_w-1,c_atm_ar,ip))+sum(dataGpr(f_h:f_w-1,c_atm_as,ip)), &
sum(dataGpr(f_h:f_w-1,c_lnd_lr,ip))+sum(dataGpr(f_h:f_w-1,c_lnd_ls,ip)), &
sum(dataGpr(f_h:f_w-1,c_ocn_or,ip))+sum(dataGpr(f_h:f_w-1,c_ocn_os,ip)), &
sum(dataGpr(f_h:f_w-1,c_inh_ir,ip))+sum(dataGpr(f_h:f_w-1,c_inh_is,ip)), &
sum(dataGpr(f_h:f_w-1,c_ish_ir,ip))+sum(dataGpr(f_h:f_w-1,c_ish_is,ip)), &
sum(dataGpr(f_h:f_w-1,c_atm_ar,ip))+sum(dataGpr(f_h:f_w-1,c_atm_as,ip))+ &
sum(dataGpr(f_h:f_w-1,c_lnd_lr,ip))+sum(dataGpr(f_h:f_w-1,c_lnd_ls,ip))+ &
sum(dataGpr(f_h:f_w-1,c_ocn_or,ip))+sum(dataGpr(f_h:f_w-1,c_ocn_os,ip))+ &
sum(dataGpr(f_h:f_w-1,c_inh_ir,ip))+sum(dataGpr(f_h:f_w-1,c_inh_is,ip))+ &
sum(dataGpr(f_h:f_w-1,c_ish_ir,ip))+sum(dataGpr(f_h:f_w-1,c_ish_is,ip))
write(logunit,*) ' '
write(logunit,FAH) subname,'NET WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
write(logunit,FA0) ' atm',' lnd',' ocn',' ice nh',' ice sh',' *SUM* '
do if = f_w, f_size
write(logunit,FA1) fname(if),dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip), &
dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip), &
dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip), &
dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip), &
dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip), &
dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip)+ &
dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip)+ &
dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip)+ &
dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip)+ &
dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip)
enddo
write(logunit,FA1) ' *SUM*',sum(dataGpr(f_w:f_size,c_atm_ar,ip))+sum(dataGpr(f_w:f_size,c_atm_as,ip)), &
sum(dataGpr(f_w:f_size,c_lnd_lr,ip))+sum(dataGpr(f_w:f_size,c_lnd_ls,ip)), &
sum(dataGpr(f_w:f_size,c_ocn_or,ip))+sum(dataGpr(f_w:f_size,c_ocn_os,ip)), &
sum(dataGpr(f_w:f_size,c_inh_ir,ip))+sum(dataGpr(f_w:f_size,c_inh_is,ip)), &
sum(dataGpr(f_w:f_size,c_ish_ir,ip))+sum(dataGpr(f_w:f_size,c_ish_is,ip)), &
sum(dataGpr(f_w:f_size,c_atm_ar,ip))+sum(dataGpr(f_w:f_size,c_atm_as,ip))+ &
sum(dataGpr(f_w:f_size,c_lnd_lr,ip))+sum(dataGpr(f_w:f_size,c_lnd_ls,ip))+ &
sum(dataGpr(f_w:f_size,c_ocn_or,ip))+sum(dataGpr(f_w:f_size,c_ocn_os,ip))+ &
sum(dataGpr(f_w:f_size,c_inh_ir,ip))+sum(dataGpr(f_w:f_size,c_inh_is,ip))+ &
sum(dataGpr(f_w:f_size,c_ish_ir,ip))+sum(dataGpr(f_w:f_size,c_ish_is,ip))
endif
write(logunit,*) ' '
! ---- doprint ---- doprint ---- doprint ----
endif ! plev > 0
enddo ! ip = 1,p_size
end subroutine seq_diag_print_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_avect_mct - print global budget diagnostics
!
! !DESCRIPTION:
! Print global diagnostics for AV/ID.
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
SUBROUTINE seq_diag_avect_mct(cdata,AV,comment) 23,6
use seq_infodata_mod
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(seq_cdata) , intent(in) :: cdata
type(mct_aVect) , intent(in) :: AV
character(len=*), intent(in), optional :: comment
!EOP
!--- local ---
type(mct_gGrid) ,pointer :: dom
type(seq_infodata_type),pointer :: infodata
type(mct_gsMap) ,pointer :: gsmap
integer(in) :: ID
logical :: bfbflag
integer(in) :: n,k ! counters
integer(in) :: npts,nptsg ! number of local/global pts in AV
integer(in) :: kflds ! number of fields in AV
real(r8),pointer :: sumbuf (:) ! sum buffer
real(r8),pointer :: minbuf (:) ! min buffer
real(r8),pointer :: maxbuf (:) ! max buffer
real(r8),pointer :: sumbufg(:) ! sum buffer reduced
real(r8),pointer :: minbufg(:) ! min buffer reduced
real(r8),pointer :: maxbufg(:) ! max buffer reduced
integer(i8),pointer :: isumbuf (:) ! integer local sum
integer(i8),pointer :: isumbufg(:) ! integer global sum
integer(i8) :: ihuge ! huge
integer(in) :: mpicom ! mpi comm
integer(in) :: iam ! pe number
integer(in) :: km,ka ! field indices
integer(in) :: ns ! size of local AV
integer(in) :: rcode ! error code
real(r8),pointer :: weight(:) ! weight
type(mct_string) :: mstring ! mct char type
character(CL) :: lcomment ! should be long enough
character(CL) :: itemc ! string converted to char
type(mct_avect) :: AV1 ! local avect with one field
type(mct_avect) :: AVr1 ! avect on root with one field
type(mct_avect) :: AVr2 ! avect on root with one field
!----- formats -----
character(*),parameter :: subName = '(seq_diag_avect_mct) '
character(*),parameter :: F00 = "('(seq_diag_avect_mct) ',4a)"
!-------------------------------------------------------------------------------
! print instantaneous budget data
!-------------------------------------------------------------------------------
call seq_cdata_setptrs(cdata,ID=ID,infodata=infodata,dom=dom,gsmap=gsmap)
call seq_comm_setptrs(ID,mpicom=mpicom,iam=iam)
call seq_infodata_GetData(infodata,bfbflag=bfbflag)
lcomment = ''
if (present(comment)) then
lcomment=trim(comment)
endif
ns = mct_aVect_lsize(AV)
npts = mct_aVect_lsize(dom%data)
if (ns /= npts) call shr_sys_abort(trim(subname)//' ERROR: size of AV,dom')
km = mct_aVect_indexRA(dom%data,'mask')
ka = mct_aVect_indexRA(dom%data,afldname)
kflds = mct_aVect_nRattr(AV)
allocate(sumbuf(kflds),sumbufg(kflds))
sumbuf = 0.0_r8
if (bfbflag) then
npts = mct_aVect_lsize(AV)
allocate(weight(npts))
weight(:) = 1.0_r8
do n = 1,npts
if (dom%data%rAttr(km,n) <= 1.0e-06_R8) then
weight(n) = 0.0_r8
else
weight(n) = dom%data%rAttr(ka,n)*shr_const_rearth*shr_const_rearth
endif
enddo
allocate(maxbuf(kflds),maxbufg(kflds))
maxbuf = 0.0_r8
do n = 1,npts
do k = 1,kflds
if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
maxbuf(k) = max(maxbuf(k),abs(AV%rAttr(k,n)*weight(n)))
endif
enddo
enddo
call shr_mpi_max(maxbuf,maxbufg,mpicom,subname,all=.true.)
call shr_mpi_sum(npts,nptsg,mpicom,subname,all=.true.)
do k = 1,kflds
if (maxbufg(k) < 1000.0*TINY(maxbufg(k)) .or. &
maxbufg(k) > HUGE(maxbufg(k))/(2.0_r8*nptsg)) then
maxbufg(k) = 0.0_r8
else
maxbufg(k) = (1.1_r8) * maxbufg(k) * nptsg
endif
enddo
allocate(isumbuf(kflds),isumbufg(kflds))
isumbuf = 0
ihuge = HUGE(isumbuf)
do n = 1,npts
do k = 1,kflds
if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
if (abs(maxbufg(k)) > 1000.0_r8 * TINY(maxbufg)) then
isumbuf(k) = isumbuf(k) + int((AV%rAttr(k,n)*weight(n)/maxbufg(k))*ihuge,i8)
endif
endif
enddo
enddo
call shr_mpi_sum(isumbuf,isumbufg,mpicom,subname)
do k = 1,kflds
sumbufg(k) = isumbufg(k)*maxbufg(k)/ihuge
enddo
deallocate(weight)
deallocate(maxbuf,maxbufg)
deallocate(isumbuf,isumbufg)
#if (1 == 0)
! OLD VERSION
call mct_aVect_init(AV1,rList='varf1',lsize=ns)
AV1%rAttr(1,1:ns) = dom%data%rAttr(km,1:ns) ! mask = AVr1
call mct_aVect_gather(AV1,AVr1,gsmap,0,mpicom,rcode)
AV1%rAttr(1,1:ns) = dom%data%rAttr(ka,1:ns) ! area = AVr2
call mct_aVect_gather(AV1,AVr2,gsmap,0,mpicom,rcode)
! --- compute weight on root pe
if (iam == 0) then
npts = mct_aVect_lsize(AVr1)
allocate(weight(npts))
weight(:) = 1.0_r8
do n = 1,npts
if (AVr1%rAttr(1,n) <= 1.0e-06_R8) then
weight(n) = 0.0_r8
else
weight(n) = AVr2%rAttr(1,n)*shr_const_rearth*shr_const_rearth
endif
enddo
! write(logunit,*) trim(subname),'tcx1 ',ns,npts,km,ka,kflds,minval(AVr2%rAttr),maxval(AVr2%rAttr)
! write(logunit,*) trim(subname),'tcx2 ',size(AVr1%rAttr),size(AVr2%rAttr),minval(AVr1%rAttr),maxval(AVr1%rAttr)
endif
! --- gather and compute stats one field at a time (to minimize memory)
do k = 1,kflds
AV1%rAttr(1,1:ns) = AV%rAttr(k,1:ns) ! fld k = AVr2
call mct_aVect_gather(AV1,AVr2,gsmap,0,mpicom,rcode)
if (iam == 0) then
npts = mct_aVect_lsize(AVr2)
! write(logunit,*) trim(subname),'tcx3 ',ns,npts,k,size(AVr2%rAttr),minval(AVr2%rAttr),maxval(AVr2%rAttr)
! write(logunit,*) trim(subname),'tcx4 ',size(weight),minval(weight),maxval(weight)
do n = 1,npts
if (AVr2%rAttr(1,n) > 1.01_r8*shr_const_spval .or. &
AVr2%rAttr(1,n) < 0.99_r8*shr_const_spval) then
sumbuf(k) = sumbuf(k) + AVr2%rAttr(1,n)*weight(n)
endif
enddo
endif
enddo
! --- local copy, relevant only on root pe
sumbufg = sumbuf
call mct_avect_clean(AV1)
if (iam == 0) then
call mct_aVect_clean(AVr1)
call mct_aVect_clean(AVr2)
endif
#endif
else
npts = mct_aVect_lsize(AV)
allocate(weight(npts))
weight(:) = 1.0_r8
do n = 1,npts
if (dom%data%rAttr(km,n) <= 1.0e-06_R8) then
weight(n) = 0.0_r8
else
weight(n) = dom%data%rAttr(ka,n)*shr_const_rearth*shr_const_rearth
endif
enddo
do n = 1,npts
do k = 1,kflds
if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
sumbuf(k) = sumbuf(k) + AV%rAttr(k,n)*weight(n)
endif
enddo
enddo
!--- global reduction ---
call shr_mpi_sum(sumbuf,sumbufg,mpicom,subname)
deallocate(weight)
endif
if (iam == 0) then
! write(logunit,*) 'sdAV: *** writing ',trim(lcomment),': k fld min/max/sum ***'
do k = 1,kflds
call mct_aVect_getRList(mstring,k,AV)
itemc = mct_string_toChar(mstring)
call mct_string_clean(mstring)
if (len_trim(lcomment) > 0) then
write(logunit,100) 'xxx','sorr',k,sumbufg(k),trim(lcomment),trim(itemc)
else
write(logunit,101) 'xxx','sorr',k,sumbufg(k),trim(itemc)
endif
enddo
call shr_sys_flush(logunit)
endif
deallocate(sumbuf,sumbufg)
100 format('comm_diag ',a3,1x,a4,1x,i3,es26.19,1x,a,1x,a)
101 format('comm_diag ',a3,1x,a4,1x,i3,es26.19,1x,a)
end subroutine seq_diag_avect_mct
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: seq_diag_avdiff_mct - print global budget diagnostics
!
! !DESCRIPTION:
! Print global diagnostics for AV/ID.
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
SUBROUTINE seq_diag_avdiff_mct(AV1,AV2,ID,comment),1
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(mct_aVect) , intent(in) :: AV1
type(mct_aVect) , intent(in) :: AV2
integer , intent(in) :: ID
character(len=*), intent(in), optional :: comment
!EOP
!--- local ---
integer(in) :: n,k,n1,k1,n2,k2 ! counters
integer(in) :: iam ! pe number
integer(in) :: cnt ! counter
real(r8) :: adiff,rdiff ! diff values
type(mct_string) :: mstring ! mct char type
character(len=64):: lcomment ! should be long enough
!----- formats -----
character(*),parameter :: subName = '(seq_diag_avdiff_mct) '
character(*),parameter :: F00 = "('(seq_diag_avdiff_mct) ',4a)"
!-------------------------------------------------------------------------------
! print instantaneous budget data
!-------------------------------------------------------------------------------
call seq_comm_setptrs(ID,iam=iam)
lcomment = ''
if (present(comment)) then
lcomment=trim(comment)
endif
n1 = mct_aVect_lsize(AV1)
k1 = mct_aVect_nRattr(AV1)
n2 = mct_aVect_lsize(AV2)
k2 = mct_aVect_nRattr(AV2)
if (n1 /= n2 .or. k1 /= k2) then
write(s_logunit,*) subname,trim(lcomment),' AV sizes different ',n1,n2,k1,k2
return
endif
do k = 1,k1
cnt = 0
adiff = 0.
rdiff = 0.
do n = 1,n1
if (AV1%rAttr(k,n) /= AV2%rAttr(k,n)) then
cnt = cnt + 1
adiff = max(adiff, abs(AV1%rAttr(k,n)-AV2%rAttr(k,n)))
rdiff = max(rdiff, abs(AV1%rAttr(k,n)-AV2%rAttr(k,n))/(abs(AV1%rAttr(k,n))+abs(AV2%rAttr(k,n))))
endif
enddo
if (cnt > 0) then
call mct_aVect_getRList(mstring,k,AV1)
write(s_logunit,*) subname,trim(lcomment),' AVs fld k diff ', &
iam,mct_string_toChar(mstring),cnt,adiff,rdiff, &
minval(AV1%rAttr(k,:)),minval(AV1%rAttr(k,:)), &
maxval(AV1%rAttr(k,:)),maxval(AV2%rAttr(k,:))
call mct_string_clean(mstring)
endif
enddo
end subroutine seq_diag_avdiff_mct
!===============================================================================
end module seq_diag_mct