CESM_Driver: seq_diag_mct.F90 Source File

00001 !===============================================================================
00002 ! SVN $Id: seq_diag_mct.F90 19586 2009-11-23 21:50:15Z kauff $
00003 ! 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 $
00004 !===============================================================================
00005 !BOP ===========================================================================
00006 !
00007 ! !MODULE: seq_diag_mod -- computes spatial \& time averages of fluxed quatities
00008 !
00009 ! !DESCRIPTION:
00010 !    The coupler is required to do certain diagnostics, those calculations are
00011 !    located in this module.
00012 !
00013 ! !REMARKS:
00014 !    Sign convention:
00015 !       positive value <=> the model is gaining water, heat, momentum, etc.
00016 !    Unit convention:
00017 !       heat flux     ~ W/m^2
00018 !       momentum flux ~ N/m^2
00019 !       water flux    ~ (kg/s)/m^2
00020 !       salt  flux    ~ (kg/s)/m^2
00021 !
00022 ! !REVISION HISTORY:
00023 !    2008-jul-10 - T. Craig    - updated budget implementation
00024 !    2007-may-07 - B. Kauffman - initial port to cpl7. 
00025 !    2002-nov-21 - R. Jacob    - initial port to cpl6. 
00026 !    199x-mmm-dd - B. Kauffman - original version in cpl4.
00027 !
00028 ! !INTERFACE: ------------------------------------------------------------------
00029 
00030 module seq_diag_mct
00031   
00032 ! !USES:
00033 
00034    use shr_kind_mod, only: r8 => shr_kind_r8, in=>shr_kind_in
00035    use shr_kind_mod, only: i8 => shr_kind_i8,  cl=>shr_kind_cl
00036    use shr_sys_mod       ! system calls
00037    use shr_mpi_mod       ! mpi wrappers
00038    use shr_const_mod     ! shared constants
00039    use mct_mod           ! mct wrappers
00040    use esmf_mod
00041 
00042    use seq_flds_indices  ! shared field indicies
00043    use seq_comm_mct  ! mpi comm groups & related
00044    use seq_cdata_mod
00045    use seq_timemgr_mod
00046 
00047 
00048    implicit none
00049    private
00050 
00051 ! !PUBLIC TYPES:
00052 
00053    ! none
00054 
00055 !PUBLIC MEMBER FUNCTIONS:
00056 
00057    public seq_diag_zero_mct
00058    public seq_diag_atm_mct
00059    public seq_diag_lnd_mct
00060    public seq_diag_rtm_mct
00061    public seq_diag_ocn_mct
00062    public seq_diag_ice_mct
00063    public seq_diag_accum_mct
00064    public seq_diag_sum0_mct
00065    public seq_diag_print_mct
00066    public seq_diag_avect_mct
00067    public seq_diag_avdiff_mct
00068 
00069 !EOP
00070 
00071    !----------------------------------------------------------------------------
00072    ! Local data
00073    !----------------------------------------------------------------------------
00074 
00075    !----- local constants -----
00076    real(r8),parameter :: HFLXtoWFLX =  ! water flux implied by latent heat of fusion
00077      - (shr_const_ocn_ref_sal-shr_const_ice_ref_sal) / &
00078        (shr_const_ocn_ref_sal*shr_const_latice)
00079 
00080    !--- C for component ---
00081    !--- "r" is recieve in the coupler, "s" is send from the coupler
00082 
00083    integer(in),parameter :: c_size = 18
00084 
00085    integer(in),parameter :: c_atm_as   = 1 ! model index: atm
00086    integer(in),parameter :: c_atm_ar   = 2 ! model index: atm
00087    integer(in),parameter :: c_inh_is   = 3 ! model index: ice, northern
00088    integer(in),parameter :: c_inh_ir   = 4 ! model index: ice, northern
00089    integer(in),parameter :: c_ish_is   = 5 ! model index: ice, southern
00090    integer(in),parameter :: c_ish_ir   = 6 ! model index: ice, southern
00091    integer(in),parameter :: c_lnd_ls   = 7 ! model index: lnd
00092    integer(in),parameter :: c_lnd_lr   = 8 ! model index: lnd
00093    integer(in),parameter :: c_ocn_os   = 9 ! model index: ocn
00094    integer(in),parameter :: c_ocn_or   =10 ! model index: ocn
00095    ! --- on atm grid ---
00096    integer(in),parameter :: c_inh_as   =11 ! model index: ice, northern
00097    integer(in),parameter :: c_inh_ar   =12 ! model index: ice, northern
00098    integer(in),parameter :: c_ish_as   =13 ! model index: ice, southern
00099    integer(in),parameter :: c_ish_ar   =14 ! model index: ice, southern
00100    integer(in),parameter :: c_lnd_as   =15 ! model index: lnd
00101    integer(in),parameter :: c_lnd_ar   =16 ! model index: lnd
00102    integer(in),parameter :: c_ocn_as   =17 ! model index: ocn
00103    integer(in),parameter :: c_ocn_ar   =18 ! model index: ocn
00104 
00105    character(len=8),parameter :: cname(c_size) = 
00106       (/' c2a_atm',' a2c_atm',' c2i_inh',' i2c_inh',' c2i_ish',' i2c_ish', 
00107         ' c2l_lnd',' l2c_lnd',' c2o_ocn',' o2c_ocn', 
00108         ' c2a_inh',' a2c_inh',' c2a_ish',' a2c_ish', 
00109         ' c2a_lnd',' a2c_lnd',' c2a_ocn',' a2c_ocn' /)
00110 
00111    !--- F for field ---
00112 
00113    integer(in),parameter :: f_size = 17
00114    integer(in),parameter :: f_a    = 1    ! index for area
00115    integer(in),parameter :: f_h    = 2    ! 1st index for heat
00116    integer(in),parameter :: f_w    = 11   ! 1st index for water
00117 
00118    integer(in),parameter :: f_area    = 1 ! area (wrt to unit sphere)
00119    integer(in),parameter :: f_hfrz    = 2 ! heat : latent, freezing
00120    integer(in),parameter :: f_hmelt   = 3 ! heat : latent, melting
00121    integer(in),parameter :: f_hswnet  = 4 ! heat : short wave, net
00122    integer(in),parameter :: f_hlwdn   = 5 ! heat : longwave down
00123    integer(in),parameter :: f_hlwup   = 6 ! heat : longwave up
00124    integer(in),parameter :: f_hlatv   = 7 ! heat : latent, vaporization
00125    integer(in),parameter :: f_hlatf   = 8 ! heat : latent, fusion, snow       
00126    integer(in),parameter :: f_hioff   = 9 ! heat : latent, fusion, frozen runoff
00127    integer(in),parameter :: f_hsen    =10 ! heat : sensible
00128    integer(in),parameter :: f_wfrz    =11 ! water: freezing
00129    integer(in),parameter :: f_wmelt   =12 ! water: melting
00130    integer(in),parameter :: f_wrain   =13 ! water: precip, liquid
00131    integer(in),parameter :: f_wsnow   =14 ! water: precip, frozen
00132    integer(in),parameter :: f_wevap   =15 ! water: evaporation
00133    integer(in),parameter :: f_wroff   =16 ! water: runoff
00134    integer(in),parameter :: f_wioff   =17 ! water: frozen runoff
00135 
00136    character(len=8),parameter :: fname(f_size) = 
00137       (/'    area',' hfreeze','   hmelt','  hnetsw','   hlwdn', 
00138         '   hlwup',' hlatvap',' hlatfus','  hiroff','    hsen', 
00139         ' wfreeze','   wmelt','   wrain','   wsnow', 
00140         '   wevap',' wrunoff',' wfrzrof' /)
00141 
00142    !--- P for period ---
00143 
00144    integer(in),parameter :: p_size = 5
00145 
00146    integer(in),parameter :: p_inst = 1
00147    integer(in),parameter :: p_day  = 2
00148    integer(in),parameter :: p_mon  = 3
00149    integer(in),parameter :: p_ann  = 4
00150    integer(in),parameter :: p_inf  = 5
00151 
00152    character(len=8),parameter :: pname(p_size) = 
00153       (/'    inst','   daily',' monthly','  annual','all_time' /)
00154 
00155 ! !PUBLIC DATA MEMBERS
00156 
00157    !--- time-averaged (annual?) global budge diagnostics ---
00158    !--- note: call sum0 then save budg_dataG and budg_ns on restart from/to root pe ---
00159    real(r8),public :: budg_dataL(f_size,c_size,p_size) ! local sum, valid on all pes
00160    real(r8),public :: budg_dataG(f_size,c_size,p_size) ! global sum, valid only on root pe
00161    real(r8),public :: budg_ns   (f_size,c_size,p_size) ! counter, valid only on root pe
00162 
00163    character(len=*),parameter :: afldname  = 'aream'
00164    character(len=*),parameter :: latname   = 'lat'
00165    character(len=*),parameter :: afracname = 'afrac'
00166    character(len=*),parameter :: lfracname = 'lfrac'
00167    character(len=*),parameter :: ofracname = 'ofrac'
00168    character(len=*),parameter :: ifracname = 'ifrac'
00169    character(len=*),parameter :: ascalname = 'ascale'
00170    character(len=*),parameter :: lfrinname = 'lfrin'
00171 
00172    character(*),parameter :: modName = "(seq_diag_mct) "
00173 
00174    integer(in),parameter :: debug = 0 ! internal debug level
00175 
00176 !===============================================================================
00177 contains
00178 !===============================================================================
00179 
00180 !===============================================================================
00181 !BOP ===========================================================================
00182 !
00183 ! !IROUTINE: seq_diag_zero_mct - zero out global budget diagnostic data.
00184 !
00185 ! !DESCRIPTION:
00186 !    Zero out global budget diagnostic data.
00187 !
00188 ! !REVISION HISTORY:
00189 !    2008-jul-11 - T. Craig - update
00190 !
00191 ! !INTERFACE: ------------------------------------------------------------------
00192 
00193 subroutine seq_diag_zero_mct(EClock,mode)
00194 
00195 ! !INPUT/OUTPUT PARAMETERS:
00196 
00197    type(ESMF_Clock), intent(in),optional :: EClock
00198    character(len=*), intent(in),optional :: mode
00199 
00200 !EOP
00201 
00202    integer(IN) :: ip,yr,mon,day,sec
00203    !----- formats -----
00204    character(*),parameter :: subName = '(seq_diag_zero_mct) '
00205 
00206 !-------------------------------------------------------------------------------
00207 !
00208 !-------------------------------------------------------------------------------
00209 
00210    if (.not. present(EClock) .and. .not. present(mode)) then
00211       call shr_sys_abort(subName//' ERROR EClock or mode should be present')
00212    endif
00213 
00214    if (present(EClock)) then
00215       call seq_timemgr_EClockGetData(EClock,curr_yr=yr, &
00216          curr_mon=mon,curr_day=day,curr_tod=sec)
00217 
00218       do ip = 1,p_size
00219          if (ip == p_inst) then
00220             budg_dataL(:,:,ip) = 0.0_r8
00221             budg_dataG(:,:,ip) = 0.0_r8
00222             budg_ns(:,:,ip) = 0.0_r8
00223          endif
00224          if (ip==p_day .and. sec==0) then
00225             budg_dataL(:,:,ip) = 0.0_r8
00226             budg_dataG(:,:,ip) = 0.0_r8
00227             budg_ns(:,:,ip) = 0.0_r8
00228          endif
00229          if (ip==p_mon .and. day==1 .and. sec==0) then
00230             budg_dataL(:,:,ip) = 0.0_r8
00231             budg_dataG(:,:,ip) = 0.0_r8
00232             budg_ns(:,:,ip) = 0.0_r8
00233          endif
00234          if (ip==p_ann .and. mon==1 .and. day==1 .and. sec==0) then
00235             budg_dataL(:,:,ip) = 0.0_r8
00236             budg_dataG(:,:,ip) = 0.0_r8
00237             budg_ns(:,:,ip) = 0.0_r8
00238          endif
00239       enddo
00240    endif
00241 
00242    if (present(mode)) then
00243       if (trim(mode) == 'inst') then
00244          budg_dataL(:,:,p_inst) = 0.0_r8
00245          budg_dataG(:,:,p_inst) = 0.0_r8
00246          budg_ns(:,:,p_inst) = 0.0_r8
00247       elseif (trim(mode) == 'day') then
00248          budg_dataL(:,:,p_day) = 0.0_r8
00249          budg_dataG(:,:,p_day) = 0.0_r8
00250          budg_ns(:,:,p_day) = 0.0_r8
00251       elseif (trim(mode) == 'mon') then
00252          budg_dataL(:,:,p_mon) = 0.0_r8
00253          budg_dataG(:,:,p_mon) = 0.0_r8
00254          budg_ns(:,:,p_mon) = 0.0_r8
00255       elseif (trim(mode) == 'ann') then
00256          budg_dataL(:,:,p_ann) = 0.0_r8
00257          budg_dataG(:,:,p_ann) = 0.0_r8
00258          budg_ns(:,:,p_ann) = 0.0_r8
00259       elseif (trim(mode) == 'inf') then
00260          budg_dataL(:,:,p_inf) = 0.0_r8
00261          budg_dataG(:,:,p_inf) = 0.0_r8
00262          budg_ns(:,:,p_inf) = 0.0_r8
00263       elseif (trim(mode) == 'all') then
00264          budg_dataL(:,:,:) = 0.0_r8
00265          budg_dataG(:,:,:) = 0.0_r8
00266          budg_ns(:,:,:) = 0.0_r8
00267       else
00268          call shr_sys_abort(subname//' ERROR in mode '//trim(mode))
00269       endif
00270    endif
00271 
00272 end subroutine seq_diag_zero_mct
00273 
00274 !===============================================================================
00275 !BOP ===========================================================================
00276 !
00277 ! !IROUTINE: seq_diag_accum_mct - accum out global budget diagnostic data.
00278 !
00279 ! !DESCRIPTION:
00280 !    Accum out global budget diagnostic data.
00281 !
00282 ! !REVISION HISTORY:
00283 !    2008-jul-11 - T. Craig - update
00284 !
00285 ! !INTERFACE: ------------------------------------------------------------------
00286 
00287 subroutine seq_diag_accum_mct()
00288 
00289 ! !INPUT/OUTPUT PARAMETERS:
00290 
00291 !EOP
00292 
00293    integer(in) :: ip
00294 
00295    !----- formats -----
00296    character(*),parameter :: subName = '(seq_diag_accum_mct) '
00297 
00298 !-------------------------------------------------------------------------------
00299 !
00300 !-------------------------------------------------------------------------------
00301 
00302    do ip = p_inst+1,p_size
00303       budg_dataL(:,:,ip) = budg_dataL(:,:,ip) + budg_dataL(:,:,p_inst)
00304    enddo
00305    budg_ns(:,:,:) = budg_ns(:,:,:) + 1.0_r8
00306 
00307 end subroutine seq_diag_accum_mct
00308 
00309 !===============================================================================
00310 !BOP ===========================================================================
00311 !
00312 ! !IROUTINE: seq_diag_sum0_mct - sum local to global on root
00313 !
00314 ! !DESCRIPTION:
00315 !    Sum local values to global on root
00316 !
00317 ! !REVISION HISTORY:
00318 !    2008-jul-19 - T. Craig - update
00319 !
00320 ! !INTERFACE: ------------------------------------------------------------------
00321 
00322 subroutine seq_diag_sum0_mct()
00323 
00324 ! !INPUT/OUTPUT PARAMETERS:
00325 
00326 !EOP
00327 
00328    real(r8) :: budg_dataGtmp(f_size,c_size,p_size) ! temporary sum
00329    integer(in)      :: mpicom      ! mpi comm
00330    !----- formats -----
00331    character(*),parameter :: subName = '(seq_diag_sum0_mct) '
00332 
00333 !-------------------------------------------------------------------------------
00334 !
00335 !-------------------------------------------------------------------------------
00336 
00337    call seq_comm_setptrs(CPLID,mpicom=mpicom)
00338    budg_dataGtmp = 0.0_r8
00339    call shr_mpi_sum(budg_dataL,budg_dataGtmp,mpicom,subName)
00340    budg_dataG = budg_dataG + budg_dataGtmp
00341    budg_dataL = 0.0_r8
00342 
00343 end subroutine seq_diag_sum0_mct
00344 
00345 !===============================================================================
00346 !BOP ===========================================================================
00347 !
00348 ! !IROUTINE: seq_diag_atm_mct - compute global atm input/output flux diagnostics
00349 !
00350 ! !DESCRIPTION:
00351 !     Compute global atm input/output flux diagnostics
00352 !
00353 ! !REVISION HISTORY:
00354 !    2008-jul-10 - T. Craig - update
00355 !
00356 ! !INTERFACE: ------------------------------------------------------------------
00357 
00358 subroutine seq_diag_atm_mct( dom_a, frac_a, a2x_a, x2a_a )
00359 
00360 ! !INPUT/OUTPUT PARAMETERS:
00361 
00362    type(mct_gGrid),intent(in)          ::  dom_a ! model domain
00363    type(mct_aVect),intent(in)          :: frac_a ! domain fractions
00364    type(mct_aVect),intent(in),optional ::  a2x_a ! model to drv bundle
00365    type(mct_aVect),intent(in),optional ::  x2a_a ! drv to model bundle
00366 
00367 !EOP
00368 
00369    !----- local -----
00370    integer(in)      :: k,n,ic,if,ip      ! generic index
00371    integer(in)      :: kArea             ! index of area field in aVect
00372    integer(in)      :: kLat              ! index of lat field in aVect
00373    integer(in)      :: kl,ka,ko,ki       ! fraction indices
00374    integer(in)      :: lSize             ! size of aVect
00375    real(r8)         :: da,di,do,dl       ! area of a grid cell
00376 
00377    !----- formats -----
00378    character(*),parameter :: subName = '(seq_diag_atm_mct) '
00379 
00380 !-------------------------------------------------------------------------------
00381 !
00382 !-------------------------------------------------------------------------------
00383 
00384    if (.not. present(a2x_a) .and. .not. present(x2a_a)) then
00385       call shr_sys_abort(subName//"ERROR: must input a bundle")
00386    end if
00387 
00388    kArea = mct_aVect_indexRA(dom_a%data,afldname)
00389    kLat  = mct_aVect_indexRA(dom_a%data,latname)
00390    ka    = mct_aVect_indexRA(frac_a,afracname)
00391    kl    = mct_aVect_indexRA(frac_a,lfracname)
00392    ko    = mct_aVect_indexRA(frac_a,ofracname)
00393    ki    = mct_aVect_indexRA(frac_a,ifracname)
00394 
00395    !---------------------------------------------------------------------------
00396    ! add values found in this bundle to the budget table
00397    !---------------------------------------------------------------------------
00398 
00399    ip = p_inst
00400 
00401    if (present(a2x_a)) then
00402       lSize = mct_avect_lSize(a2x_a)
00403       do n=1,lSize
00404       do k=1,4
00405 
00406          if (k == 1) then
00407             ic = c_atm_ar
00408             da = -dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ka,n)
00409          elseif (k == 2) then
00410             ic = c_lnd_ar
00411             da =  dom_a%data%rAttr(kArea,n) * frac_a%rAttr(kl,n)
00412          elseif (k == 3) then
00413             ic = c_ocn_ar
00414             da =  dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ko,n)
00415          elseif (k == 4) then
00416             if (dom_a%data%rAttr(kLat,n) > 0.0_r8) then
00417                ic = c_inh_ar
00418             else
00419                ic = c_ish_ar
00420             endif
00421             da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ki,n)
00422          endif
00423 
00424          if = f_area  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da
00425          if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_swnet,n)
00426          if = f_hlwdn ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_lwdn,n)
00427          if = f_wrain ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_rainc,n) &
00428                                                                     + da*a2x_a%rAttr(index_a2x_Faxa_rainl,n)
00429          if = f_wsnow ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*a2x_a%rAttr(index_a2x_Faxa_snowc,n) &
00430                                                                     + da*a2x_a%rAttr(index_a2x_Faxa_snowl,n)
00431       enddo
00432       enddo
00433       ! --- heat implied by snow flux ---
00434       ic = c_atm_ar;  budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00435       ic = c_lnd_ar;  budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00436       ic = c_ocn_ar;  budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00437       ic = c_inh_ar;  budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00438       ic = c_ish_ar;  budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00439    end if
00440 
00441    if (present(x2a_a)) then
00442       lSize = mct_avect_lSize(x2a_a)
00443       do n=1,lSize
00444       do k=1,4
00445 
00446          if (k == 1) then
00447             ic = c_atm_as
00448             da = -dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ka,n)
00449          elseif (k == 2) then
00450             ic = c_lnd_as
00451             da =  dom_a%data%rAttr(kArea,n) * frac_a%rAttr(kl,n)
00452          elseif (k == 3) then
00453             ic = c_ocn_as
00454             da =  dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ko,n)
00455          elseif (k == 4) then
00456             if (dom_a%data%rAttr(kLat,n) > 0.0_r8) then
00457                ic = c_inh_as
00458             else
00459                ic = c_ish_as
00460             endif
00461             da = dom_a%data%rAttr(kArea,n) * frac_a%rAttr(ki,n)
00462          endif
00463 
00464          if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da
00465          if = f_hlwup; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_lwup,n)
00466          if = f_hlatv; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_lat,n)
00467          if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_sen,n)
00468          if = f_wevap; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + da*x2a_a%rAttr(index_x2a_Faxx_evap,n)
00469 
00470       enddo
00471       enddo
00472    end if
00473 
00474 end subroutine seq_diag_atm_mct
00475 
00476 !===============================================================================
00477 !BOP ===========================================================================
00478 !
00479 ! !IROUTINE: seq_diag_lnd_mct - compute global lnd input/output flux diagnostics
00480 !
00481 ! !DESCRIPTION:
00482 !     Compute global lnd input/output flux diagnostics
00483 !
00484 ! !REVISION HISTORY:
00485 !    2008-jul-10 - T. Craig - update
00486 !
00487 ! !INTERFACE: ------------------------------------------------------------------
00488 
00489 subroutine seq_diag_lnd_mct( dom_l, frac_l, l2x_l, x2l_l)
00490 
00491    type(mct_gGrid),intent(in)          ::  dom_l ! model domain
00492    type(mct_aVect),intent(in)          :: frac_l ! frac bundle
00493    type(mct_aVect),intent(in),optional ::  l2x_l ! model to drv bundle
00494    type(mct_aVect),intent(in),optional ::  x2l_l ! drv to model bundle
00495 
00496 !EOP
00497 
00498    !----- local -----
00499    integer(in)      :: k,n,ic,if,ip      ! generic index
00500    integer(in)      :: kArea             ! index of area field in aVect
00501    integer(in)      :: kLat              ! index of lat field in aVect
00502    integer(in)      :: kl,ka,ko,ki       ! fraction indices
00503    integer(in)      :: lSize             ! size of aVect
00504    real(r8)         :: da,di,do,dl       ! area of a grid cell
00505 
00506    !----- formats -----
00507    character(*),parameter :: subName = '(seq_diag_lnd_mct) '
00508 
00509 !-------------------------------------------------------------------------------
00510 !
00511 !-------------------------------------------------------------------------------
00512 
00513    if (.not. present(l2x_l) .and. .not. present(x2l_l)) then
00514       call shr_sys_abort(subName//"ERROR: must input a bundle")
00515    end if
00516 
00517    !---------------------------------------------------------------------------
00518    ! add values found in this bundle to the budget table
00519    !---------------------------------------------------------------------------
00520 
00521    ip = p_inst
00522 
00523    kArea = mct_aVect_indexRA(dom_l%data,afldname)
00524    kl    = mct_aVect_indexRA(frac_l,lfrinname)
00525    ka    = mct_aVect_indexRA(dom_l%data,ascalname)
00526 
00527    if (present(l2x_l)) then
00528       lSize = mct_avect_lSize(l2x_l)
00529       ic = c_lnd_lr
00530       do n=1,lSize
00531          dl =  dom_l%data%rAttr(kArea,n) * frac_l%rAttr(kl,n) * dom_l%data%rAttr(ka,n)
00532          if = f_area  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl
00533          if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_swnet,n)
00534          if = f_hlwup ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_lwup,n)
00535          if = f_hlatv ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_lat,n)
00536          if = f_hsen  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_sen,n)
00537          if = f_wevap ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*l2x_l%rAttr(index_l2x_Fall_evap,n)
00538       end do
00539    end if
00540 
00541    if (present(x2l_l)) then
00542       lSize = mct_avect_lSize(x2l_l)
00543       ic = c_lnd_ls
00544       do n=1,lSize
00545          dl =  dom_l%data%rAttr(kArea,n) * frac_l%rAttr(kl,n) * dom_l%data%rAttr(ka,n)
00546          if = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl
00547          if = f_hlwdn; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_lwdn,n)
00548          if = f_wrain; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_rainc,n) &
00549                                                                    + dl*x2l_l%rAttr(index_x2l_Faxa_rainl,n)
00550          if = f_wsnow; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + dl*x2l_l%rAttr(index_x2l_Faxa_snowc,n) &
00551                                                                    + dl*x2l_l%rAttr(index_x2l_Faxa_snowl,n)
00552       end do
00553       budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00554    end if
00555 
00556 end subroutine seq_diag_lnd_mct
00557 
00558 !===============================================================================
00559 !BOP ===========================================================================
00560 !
00561 ! !IROUTINE: seq_diag_rtm_mct - compute global rtm input/output flux diagnostics
00562 !
00563 ! !DESCRIPTION:
00564 !     Compute global rtm input/output flux diagnostics
00565 !
00566 ! !REVISION HISTORY:
00567 !    2008-jul-10 - T. Craig - update
00568 !
00569 ! !INTERFACE: ------------------------------------------------------------------
00570 
00571 subroutine seq_diag_rtm_mct( dom_r, r2x_r)
00572 
00573    type(mct_gGrid),intent(in)          ::  dom_r ! model domain
00574    type(mct_aVect),intent(in)          ::  r2x_r ! model to drv bundle
00575 
00576 !EOP
00577 
00578    !----- local -----
00579    integer(in)      :: k,n,ic,if,ip      ! generic index
00580    integer(in)      :: kArea             ! index of area field in aVect
00581    integer(in)      :: kLat              ! index of lat field in aVect
00582    integer(in)      :: kl,ka,ko,ki       ! fraction indices
00583    integer(in)      :: lSize             ! size of aVect
00584    real(r8)         :: da,di,do,dl       ! area of a grid cell
00585 
00586    !----- formats -----
00587    character(*),parameter :: subName = '(seq_diag_rtm_mct) '
00588 
00589 !-------------------------------------------------------------------------------
00590 !
00591 !-------------------------------------------------------------------------------
00592 
00593    !---------------------------------------------------------------------------
00594    ! add values found in this bundle to the budget table
00595    !---------------------------------------------------------------------------
00596 
00597    ip = p_inst
00598    ic = c_lnd_lr
00599    kArea = mct_aVect_indexRA(dom_r%data,afldname)
00600    lSize = mct_avect_lSize(r2x_r)
00601    do n=1,lSize
00602       dl =  dom_r%data%rAttr(kArea,n)
00603       if = f_wroff; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - dl*r2x_r%rAttr(index_r2x_Forr_roff,n)
00604       if = f_wioff; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - dl*r2x_r%rAttr(index_r2x_Forr_ioff,n)
00605    end do
00606    budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice
00607 
00608 end subroutine seq_diag_rtm_mct
00609 
00610 !BOP ===========================================================================
00611 !
00612 ! !IROUTINE: seq_diag_ocn_mct - compute global ocn input/output flux diagnostics
00613 !
00614 ! !DESCRIPTION:
00615 !     Compute global ocn input/output flux diagnostics
00616 !
00617 ! !REVISION HISTORY:
00618 !    2008-jul-10 - T. Craig - update
00619 !
00620 ! !INTERFACE: ------------------------------------------------------------------
00621 
00622 subroutine seq_diag_ocn_mct( dom_o, frac_o, o2x_o, x2o_o, xao_o, r2x_o)
00623 
00624    type(mct_gGrid),intent(in)          ::  dom_o ! model domain
00625    type(mct_aVect),intent(in)          :: frac_o ! frac bundle
00626    type(mct_aVect),intent(in),optional ::  o2x_o ! model to drv bundle
00627    type(mct_aVect),intent(in),optional ::  x2o_o ! drv to model bundle
00628    type(mct_aVect),intent(in),optional ::  xao_o ! drv to model bundle
00629    type(mct_aVect),intent(in),optional ::  r2x_o ! roff to drv bundle
00630 
00631 !EOP
00632 
00633    !----- local -----
00634    integer(in)      :: k,n,if,ic,ip      ! generic index
00635    integer(in)      :: kArea             ! index of area field in aVect
00636    integer(in)      :: kLat              ! index of lat field in aVect
00637    integer(in)      :: kl,ka,ko,ki       ! fraction indices
00638    integer(in)      :: lSize             ! size of aVect
00639    real(r8)         :: da,di,do,dl       ! area of a grid cell
00640 
00641    !----- formats -----
00642    character(*),parameter :: subName = '(seq_diag_ocn_mct) '
00643 
00644 !-------------------------------------------------------------------------------
00645 !
00646 !-------------------------------------------------------------------------------
00647 
00648    if (.not. present(o2x_o) .and. .not. present(x2o_o) .and. .not. present(xao_o)) then
00649       call shr_sys_abort(subName//"ERROR: must input a bundle")
00650    end if
00651 
00652    !---------------------------------------------------------------------------
00653    ! add values found in this bundle to the budget table
00654    !---------------------------------------------------------------------------
00655 
00656    ip = p_inst
00657 
00658    kArea = mct_aVect_indexRA(dom_o%data,afldname)
00659    ko    = mct_aVect_indexRA(frac_o,ofracname)
00660    ki    = mct_aVect_indexRA(frac_o,ifracname)
00661 
00662    if (present(o2x_o)) then
00663       lSize = mct_avect_lSize(o2x_o)
00664       ic = c_ocn_or
00665       do n=1,lSize
00666          do =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
00667          di =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
00668          if = f_area; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do
00669          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))
00670       end do
00671       budg_dataL(f_wfrz,ic,ip) = budg_dataL(f_hfrz,ic,ip) * HFLXtoWFLX
00672    end if
00673 
00674    if (present(xao_o)) then
00675       lSize = mct_avect_lSize(xao_o)
00676       ic = c_ocn_or
00677       do n=1,lSize
00678          do =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
00679          if = f_hlwup; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_lwup,n)
00680          if = f_hlatv; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_lat,n)
00681          if = f_hsen ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_sen,n)
00682          if = f_wevap; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do*xao_o%rAttr(index_xao_Faox_evap,n)
00683       end do
00684    end if
00685 
00686    if (present(x2o_o)) then
00687       lSize = mct_avect_lSize(x2o_o)
00688       ic = c_ocn_os
00689       do n=1,lSize
00690          do =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
00691          di =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
00692          if = f_area  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + do
00693          if = f_hmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_melth,n)
00694          if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_swnet,n)
00695          if = f_hlwdn ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_lwdn,n)
00696          if = f_wmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_meltw,n)
00697          if = f_wrain ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_rain,n)
00698          if = f_wsnow ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Foxx_snow,n)
00699       !  if = f_wroff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Forr_roff,n)
00700       !  if = f_wioff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*x2o_o%rAttr(index_x2o_Forr_ioff,n)
00701       end do
00702       budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00703    end if
00704 
00705    if (present(r2x_o)) then
00706       lSize = mct_avect_lSize(r2x_o)
00707       ic = c_ocn_os
00708       do n=1,lSize
00709          do =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
00710          di =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
00711          if = f_wroff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*r2x_o%rAttr(index_r2x_Forr_roff,n)
00712          if = f_wioff ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + (do+di)*r2x_o%rAttr(index_r2x_Forr_ioff,n)
00713       end do
00714       budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice
00715    end if
00716 
00717 end subroutine seq_diag_ocn_mct
00718 
00719 !===============================================================================
00720 !BOP ===========================================================================
00721 !
00722 ! !IROUTINE: seq_diag_ice_mct - compute global ice input/output flux diagnostics
00723 !
00724 ! !DESCRIPTION:
00725 !     Compute global ice input/output flux diagnostics
00726 !
00727 ! !REVISION HISTORY:
00728 !    2008-jul-10 - T. Craig - update
00729 !
00730 ! !INTERFACE: ------------------------------------------------------------------
00731 
00732 subroutine seq_diag_ice_mct( dom_i, frac_i, i2x_i, x2i_i)
00733 
00734    type(mct_gGrid),intent(in)          ::  dom_i ! model domain
00735    type(mct_aVect),intent(in)          :: frac_i ! frac bundle
00736    type(mct_aVect),intent(in),optional ::  i2x_i ! model to drv bundle
00737    type(mct_aVect),intent(in),optional ::  x2i_i ! drv to model bundle
00738 
00739 !EOP
00740 
00741    !----- local -----
00742    integer(in)      :: k,n,ic,if,ip      ! generic index
00743    integer(in)      :: kArea             ! index of area field in aVect
00744    integer(in)      :: kLat              ! index of lat field in aVect
00745    integer(in)      :: kl,ka,ko,ki       ! fraction indices
00746    integer(in)      :: lSize             ! size of aVect
00747    real(r8)         :: da,di,do,dl       ! area of a grid cell
00748 
00749    !----- formats -----
00750    character(*),parameter :: subName = '(seq_diag_ice_mct) '
00751 
00752 !-------------------------------------------------------------------------------
00753 !
00754 !-------------------------------------------------------------------------------
00755 
00756    if (.not. present(i2x_i) .and. .not. present(x2i_i)) then
00757       call shr_sys_abort(subName//"ERROR: must input a bundle")
00758    end if
00759 
00760    !---------------------------------------------------------------------------
00761    ! add values found in this bundle to the budget table
00762    !---------------------------------------------------------------------------
00763 
00764    ip = p_inst
00765 
00766    kArea = mct_aVect_indexRA(dom_i%data,afldname)
00767    kLat  = mct_aVect_indexRA(dom_i%data,latname)
00768    ki    = mct_aVect_indexRA(frac_i,ifracname)
00769    ko    = mct_aVect_indexRA(frac_i,ofracname)
00770 
00771    if (present(i2x_i)) then
00772       lSize = mct_avect_lSize(i2x_i)
00773       do n=1,lSize
00774          if (dom_i%data%rAttr(kLat,n) > 0.0_r8) then
00775             ic = c_inh_ir
00776          else
00777             ic = c_ish_ir
00778          endif
00779          di =  dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ki,n)
00780          if = f_area  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di
00781          if = f_hmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - di*i2x_i%rAttr(index_i2x_Fioi_melth,n)
00782          if = f_wmelt ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) - di*i2x_i%rAttr(index_i2x_Fioi_meltw,n)
00783          if = f_hswnet; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_swnet,n) &
00784                                                                     - di*i2x_i%rAttr(index_i2x_Fioi_swpen,n)
00785          if = f_hlwup ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_lwup,n)
00786          if = f_hlatv ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_lat,n)
00787          if = f_hsen  ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_sen,n)
00788          if = f_wevap ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*i2x_i%rAttr(index_i2x_Faii_evap,n)
00789       end do
00790    end if
00791 
00792    if (present(x2i_i)) then
00793       lSize = mct_avect_lSize(x2i_i)
00794       do n=1,lSize
00795          if (dom_i%data%rAttr(kLat,n) > 0.0_r8) then
00796             ic = c_inh_is
00797          else
00798             ic = c_ish_is
00799          endif
00800          do =  dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ko,n)
00801          di =  dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ki,n)
00802          if  = f_area ; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di
00803          if  = f_hlwdn; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_lwdn,n)
00804          if  = f_wrain; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_rain,n)
00805          if  = f_wsnow; budg_dataL(if,ic,ip) = budg_dataL(if,ic,ip) + di*x2i_i%rAttr(index_x2i_Faxa_snow,n)
00806          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))
00807       end do
00808       ic = c_inh_is  
00809       budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00810       budg_dataL(f_wfrz ,ic,ip) =  budg_dataL(f_hfrz ,ic,ip)*HFLXtoWFLX
00811       ic = c_ish_is
00812       budg_dataL(f_hlatf,ic,ip) = -budg_dataL(f_wsnow,ic,ip)*shr_const_latice
00813       budg_dataL(f_wfrz ,ic,ip) =  budg_dataL(f_hfrz ,ic,ip)*HFLXtoWFLX
00814    end if
00815 
00816 end subroutine seq_diag_ice_mct
00817 
00818 !===============================================================================
00819 !BOP ===========================================================================
00820 !
00821 ! !IROUTINE: seq_diag_print_mct - print global budget diagnostics
00822 !
00823 ! !DESCRIPTION:
00824 !   Print global budget diagnostics.
00825 !
00826 ! !REVISION HISTORY:
00827 !
00828 ! !INTERFACE: ------------------------------------------------------------------
00829 
00830 SUBROUTINE seq_diag_print_mct(EClock,stop_alarm, &
00831    budg_print_inst, budg_print_daily, budg_print_month, &
00832    budg_print_ann, budg_print_ltann, budg_print_ltend)
00833 
00834    implicit none
00835 
00836 ! !INPUT/OUTPUT PARAMETERS:
00837 
00838    type(ESMF_Clock), intent(in) :: EClock
00839    logical,intent(in)           :: stop_alarm
00840    integer,intent(in)           :: budg_print_inst
00841    integer,intent(in)           :: budg_print_daily
00842    integer,intent(in)           :: budg_print_month
00843    integer,intent(in)           :: budg_print_ann
00844    integer,intent(in)           :: budg_print_ltann
00845    integer,intent(in)           :: budg_print_ltend
00846 
00847 !EOP
00848 
00849    !--- local ---
00850    integer(in)      :: ic,if,ip    ! data array indicies
00851    integer(in)      :: ica,icl,icn,ics,ico
00852    integer(in)      :: icar,icxs,icxr,icas
00853    integer(in)      :: n           ! loop counter
00854    integer(in)      :: nday        ! number of days in time avg
00855    integer(in)      :: cdate,sec   ! coded date, seconds
00856    integer(in)      :: yr,mon,day  ! date
00857    integer(in)      :: iam         ! pe number
00858    integer(in)      :: plev        ! print level
00859    logical          :: sumdone     ! has a sum been computed yet
00860    character(len=40):: str         ! string
00861    real(r8) :: dataGpr (f_size,c_size,p_size) ! values to print, scaled and such
00862 
00863    !----- formats -----
00864    character(*),parameter :: subName = '(seq_diag_print_mct) '
00865    character(*),parameter :: F00   = "('(seq_diag_print_mct) ',4a)"
00866 
00867    !----- formats -----
00868    character(*),parameter :: FAH="(4a,i9,i6)"
00869    character(*),parameter :: FA0="('    ',8x,6(6x,a8,1x))"
00870    character(*),parameter :: FA1="('    ',a8,6f15.8)"
00871 
00872 !-------------------------------------------------------------------------------
00873 ! print instantaneous budget data
00874 !-------------------------------------------------------------------------------
00875 
00876    sumdone = .false.
00877    call seq_comm_setptrs(CPLID,iam=iam)
00878    call seq_timemgr_EClockGetData(EClock,curr_yr=yr, &
00879       curr_mon=mon,curr_day=day,curr_tod=sec)
00880    cdate = yr*10000+mon*100+day
00881 
00882    do ip = 1,p_size
00883       plev = 0
00884       if (ip == p_inst) then
00885          plev = max(plev,budg_print_inst)
00886       endif
00887       if (ip==p_day .and. sec==0) then
00888          plev = max(plev,budg_print_daily)
00889       endif
00890       if (ip==p_mon .and. day==1 .and. sec==0) then
00891          plev = max(plev,budg_print_month)
00892       endif
00893       if (ip==p_ann .and. mon==1 .and. day==1 .and. sec==0) then
00894          plev = max(plev,budg_print_ann)
00895       endif
00896       if (ip==p_inf .and. mon==1 .and. day==1 .and. sec==0) then
00897          plev = max(plev,budg_print_ltann)
00898       endif
00899       if (ip==p_inf .and. stop_alarm) then
00900          plev = max(plev,budg_print_ltend)
00901       endif
00902 
00903    if (plev > 0) then
00904 ! ---- doprint ---- doprint ---- doprint ----
00905 
00906    if (.not.sumdone) then
00907       call seq_diag_sum0_mct()
00908       dataGpr = budg_dataG
00909       sumdone = .true.
00910 
00911    !  old budget normalizations (global area and 1e6 for water)
00912       dataGpr = dataGpr/(4.0_r8*shr_const_pi)
00913       dataGpr(f_w:f_size,:,:) = dataGpr(f_w:f_size,:,:) * 1.0e6_r8
00914       dataGpr = dataGpr/budg_ns
00915 
00916       if (iam /= 0) return
00917    endif
00918 
00919    ! ---------------------------------------------------------
00920    ! ---- detail atm budgets and breakdown into components ---
00921    ! ---------------------------------------------------------
00922 
00923    if (plev >= 3) then
00924    do ic = 1,2
00925       if (ic == 1) then
00926          ica = c_atm_ar
00927          icl = c_lnd_ar
00928          icn = c_inh_ar
00929          ics = c_ish_ar
00930          ico = c_ocn_ar
00931          str = "ATM_to_CPL"
00932       elseif (ic == 2) then
00933          ica = c_atm_as
00934          icl = c_lnd_as
00935          icn = c_inh_as
00936          ics = c_ish_as
00937          ico = c_ocn_as
00938          str = "CPL_TO_ATM"
00939       else
00940          call shr_sys_abort(subname//' ERROR in ic index code 411')
00941       endif
00942 
00943       write(logunit,*) ' '
00944       write(logunit,FAH) subname,trim(str)//' AREA BUDGET (m2/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
00945       write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM*  '
00946       do if = f_a, f_h-1
00947          write(logunit,FA1)    fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
00948                    dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
00949                                          dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
00950                    dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip) 
00951       enddo
00952 
00953       write(logunit,*) ' '
00954       write(logunit,FAH) subname,trim(str)//' HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
00955       write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM*  '
00956       do if = f_h, f_w-1
00957          write(logunit,FA1)    fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
00958                    dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
00959                                          dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
00960                    dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip) 
00961       enddo
00962       write(logunit,FA1)    '   *SUM*',sum(dataGpr(f_h:f_w-1,ica,ip)),sum(dataGpr(f_h:f_w-1,icl,ip)), &
00963          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)), &
00964                                        sum(dataGpr(f_h:f_w-1,ica,ip))+sum(dataGpr(f_h:f_w-1,icl,ip))+ &
00965          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)) 
00966 
00967       write(logunit,*) ' '
00968       write(logunit,FAH) subname,trim(str)//' WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
00969       write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM*  '
00970       do if = f_w, f_size
00971          write(logunit,FA1)    fname(if),dataGpr(if,ica,ip),dataGpr(if,icl,ip), &
00972                    dataGpr(if,icn,ip),dataGpr(if,ics,ip),dataGpr(if,ico,ip), &
00973                                          dataGpr(if,ica,ip)+dataGpr(if,icl,ip)+ &
00974                    dataGpr(if,icn,ip)+dataGpr(if,ics,ip)+dataGpr(if,ico,ip) 
00975       enddo
00976       write(logunit,FA1)    '   *SUM*',sum(dataGpr(f_w:f_size,ica,ip)),sum(dataGpr(f_w:f_size,icl,ip)), &
00977          sum(dataGpr(f_w:f_size,icn,ip)),sum(dataGpr(f_w:f_size,ics,ip)),sum(dataGpr(f_w:f_size,ico,ip)), &
00978                                        sum(dataGpr(f_w:f_size,ica,ip))+sum(dataGpr(f_w:f_size,icl,ip))+ &
00979          sum(dataGpr(f_w:f_size,icn,ip))+sum(dataGpr(f_w:f_size,ics,ip))+sum(dataGpr(f_w:f_size,ico,ip)) 
00980    enddo
00981    endif   ! plev
00982 
00983    ! ---------------------------------------------------------
00984    ! ---- detail lnd/ocn/ice component budgets ----
00985    ! ---------------------------------------------------------
00986 
00987    if (plev >= 2) then
00988    do ic = 1,4
00989       if (ic == 1) then
00990          icar = c_lnd_ar
00991          icxs = c_lnd_ls
00992          icxr = c_lnd_lr
00993          icas = c_lnd_as
00994          str = "LND"
00995       elseif (ic == 2) then
00996          icar = c_ocn_ar
00997          icxs = c_ocn_os
00998          icxr = c_ocn_or
00999          icas = c_ocn_as
01000          str = "OCN"
01001       elseif (ic == 3) then
01002          icar = c_inh_ar
01003          icxs = c_inh_is
01004          icxr = c_inh_ir
01005          icas = c_inh_as
01006          str = "ICE_NH"
01007       elseif (ic == 4) then
01008          icar = c_ish_ar
01009          icxs = c_ish_is
01010          icxr = c_ish_ir
01011          icas = c_ish_as
01012          str = "ICE_SH"
01013       else
01014          call shr_sys_abort(subname//' ERROR in ic index code 412')
01015       endif
01016 
01017       write(logunit,*) ' '
01018       write(logunit,FAH) subname,trim(str)//' HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
01019       write(logunit,FA0) cname(icar),cname(icxs),cname(icxr),cname(icas),' *SUM*  '
01020       do if = f_h, f_w-1
01021          write(logunit,FA1)    fname(if),-dataGpr(if,icar,ip),dataGpr(if,icxs,ip), &
01022                                           dataGpr(if,icxr,ip),-dataGpr(if,icas,ip), &
01023                                          -dataGpr(if,icar,ip)+dataGpr(if,icxs,ip)+ &
01024                                           dataGpr(if,icxr,ip)-dataGpr(if,icas,ip)
01025       enddo
01026       write(logunit,FA1)    '   *SUM*',-sum(dataGpr(f_h:f_w-1,icar,ip)),sum(dataGpr(f_h:f_w-1,icxs,ip)), &
01027                                        sum(dataGpr(f_h:f_w-1,icxr,ip)),-sum(dataGpr(f_h:f_w-1,icas,ip)), &
01028                                        -sum(dataGpr(f_h:f_w-1,icar,ip))+sum(dataGpr(f_h:f_w-1,icxs,ip))+ &
01029                                        sum(dataGpr(f_h:f_w-1,icxr,ip))-sum(dataGpr(f_h:f_w-1,icas,ip))
01030 
01031       write(logunit,*) ' '
01032       write(logunit,FAH) subname,trim(str)//' WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
01033       write(logunit,FA0) cname(icar),cname(icxs),cname(icxr),cname(icas),' *SUM*  '
01034       do if = f_w, f_size
01035          write(logunit,FA1)    fname(if),-dataGpr(if,icar,ip),dataGpr(if,icxs,ip), &
01036                                          dataGpr(if,icxr,ip),-dataGpr(if,icas,ip), &
01037                                          -dataGpr(if,icar,ip)+dataGpr(if,icxs,ip)+ &
01038                                          dataGpr(if,icxr,ip)-dataGpr(if,icas,ip)
01039       enddo
01040       write(logunit,FA1)    '   *SUM*',-sum(dataGpr(f_w:f_size,icar,ip)),sum(dataGpr(f_w:f_size,icxs,ip)), &
01041                                        sum(dataGpr(f_w:f_size,icxr,ip)),-sum(dataGpr(f_w:f_size,icas,ip)), &
01042                                        -sum(dataGpr(f_w:f_size,icar,ip))+sum(dataGpr(f_w:f_size,icxs,ip))+ &
01043                                        sum(dataGpr(f_w:f_size,icxr,ip))-sum(dataGpr(f_w:f_size,icas,ip))
01044 
01045    enddo
01046    endif   ! plev
01047 
01048    ! ---------------------------------------------------------
01049    ! ---- net summary budgets ----
01050    ! ---------------------------------------------------------
01051 
01052    if (plev >= 1) then
01053 
01054       write(logunit,*) ' '
01055       write(logunit,FAH) subname,'NET AREA BUDGET (m2/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
01056       write(logunit,FA0) '     atm','     lnd','     ocn','  ice nh','  ice sh',' *SUM*  '
01057       do if = 1,f_h-1
01058          write(logunit,FA1)    fname(if),dataGpr(if,c_atm_ar,ip), &
01059                                          dataGpr(if,c_lnd_lr,ip), &
01060                                          dataGpr(if,c_ocn_or,ip), &
01061                                          dataGpr(if,c_inh_ir,ip), &
01062                                          dataGpr(if,c_ish_ir,ip), &
01063                                          dataGpr(if,c_atm_ar,ip)+ &
01064                                          dataGpr(if,c_lnd_lr,ip)+ &
01065                                          dataGpr(if,c_ocn_or,ip)+ &
01066                                          dataGpr(if,c_inh_ir,ip)+ &
01067                                          dataGpr(if,c_ish_ir,ip)
01068       enddo
01069 
01070       write(logunit,*) ' '
01071       write(logunit,FAH) subname,'NET HEAT BUDGET (W/m2): period = ',trim(pname(ip)),': date = ',cdate,sec
01072       write(logunit,FA0) '     atm','     lnd','     ocn','  ice nh','  ice sh',' *SUM*  '
01073       do if = f_h, f_w-1
01074          write(logunit,FA1)    fname(if),dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip), &
01075                                          dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip), &
01076                                          dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip), &
01077                                          dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip), &
01078                                          dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip), &
01079                                          dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip)+ &
01080                                          dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip)+ &
01081                                          dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip)+ &
01082                                          dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip)+ &
01083                                          dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip)
01084       enddo
01085       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)), &
01086                                     sum(dataGpr(f_h:f_w-1,c_lnd_lr,ip))+sum(dataGpr(f_h:f_w-1,c_lnd_ls,ip)), &
01087                                     sum(dataGpr(f_h:f_w-1,c_ocn_or,ip))+sum(dataGpr(f_h:f_w-1,c_ocn_os,ip)), &
01088                                     sum(dataGpr(f_h:f_w-1,c_inh_ir,ip))+sum(dataGpr(f_h:f_w-1,c_inh_is,ip)), &
01089                                     sum(dataGpr(f_h:f_w-1,c_ish_ir,ip))+sum(dataGpr(f_h:f_w-1,c_ish_is,ip)), &
01090                                     sum(dataGpr(f_h:f_w-1,c_atm_ar,ip))+sum(dataGpr(f_h:f_w-1,c_atm_as,ip))+ &
01091                                     sum(dataGpr(f_h:f_w-1,c_lnd_lr,ip))+sum(dataGpr(f_h:f_w-1,c_lnd_ls,ip))+ &
01092                                     sum(dataGpr(f_h:f_w-1,c_ocn_or,ip))+sum(dataGpr(f_h:f_w-1,c_ocn_os,ip))+ &
01093                                     sum(dataGpr(f_h:f_w-1,c_inh_ir,ip))+sum(dataGpr(f_h:f_w-1,c_inh_is,ip))+ &
01094                                     sum(dataGpr(f_h:f_w-1,c_ish_ir,ip))+sum(dataGpr(f_h:f_w-1,c_ish_is,ip))
01095 
01096       write(logunit,*) ' '
01097       write(logunit,FAH) subname,'NET WATER BUDGET (kg/m2s*1e6): period = ',trim(pname(ip)),': date = ',cdate,sec
01098       write(logunit,FA0) '     atm','     lnd','     ocn','  ice nh','  ice sh',' *SUM*  '
01099       do if = f_w, f_size
01100          write(logunit,FA1)    fname(if),dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip), &
01101                                          dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip), &
01102                                          dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip), &
01103                                          dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip), &
01104                                          dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip), &
01105                                          dataGpr(if,c_atm_ar,ip)+dataGpr(if,c_atm_as,ip)+ &
01106                                          dataGpr(if,c_lnd_lr,ip)+dataGpr(if,c_lnd_ls,ip)+ &
01107                                          dataGpr(if,c_ocn_or,ip)+dataGpr(if,c_ocn_os,ip)+ &
01108                                          dataGpr(if,c_inh_ir,ip)+dataGpr(if,c_inh_is,ip)+ &
01109                                          dataGpr(if,c_ish_ir,ip)+dataGpr(if,c_ish_is,ip)
01110       enddo
01111       write(logunit,FA1) '   *SUM*',sum(dataGpr(f_w:f_size,c_atm_ar,ip))+sum(dataGpr(f_w:f_size,c_atm_as,ip)), &
01112                                     sum(dataGpr(f_w:f_size,c_lnd_lr,ip))+sum(dataGpr(f_w:f_size,c_lnd_ls,ip)), &
01113                                     sum(dataGpr(f_w:f_size,c_ocn_or,ip))+sum(dataGpr(f_w:f_size,c_ocn_os,ip)), &
01114                                     sum(dataGpr(f_w:f_size,c_inh_ir,ip))+sum(dataGpr(f_w:f_size,c_inh_is,ip)), &
01115                                     sum(dataGpr(f_w:f_size,c_ish_ir,ip))+sum(dataGpr(f_w:f_size,c_ish_is,ip)), &
01116                                     sum(dataGpr(f_w:f_size,c_atm_ar,ip))+sum(dataGpr(f_w:f_size,c_atm_as,ip))+ &
01117                                     sum(dataGpr(f_w:f_size,c_lnd_lr,ip))+sum(dataGpr(f_w:f_size,c_lnd_ls,ip))+ &
01118                                     sum(dataGpr(f_w:f_size,c_ocn_or,ip))+sum(dataGpr(f_w:f_size,c_ocn_os,ip))+ &
01119                                     sum(dataGpr(f_w:f_size,c_inh_ir,ip))+sum(dataGpr(f_w:f_size,c_inh_is,ip))+ &
01120                                     sum(dataGpr(f_w:f_size,c_ish_ir,ip))+sum(dataGpr(f_w:f_size,c_ish_is,ip))
01121 
01122    endif
01123 
01124    write(logunit,*) ' '
01125 ! ---- doprint ---- doprint ---- doprint ----
01126    endif  ! plev > 0
01127    enddo  ! ip = 1,p_size
01128 
01129 end subroutine seq_diag_print_mct
01130 
01131 !===============================================================================
01132 !BOP ===========================================================================
01133 !
01134 ! !IROUTINE: seq_diag_avect_mct - print global budget diagnostics
01135 !
01136 ! !DESCRIPTION:
01137 !   Print global diagnostics for AV/ID.
01138 !
01139 ! !REVISION HISTORY:
01140 !
01141 ! !INTERFACE: ------------------------------------------------------------------
01142 
01143 SUBROUTINE seq_diag_avect_mct(cdata,AV,comment)
01144 
01145    use seq_infodata_mod
01146 
01147    implicit none
01148 
01149 ! !INPUT/OUTPUT PARAMETERS:
01150 
01151    type(seq_cdata) , intent(in) :: cdata
01152    type(mct_aVect) , intent(in) :: AV
01153    character(len=*), intent(in), optional :: comment
01154 
01155 !EOP
01156 
01157    !--- local ---
01158    type(mct_gGrid)   ,pointer :: dom
01159    type(seq_infodata_type),pointer :: infodata
01160    type(mct_gsMap)   ,pointer :: gsmap
01161    integer(in)      :: ID
01162    logical          :: bfbflag
01163    integer(in)      :: n,k         ! counters
01164    integer(in)      :: npts,nptsg  ! number of local/global pts in AV
01165    integer(in)      :: kflds       ! number of fields in AV
01166    real(r8),pointer :: sumbuf (:)  ! sum buffer
01167    real(r8),pointer :: minbuf (:)  ! min buffer
01168    real(r8),pointer :: maxbuf (:)  ! max buffer
01169    real(r8),pointer :: sumbufg(:)  ! sum buffer reduced
01170    real(r8),pointer :: minbufg(:)  ! min buffer reduced
01171    real(r8),pointer :: maxbufg(:)  ! max buffer reduced
01172    integer(i8),pointer :: isumbuf (:) ! integer local sum
01173    integer(i8),pointer :: isumbufg(:) ! integer global sum
01174    integer(i8)      :: ihuge       ! huge
01175    integer(in)      :: mpicom      ! mpi comm
01176    integer(in)      :: iam         ! pe number
01177    integer(in)      :: km,ka       ! field indices
01178    integer(in)      :: ns          ! size of local AV
01179    integer(in)      :: rcode       ! error code
01180    real(r8),pointer :: weight(:)   ! weight
01181    type(mct_string) :: mstring     ! mct char type
01182    character(CL)    :: lcomment    ! should be long enough
01183    character(CL)    :: itemc       ! string converted to char
01184 
01185    type(mct_avect)  :: AV1         ! local avect with one field
01186    type(mct_avect)  :: AVr1        ! avect on root with one field
01187    type(mct_avect)  :: AVr2        ! avect on root with one field
01188 
01189    !----- formats -----
01190    character(*),parameter :: subName = '(seq_diag_avect_mct) '
01191    character(*),parameter :: F00   = "('(seq_diag_avect_mct) ',4a)"
01192 
01193 !-------------------------------------------------------------------------------
01194 ! print instantaneous budget data
01195 !-------------------------------------------------------------------------------
01196 
01197    call seq_cdata_setptrs(cdata,ID=ID,infodata=infodata,dom=dom,gsmap=gsmap)
01198    call seq_comm_setptrs(ID,mpicom=mpicom,iam=iam)
01199    call seq_infodata_GetData(infodata,bfbflag=bfbflag)
01200 
01201    lcomment = ''
01202    if (present(comment)) then
01203       lcomment=trim(comment)
01204    endif
01205 
01206    ns = mct_aVect_lsize(AV)
01207    npts = mct_aVect_lsize(dom%data)
01208    if (ns /= npts) call shr_sys_abort(trim(subname)//' ERROR: size of AV,dom')
01209    km = mct_aVect_indexRA(dom%data,'mask')
01210    ka = mct_aVect_indexRA(dom%data,afldname)
01211    kflds = mct_aVect_nRattr(AV)
01212    allocate(sumbuf(kflds),sumbufg(kflds))
01213 
01214    sumbuf =       0.0_r8
01215 
01216    if (bfbflag) then
01217 
01218       npts = mct_aVect_lsize(AV)
01219       allocate(weight(npts))
01220       weight(:) = 1.0_r8 
01221       do n = 1,npts
01222          if (dom%data%rAttr(km,n) <= 1.0e-06_R8) then
01223             weight(n) = 0.0_r8
01224          else
01225             weight(n) = dom%data%rAttr(ka,n)*shr_const_rearth*shr_const_rearth
01226          endif
01227       enddo
01228 
01229       allocate(maxbuf(kflds),maxbufg(kflds))
01230       maxbuf = 0.0_r8
01231 
01232       do n = 1,npts
01233       do k = 1,kflds
01234          if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
01235              AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
01236              maxbuf(k) = max(maxbuf(k),abs(AV%rAttr(k,n)*weight(n)))
01237          endif
01238       enddo
01239       enddo
01240 
01241       call shr_mpi_max(maxbuf,maxbufg,mpicom,subname,all=.true.)
01242       call shr_mpi_sum(npts,nptsg,mpicom,subname,all=.true.)
01243 
01244       do k = 1,kflds
01245          if (maxbufg(k) < 1000.0*TINY(maxbufg(k)) .or. &
01246              maxbufg(k) > HUGE(maxbufg(k))/(2.0_r8*nptsg)) then
01247             maxbufg(k) = 0.0_r8
01248          else
01249             maxbufg(k) = (1.1_r8) * maxbufg(k) * nptsg
01250          endif
01251       enddo
01252 
01253       allocate(isumbuf(kflds),isumbufg(kflds))
01254       isumbuf = 0
01255       ihuge = HUGE(isumbuf)
01256 
01257       do n = 1,npts
01258       do k = 1,kflds
01259          if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
01260              AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
01261              if (abs(maxbufg(k)) > 1000.0_r8 * TINY(maxbufg)) then
01262                 isumbuf(k) = isumbuf(k) + int((AV%rAttr(k,n)*weight(n)/maxbufg(k))*ihuge,i8)
01263              endif
01264          endif
01265       enddo
01266       enddo
01267 
01268       call shr_mpi_sum(isumbuf,isumbufg,mpicom,subname)
01269 
01270       do k = 1,kflds
01271          sumbufg(k) = isumbufg(k)*maxbufg(k)/ihuge
01272       enddo
01273 
01274       deallocate(weight)
01275       deallocate(maxbuf,maxbufg)
01276       deallocate(isumbuf,isumbufg)
01277 
01278 #if (1 == 0)
01279 ! OLD VERSION
01280       call mct_aVect_init(AV1,rList='varf1',lsize=ns)
01281 
01282       AV1%rAttr(1,1:ns) = dom%data%rAttr(km,1:ns)           ! mask = AVr1
01283       call mct_aVect_gather(AV1,AVr1,gsmap,0,mpicom,rcode)
01284       AV1%rAttr(1,1:ns) = dom%data%rAttr(ka,1:ns)           ! area = AVr2
01285       call mct_aVect_gather(AV1,AVr2,gsmap,0,mpicom,rcode)
01286 
01287       ! --- compute weight on root pe
01288 
01289       if (iam == 0) then
01290          npts = mct_aVect_lsize(AVr1)
01291          allocate(weight(npts))
01292          weight(:) = 1.0_r8 
01293          do n = 1,npts
01294             if (AVr1%rAttr(1,n) <= 1.0e-06_R8) then
01295                weight(n) = 0.0_r8
01296             else
01297                weight(n) = AVr2%rAttr(1,n)*shr_const_rearth*shr_const_rearth
01298             endif
01299          enddo
01300 !         write(logunit,*) trim(subname),'tcx1 ',ns,npts,km,ka,kflds,minval(AVr2%rAttr),maxval(AVr2%rAttr)
01301 !         write(logunit,*) trim(subname),'tcx2 ',size(AVr1%rAttr),size(AVr2%rAttr),minval(AVr1%rAttr),maxval(AVr1%rAttr)
01302       endif
01303 
01304       ! --- gather and compute stats one field at a time (to minimize memory)
01305 
01306       do k = 1,kflds
01307          AV1%rAttr(1,1:ns) = AV%rAttr(k,1:ns)           ! fld k = AVr2
01308          call mct_aVect_gather(AV1,AVr2,gsmap,0,mpicom,rcode)
01309          if (iam == 0) then
01310             npts = mct_aVect_lsize(AVr2)
01311 !            write(logunit,*) trim(subname),'tcx3 ',ns,npts,k,size(AVr2%rAttr),minval(AVr2%rAttr),maxval(AVr2%rAttr)
01312 !            write(logunit,*) trim(subname),'tcx4 ',size(weight),minval(weight),maxval(weight)
01313             do n = 1,npts
01314                if (AVr2%rAttr(1,n) > 1.01_r8*shr_const_spval .or. &
01315                    AVr2%rAttr(1,n) < 0.99_r8*shr_const_spval) then
01316                    sumbuf(k) = sumbuf(k) + AVr2%rAttr(1,n)*weight(n)
01317                endif
01318             enddo
01319          endif
01320       enddo
01321 
01322       ! --- local copy, relevant only on root pe
01323       sumbufg = sumbuf
01324 
01325       call mct_avect_clean(AV1)
01326       if (iam == 0) then
01327          call mct_aVect_clean(AVr1)
01328          call mct_aVect_clean(AVr2)
01329       endif
01330 #endif
01331 
01332    else
01333 
01334       npts = mct_aVect_lsize(AV)
01335       allocate(weight(npts))
01336       weight(:) = 1.0_r8 
01337       do n = 1,npts
01338          if (dom%data%rAttr(km,n) <= 1.0e-06_R8) then
01339             weight(n) = 0.0_r8
01340          else
01341             weight(n) = dom%data%rAttr(ka,n)*shr_const_rearth*shr_const_rearth
01342          endif
01343       enddo
01344 
01345       do n = 1,npts
01346       do k = 1,kflds
01347          if (AV%rAttr(k,n) > 1.01_r8*shr_const_spval .or. &
01348              AV%rAttr(k,n) < 0.99_r8*shr_const_spval) then
01349              sumbuf(k) = sumbuf(k) + AV%rAttr(k,n)*weight(n)
01350          endif
01351       enddo
01352       enddo
01353 
01354       !--- global reduction ---
01355       call shr_mpi_sum(sumbuf,sumbufg,mpicom,subname)
01356 
01357       deallocate(weight)
01358 
01359    endif
01360 
01361    if (iam == 0) then
01362 !      write(logunit,*) 'sdAV: *** writing ',trim(lcomment),': k fld min/max/sum ***'
01363       do k = 1,kflds
01364          call mct_aVect_getRList(mstring,k,AV)
01365          itemc = mct_string_toChar(mstring)
01366          call mct_string_clean(mstring)
01367          if (len_trim(lcomment) > 0) then
01368             write(logunit,100) 'xxx','sorr',k,sumbufg(k),trim(lcomment),trim(itemc)
01369          else
01370             write(logunit,101) 'xxx','sorr',k,sumbufg(k),trim(itemc)
01371          endif
01372       enddo
01373       call shr_sys_flush(logunit)
01374    endif
01375 
01376    deallocate(sumbuf,sumbufg)
01377 
01378 100  format('comm_diag ',a3,1x,a4,1x,i3,es26.19,1x,a,1x,a)
01379 101  format('comm_diag ',a3,1x,a4,1x,i3,es26.19,1x,a)
01380 
01381 end subroutine seq_diag_avect_mct
01382 
01383 !===============================================================================
01384 !BOP ===========================================================================
01385 !
01386 ! !IROUTINE: seq_diag_avdiff_mct - print global budget diagnostics
01387 !
01388 ! !DESCRIPTION:
01389 !   Print global diagnostics for AV/ID.
01390 !
01391 ! !REVISION HISTORY:
01392 !
01393 ! !INTERFACE: ------------------------------------------------------------------
01394 
01395 SUBROUTINE seq_diag_avdiff_mct(AV1,AV2,ID,comment)
01396 
01397    implicit none
01398 
01399 ! !INPUT/OUTPUT PARAMETERS:
01400 
01401    type(mct_aVect) , intent(in) :: AV1
01402    type(mct_aVect) , intent(in) :: AV2
01403    integer         , intent(in) :: ID
01404    character(len=*), intent(in), optional :: comment
01405 
01406 !EOP
01407 
01408    !--- local ---
01409    integer(in)      :: n,k,n1,k1,n2,k2         ! counters
01410    integer(in)      :: iam         ! pe number
01411    integer(in)      :: cnt         ! counter
01412    real(r8)         :: adiff,rdiff ! diff values
01413    type(mct_string) :: mstring     ! mct char type
01414    character(len=64):: lcomment    ! should be long enough
01415 
01416    !----- formats -----
01417    character(*),parameter :: subName = '(seq_diag_avdiff_mct) '
01418    character(*),parameter :: F00   = "('(seq_diag_avdiff_mct) ',4a)"
01419 
01420 !-------------------------------------------------------------------------------
01421 ! print instantaneous budget data
01422 !-------------------------------------------------------------------------------
01423 
01424    call seq_comm_setptrs(ID,iam=iam)
01425 
01426    lcomment = ''
01427    if (present(comment)) then
01428       lcomment=trim(comment)
01429    endif
01430 
01431    n1 = mct_aVect_lsize(AV1)
01432    k1 = mct_aVect_nRattr(AV1)
01433    n2 = mct_aVect_lsize(AV2)
01434    k2 = mct_aVect_nRattr(AV2)
01435 
01436    if (n1 /= n2 .or. k1 /= k2) then
01437       write(s_logunit,*) subname,trim(lcomment),' AV sizes different ',n1,n2,k1,k2
01438       return
01439    endif
01440 
01441    do k = 1,k1
01442       cnt = 0
01443       adiff = 0.
01444       rdiff = 0.
01445       do n = 1,n1
01446          if (AV1%rAttr(k,n) /= AV2%rAttr(k,n)) then
01447             cnt = cnt + 1
01448             adiff = max(adiff, abs(AV1%rAttr(k,n)-AV2%rAttr(k,n)))
01449             rdiff = max(rdiff, abs(AV1%rAttr(k,n)-AV2%rAttr(k,n))/(abs(AV1%rAttr(k,n))+abs(AV2%rAttr(k,n))))
01450          endif
01451       enddo
01452       if (cnt > 0) then
01453          call mct_aVect_getRList(mstring,k,AV1)
01454          write(s_logunit,*) subname,trim(lcomment),' AVs fld k diff ', &
01455             iam,mct_string_toChar(mstring),cnt,adiff,rdiff, &
01456             minval(AV1%rAttr(k,:)),minval(AV1%rAttr(k,:)), &
01457             maxval(AV1%rAttr(k,:)),maxval(AV2%rAttr(k,:))
01458          call mct_string_clean(mstring)       
01459       endif
01460    enddo
01461 
01462 end subroutine seq_diag_avdiff_mct
01463 
01464 !===============================================================================
01465 end module seq_diag_mct