CESM_Driver: seq_frac_mct.F90 Source File

00001 !===============================================================================
00002 ! SVN $Id: seq_frac_mct.F90 19875 2009-12-14 21:30:58Z tcraig $
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_frac_mct.F90 $
00004 !===============================================================================
00005 !BOP ===========================================================================
00006 !
00007 ! !MODULE: seq_frac_mct -- handles surface fractions.
00008 !
00009 !  Fraction Notes: tcraig, august 2008
00010 !
00011 !  the fractions fields are now afrac, ifrac, ofrac, lfrac, and lfrin.
00012 !    afrac = fraction of atm on a grid
00013 !    lfrac = fraction of lnd on a grid
00014 !    ifrac = fraction of ice on a grid
00015 !    ofrac = fraction of ocn on a grid
00016 !    lfrin = land fraction defined by the land model
00017 !    ifrad = fraction of ocn on a grid at last radiation time
00018 !    ofrad = fraction of ice on a grid at last radiation time
00019 !      afrac, lfrac, ifrac, and ofrac are the self-consistent values in the 
00020 !      system.  lfrin is the fraction on the land grid and is allowed to
00021 !      vary from the self-consistent value as descibed below.  ifrad
00022 !      and ofrad are needed for the swnet calculation.
00023 !  the fractions fields are defined for each grid in the fraction bundles as 
00024 !    needed as follows.
00025 !    character(*),parameter :: fraclist_a = 'afrac:ifrac:ofrac:lfrac:lfrin'
00026 !    character(*),parameter :: fraclist_o = 'afrac:ifrac:ofrac:ifrad:ofrad'
00027 !    character(*),parameter :: fraclist_i = 'afrac:ifrac:ofrac'
00028 !    character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin'
00029 !    character(*),parameter :: fraclist_g = 'gfrac'
00030 !
00031 !  we assume ocean and ice are on the same grids, same masks
00032 !  we assume ocn2atm and ice2atm are masked maps
00033 !  we assume lnd2atm is a global map
00034 !  we assume that the ice fraction evolves in time but that
00035 !    the land model fraction does not.  the ocean fraction then
00036 !    is just the complement of the ice fraction over the region
00037 !    of the ocean/ice mask.
00038 !  we assume that component domains are filled with the total 
00039 !    potential mask/fraction on that grid, but that the fractions
00040 !    sent at run time are always the relative fraction covered.
00041 !    for example, if an ice cell can be up to 50% covered in 
00042 !    ice and 50% land, then the ice domain should have a fraction
00043 !    value of 0.5 at that grid cell.  at run time though, the ice
00044 !    fraction will be between 0.0 and 1.0 meaning that grid cells
00045 !    is covered with between 0.0 and 0.5 by ice.  the "relative" fractions 
00046 !    sent at run-time are corrected by the model to be total fractions
00047 !    such that
00048 !  in general, on every grid,
00049 !              fractions_*(afrac) = 1.0
00050 !              fractions_*(ifrac) + fractions_*(ofrac) + fractions_*(lfrac) = 1.0
00051 !  where fractions_* are a bundle of fractions on a particular grid and 
00052 !    *frac (ie afrac) is the fraction of a particular component in the bundle.
00053 !
00054 !  the fractions are computed fundamentally as follows (although the 
00055 !    detailed implementation might be slightly different)
00056 !  initialization (frac_init):
00057 !    afrac is set on all grids
00058 !      fractions_a(afrac) = 1.0
00059 !      fractions_o(afrac) = mapa2o(fractions_a(afrac))
00060 !      fractions_i(afrac) = mapa2i(fractions_a(afrac))
00061 !      fractions_l(afrac) = mapa2l(fractions_a(afrac))
00062 !    initially assume ifrac on all grids is zero
00063 !      fractions_*(ifrac) = 0.0
00064 !    fractions/masks provided by surface components
00065 !      fractions_o(ofrac) = dom_o(frac)  ! ocean "mask"
00066 !      fractions_l(lfrin) = dom_l(frac)  ! land model fraction
00067 !    then mapped to the atm model
00068 !      fractions_a(ofrac) = mapo2a(fractions_o(ofrac))
00069 !      fractions_a(lfrin) = mapl2a(fractions_l(lfrin))
00070 !    and a few things are then derived
00071 !      fractions_a(lfrac) = 1.0 - fractions_a(ofrac)
00072 !        this is truncated to zero for very small values (< 0.001)
00073 !        to attempt to preserve non-land gridcells.
00074 !      fractions_l(lfrac) = mapa2l(fractions_a(lfrac))
00075 !    one final term is computed
00076 !      dom_a(ascale) = fractions_a(lfrac)/fractions_a(lfrin)
00077 !      dom_l(ascale) = mapa2l(dom_a(ascale))
00078 !        these are used to correct land fluxes in budgets and lnd2rtm coupling
00079 !        and are particularly important when the land model is running on 
00080 !        a different grid than the atm model.  in the old system, this term
00081 !        was treated as effectively 1.0 since there was always a check that
00082 !        fractions_a(lfrac) ~ fractions_a(lfrin), namely that the land model
00083 !        provided a land frac that complemented the ocean grid.  this is
00084 !        no longer a requirement in this new system and as a result, the
00085 !        ascale term can be thought of as a rescaling of the land fractions
00086 !        in the land model to be exactly complementary to the ocean model
00087 !        on whatever grid it may be running.
00088 !  run-time (frac_set):
00089 !    update fractions on ice grid
00090 !      fractions_i(ifrac) = i2x_i(Si_ifrac)  ! ice frac from ice model
00091 !      fractions_i(ofrac) = 1.0 - fractions_i(ifrac)
00092 !        note: the relative fractions are corrected to total fractions
00093 !      fractions_o(ifrac) = mapi2o(fractions_i(ifrac))
00094 !      fractions_o(ofrac) = mapi2o(fractions_i(ofrac))
00095 !      fractions_a(ifrac) = mapi2a(fractions_i(ifrac))
00096 !      fractions_a(ofrac) = mapi2a(fractions_i(ofrac))
00097 ! 
00098 !  fractions used in merging are as follows
00099 !    mrg_x2a uses fractions_a(lfrac,ofrac,ifrac)
00100 !    mrg_x2o needs to use fractions_o(ofrac,ifrac) normalized to one
00101 !      normalization happens in mrg routine
00102 !
00103 !  fraction corrections in mapping are as follows
00104 !    mapo2a uses *fractions_o(ofrac) and /fractions_a(ofrac)
00105 !    mapi2a uses *fractions_i(ifrac) and /fractions_a(ifrac)
00106 !    mapl2a uses *fractions_l(lfrin) and /fractions_a(lfrin)
00107 !    mapa2* should use *fractions_a(afrac) and /fractions_*(afrac) but this
00108 !      has been defered since the ratio always close to 1.0
00109 !
00110 !  budgets use the standard afrac, ofrac, ifrac, and lfrac to compute
00111 !    quantities except in the land budget which uses lfrin multiplied
00112 !    by the scale factor, dom_l(ascale) to compute budgets.
00113 !
00114 !  fraction and domain checks
00115 !    initialization:
00116 !      dom_i = mapo2i(dom_o)  ! lat, lon, mask, area
00117 !      where fractions_a(lfrac) > 0.0, fractions_a(lfrin) is also > 0.0
00118 !         this ensures the land will provide data everywhere the atm needs it
00119 !         and allows the land frac to be subtlely different from the
00120 !         land fraction specified in the atm.
00121 !      dom_a = mapl2a(dom_l)  ! if atm/lnd same grids
00122 !      dom_a = mapo2a(dom_o)  ! if atm/ocn same grids
00123 !      dom_a = mapi2a(dom_i)  ! if atm/ocn same grids
00124 !      0.0-eps < fractions_*(*) < 1.0+eps
00125 !      fractions_l(lfrin) = fractions_l(lfrac) 
00126 !        only if atm/lnd same grids (but this is not formally required)
00127 !        this is needed until dom_l(ascale) is sent to the land model
00128 !        as an additional field for use in l2r mapping.
00129 !    run time:
00130 !      fractions_a(lfrac) + fractions_a(ofrac) + fractions_a(ifrac) ~ 1.0
00131 !      0.0-eps < fractions_*(*) < 1.0+eps
00132 !    
00133 !!
00134 !
00135 ! !REVISION HISTORY:
00136 !    2007-may-07 - M. Vertenstein - initial port to cpl7.
00137 !
00138 ! !INTERFACE: ------------------------------------------------------------------
00139                                                                                       
00140 module seq_frac_mct
00141 
00142 ! !USES:
00143 
00144   use shr_kind_mod     , only: R8 => SHR_KIND_R8
00145   use shr_sys_mod
00146   use shr_const_mod
00147 
00148   use mct_mod
00149   use seq_cdata_mod
00150   use map_iceocn_mct
00151   use map_atmice_mct
00152   use map_atmocn_mct
00153   use map_atmlnd_mct
00154   use seq_comm_mct, only: logunit, loglevel
00155 
00156   implicit none
00157   private
00158   save
00159 
00160 ! !PUBLIC TYPES:
00161 
00162   ! none
00163 
00164 ! !PUBLIC MEMBER FUNCTIONS:
00165 
00166   public seq_frac_init
00167   public seq_frac_set
00168 
00169 ! !PUBLIC DATA MEMBERS:
00170 
00171 !EOP
00172 
00173 ! !LOCAL DATA
00174 
00175 #ifdef CPP_VECTOR
00176   logical :: usevector = .true.
00177 #else
00178   logical :: usevector = .false.
00179 #endif
00180   
00181 #ifdef SYSUNICOS
00182   logical :: usealltoall = .true.
00183 #else
00184   logical :: usealltoall = .false.
00185 #endif
00186 
00187   integer, private :: seq_frac_debug = 1
00188   logical, private :: seq_frac_abort = .true.
00189   logical, private :: seq_frac_dead
00190 
00191   !--- standard ---
00192   real(r8),parameter :: eps_fracsum = 1.0e-02   ! allowed error in sum of fracs 
00193   real(r8),parameter :: eps_fracval = 1.0e-02   ! allowed error in any frac +- 0,1
00194   real(r8),parameter :: eps_fraclim = 1.0e-03   ! truncation limit in fractions_a(lfrac)
00195   logical ,parameter :: atm_frac_correct = .false. ! turn on frac correction on atm grid
00196   !--- standard plus atm fraction consistency ---
00197 !  real(r8),parameter :: eps_fracsum = 1.0e-12   ! allowed error in sum of fracs 
00198 !  real(r8),parameter :: eps_fracval = 1.0e-02   ! allowed error in any frac +- 0,1
00199 !  real(r8),parameter :: eps_fraclim = 1.0e-03   ! truncation limit in fractions_a(lfrac)
00200 !  logical ,parameter :: atm_frac_correct = .true. ! turn on frac correction on atm grid
00201   !--- unconstrained and area conserving? ---
00202 !  real(r8),parameter :: eps_fracsum = 1.0e-12   ! allowed error in sum of fracs 
00203 !  real(r8),parameter :: eps_fracval = 1.0e-02   ! allowed error in any frac +- 0,1
00204 !  real(r8),parameter :: eps_fraclim = 1.0e-20   ! truncation limit in fractions_a(lfrac)
00205 !  logical ,parameter :: atm_frac_correct = .true. ! turn on frac correction on atm grid
00206 
00207   private seq_frac_check
00208 
00209 !===============================================================================
00210 contains
00211 !===============================================================================
00212 
00213 !===============================================================================
00214 !BOP ===========================================================================
00215 !
00216 ! !IROUTINE: seq_frac_init 
00217 !
00218 ! !DESCRIPTION:
00219 !    Initialize fraction attribute vectors and necessary ocn/ice domain
00220 !    fraction input if appropriate
00221 !
00222 ! !REVISION HISTORY:
00223 !    2007-may-07 - M. Vertenstein - initial cpl7 version.
00224 !
00225 ! !INTERFACE: ------------------------------------------------------------------
00226 
00227 subroutine seq_frac_init(cdata_a, cdata_i, cdata_l, cdata_o, cdata_g, &
00228                          ice_present, ocn_present, lnd_present, glc_present, &
00229                          dead_comps, &
00230                          fractions_a, fractions_i, fractions_l, fractions_o, &
00231                          fractions_g)
00232 
00233 ! !INPUT/OUTPUT PARAMETERS:
00234 
00235    type(seq_cdata), intent(in)    :: cdata_a
00236    type(seq_cdata), intent(in)    :: cdata_i
00237    type(seq_cdata), intent(in)    :: cdata_l
00238    type(seq_cdata), intent(in)    :: cdata_o
00239    type(seq_cdata), intent(in)    :: cdata_g
00240    logical        , intent(in)    :: ice_present   ! .true. => ocean is present and not aqua_planet
00241    logical        , intent(in)    :: ocn_present   ! .true. => ocean is present
00242    logical        , intent(in)    :: lnd_present   ! .true. => land is present
00243    logical        , intent(in)    :: glc_present   ! .true. => glc is present
00244    logical        , intent(in)    :: dead_comps    ! .true. => dead models present
00245    type(mct_aVect), intent(inout) :: fractions_a   ! Fractions on atm grid/decomp
00246    type(mct_aVect), intent(inout) :: fractions_i   ! Fractions on ice grid/decomp
00247    type(mct_aVect), intent(inout) :: fractions_l   ! Fractions on lnd grid/decomp
00248    type(mct_aVect), intent(inout) :: fractions_o   ! Fractions on ocn grid/decomp
00249    type(mct_aVect), intent(inout) :: fractions_g   ! Fractions on glc grid/decomp
00250 
00251 !EOP
00252 
00253    !----- local -----
00254    integer :: j,n                      ! indices
00255    integer :: ka, ki, kl, ko, kf, kk   ! indices 
00256    integer :: lsize                    ! local size of ice av
00257    integer :: debug_old                ! old debug value
00258    type(mct_gGrid), pointer :: dom_a   ! pointer to component domain
00259    type(mct_gGrid), pointer :: dom_i   ! pointer to component domain
00260    type(mct_gGrid), pointer :: dom_l   ! pointer to component domain
00261    type(mct_gGrid), pointer :: dom_o   ! pointer to component domain
00262    type(mct_gGrid), pointer :: dom_g   ! pointer to component domain
00263 
00264    character(*),parameter :: fraclist_a = 'afrac:ifrac:ofrac:lfrac:lfrin'
00265    character(*),parameter :: fraclist_o = 'afrac:ifrac:ofrac:ifrad:ofrad'
00266    character(*),parameter :: fraclist_i = 'afrac:ifrac:ofrac'
00267    character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin'
00268    character(*),parameter :: fraclist_g = 'gfrac'
00269 
00270    !----- formats -----
00271    character(*),parameter :: subName = '(seq_frac_init) '
00272 
00273 !-------------------------------------------------------------------------------
00274 !
00275 !-------------------------------------------------------------------------------
00276 
00277    debug_old = seq_frac_debug
00278    seq_frac_debug = 2
00279 
00280    seq_frac_dead = dead_comps
00281 
00282    ! Initialize fractions on atm grid/decomp (initialize ice fraction to zero)
00283 
00284 !  if (atm_present) then
00285       call seq_cdata_setptrs(cdata_a, dom=dom_a)
00286 
00287       lSize = mct_aVect_lSize(dom_a%data)
00288       call mct_aVect_init(fractions_a,rList=fraclist_a,lsize=lsize)
00289       call mct_aVect_zero(fractions_a)
00290 
00291       ka = mct_aVect_indexRa(fractions_a,"afrac",perrWith=subName)
00292       fractions_a%rAttr(ka,:) = 1.0_r8
00293 !  endif
00294 
00295    ! Initialize fractions on glc grid decomp, just an initial "guess", updated later
00296 
00297    if (glc_present) then
00298       call seq_cdata_setptrs(cdata_g, dom=dom_g)
00299 
00300       lSize = mct_aVect_lSize(dom_g%data)
00301       call mct_aVect_init(fractions_g,rList=fraclist_g,lsize=lsize)
00302       call mct_aVect_zero(fractions_g)
00303       fractions_g%rAttr(:,:) = 1.0_r8
00304       
00305    end if
00306       
00307    ! Initialize fractions on land grid decomp, just an initial "guess", updated later
00308 
00309    if (lnd_present) then
00310       call seq_cdata_setptrs(cdata_l, dom=dom_l)
00311 
00312       lSize = mct_aVect_lSize(dom_l%data)
00313       call mct_aVect_init(fractions_l,rList=fraclist_l,lsize=lsize)
00314       call mct_aVect_zero(fractions_l)
00315       
00316       kk = mct_aVect_indexRA(fractions_l,"lfrin",perrWith=subName)
00317       kf = mct_aVect_indexRA(dom_l%data ,"frac" ,perrWith=subName)
00318       fractions_l%rAttr(kk,:) = dom_l%data%rAttr(kf,:)
00319 
00320       call map_lnd2atm_mct(cdata_l, fractions_l, cdata_a, fractions_a, fluxlist='lfrin')
00321       call map_atm2lnd_mct(cdata_a, fractions_a, cdata_l, fractions_l, fluxlist='afrac')
00322 
00323    end if
00324       
00325    ! Initialize fractions on ice grid/decomp (initialize ice fraction to zero)
00326 
00327    if (ice_present) then
00328       call seq_cdata_setptrs(cdata_i, dom=dom_i)
00329 
00330       lSize = mct_aVect_lSize(dom_i%data)
00331       call mct_aVect_init(fractions_i,rList=fraclist_i,lsize=lsize)
00332       call mct_aVect_zero(fractions_i)
00333       
00334       ko = mct_aVect_indexRa(fractions_i,"ofrac",perrWith=subName)      
00335       kf = mct_aVect_indexRA(dom_i%data ,"frac" ,perrWith=subName)
00336       fractions_i%rAttr(ko,:) = dom_i%data%rAttr(kf,:)
00337 
00338       call map_ice2atm_mct(cdata_i, fractions_i, cdata_a, fractions_a, fluxlist='ofrac')
00339 
00340    end if
00341 
00342    ! Initialize fractions on ocean grid/decomp (initialize ice fraction to zero)
00343    ! These are initialize the same as for ice
00344 
00345    if (ocn_present) then
00346       call seq_cdata_setptrs(cdata_o, dom=dom_o)
00347       
00348       lSize = mct_aVect_lSize(dom_o%data)
00349       call mct_aVect_init(fractions_o,rList=fraclist_o,lsize=lsize)
00350       call mct_aVect_zero(fractions_o)
00351       
00352       if (ice_present) then
00353          call map_ice2ocn_mct(cdata_i, fractions_i, cdata_o, fractions_o, fluxlist='ofrac')
00354       else
00355          ko = mct_aVect_indexRa(fractions_o,"ofrac",perrWith=subName)   
00356          kf = mct_aVect_indexRA(dom_o%data ,"frac" ,perrWith=subName)
00357          fractions_o%rAttr(ko,:) = dom_o%data%rAttr(kf,:)
00358          call map_ocn2atm_mct(cdata_o, fractions_o, cdata_a, fractions_a, fluxlist='ofrac')
00359       endif
00360 
00361       call map_atm2ocn_mct(cdata_a, fractions_a, cdata_o, fractions_o, fluxlist='afrac')
00362       if (ice_present) then
00363           ! --- this should be an atm2ice call above, but atm2ice doesn't work
00364           call map_ocn2ice_mct(cdata_o, fractions_o, cdata_i, fractions_i, fluxlist='afrac')
00365       endif
00366 
00367    end if
00368 
00369    ! --- Set ofrac and lfrac on atm grid.  These should actually be mapo2a of
00370    !     ofrac and lfrac but we can't map lfrac from o2a due to masked mapping
00371    !     weights.  So we have to settle for a residual calculation that is
00372    !     truncated to zero to try to preserve "all ocean" cells.
00373 
00374    ka = mct_aVect_indexRa(fractions_a,"afrac",perrWith=subName)
00375    ki = mct_aVect_indexRa(fractions_a,"ifrac",perrWith=subName) 
00376    kl = mct_aVect_indexRa(fractions_a,"lfrac",perrWith=subName)
00377    ko = mct_aVect_indexRa(fractions_a,"ofrac",perrWith=subName) 
00378    kk = mct_aVect_indexRa(fractions_a,"lfrin",perrWith=subName) 
00379    lSize = mct_aVect_lSize(fractions_a)
00380 
00381    if (ice_present .or. ocn_present) then
00382       do n = 1,lsize
00383          fractions_a%rAttr(kl,n) = 1.0_r8 - fractions_a%rAttr(ko,n)
00384          if (abs(fractions_a%rAttr(kl,n)) < eps_fraclim) then
00385             fractions_a%rAttr(kl,n) = 0.0_r8
00386             if (atm_frac_correct) fractions_a%rAttr(ko,n) = 1.0_r8
00387          endif
00388       enddo
00389    else if (lnd_present) then
00390       do n = 1,lsize
00391          fractions_a%rAttr(kl,n) = fractions_a%rAttr(kk,n)
00392          fractions_a%rAttr(ko,n) = 1.0_r8 - fractions_a%rAttr(kl,n)
00393          if (abs(fractions_a%rAttr(ko,n)) < eps_fraclim) then
00394             fractions_a%rAttr(ko,n) = 0.0_r8
00395             if (atm_frac_correct) fractions_a%rAttr(kl,n) = 1.0_r8
00396          endif
00397       enddo
00398    endif
00399 
00400    ! --- finally, set fractions_l(lfrac) from fractions_a(lfrac)
00401 
00402    if (lnd_present) then
00403       call map_atm2lnd_mct(cdata_a, fractions_a, cdata_l, fractions_l, fluxlist='lfrac')
00404    end if
00405 
00406    if (lnd_present) call seq_frac_check(fractions_l,'lnd init')
00407    if (glc_present) call seq_frac_check(fractions_g,'glc init')
00408    if (ice_present) call seq_frac_check(fractions_i,'ice init')
00409    if (ocn_present) call seq_frac_check(fractions_o,'ocn init')
00410    if (lnd_present.or.ice_present.or.ocn_present) &
00411                     call seq_frac_check(fractions_a,'atm init')
00412    seq_frac_debug = debug_old
00413 
00414 end subroutine seq_frac_init
00415 
00416 !===============================================================================
00417 !BOP ===========================================================================
00418 !
00419 ! !IROUTINE: seq_frac_set
00420 !
00421 ! !DESCRIPTION:
00422 !    Update surface fractions
00423 !
00424 ! !REVISION HISTORY:
00425 !    2007-may-07 - M. Vertenstein - initial cpl7 version.
00426 !
00427 ! !INTERFACE: ------------------------------------------------------------------
00428 
00429 subroutine seq_frac_set(i2x_i, cdata_a, cdata_i, cdata_l, cdata_o, cdata_g, &
00430                          ice_present, ocn_present, lnd_present, glc_present, &
00431                          fractions_a, fractions_i, fractions_l, fractions_o, &
00432                          fractions_g)
00433 
00434 ! !INPUT/OUTPUT PARAMETERS:
00435 
00436    type(mct_aVect), intent(in)    :: i2x_i
00437    type(seq_cdata), intent(in)    :: cdata_a
00438    type(seq_cdata), intent(in)    :: cdata_i
00439    type(seq_cdata), intent(in)    :: cdata_l
00440    type(seq_cdata), intent(in)    :: cdata_o
00441    type(seq_cdata), intent(in)    :: cdata_g
00442    logical        , intent(in)    :: ice_present   ! true => ice is present
00443    logical        , intent(in)    :: ocn_present   ! true => ocn is present
00444    logical        , intent(in)    :: lnd_present   ! true => lnd is present
00445    logical        , intent(in)    :: glc_present   ! true => glc is present
00446    type(mct_aVect), intent(inout) :: fractions_a   ! Fractions on atm 
00447    type(mct_aVect), intent(inout) :: fractions_i   ! Fractions on ice 
00448    type(mct_aVect), intent(inout) :: fractions_l   ! Fractions on lnd 
00449    type(mct_aVect), intent(inout) :: fractions_o   ! Fractions on ocn 
00450    type(mct_aVect), intent(inout) :: fractions_g   ! Fractions on glc
00451 
00452 !EOP
00453 
00454    !----- local -----
00455    integer         :: j, n
00456    integer         :: ki, kl, ka, ko, kf
00457    integer         :: lsize
00458    real(r8),allocatable :: fcorr(:)
00459    type(mct_gGrid), pointer :: dom_i   ! pointer to component domain
00460 
00461    !----- formats -----
00462    character(*),parameter :: subName = '(seq_frac_set) '
00463 
00464 !-------------------------------------------------------------------------------
00465 !
00466 !-------------------------------------------------------------------------------
00467 
00468    !----------------------------------------------------------------------
00469    ! Update fractions
00470    ! - Update ice fraction on ice grid first, normalize to total fraction 
00471    !   available for cover
00472    ! - Update ocn fraction on ice grid as residual
00473    ! - Map ice/ocn fractions from ice grid to ocean and atm grids
00474    !----------------------------------------------------------------------
00475 
00476    if (ice_present) then
00477       call mct_aVect_copy(i2x_i, fractions_i,"Si_ifrac","ifrac")
00478       call seq_cdata_setptrs(cdata_i, dom=dom_i)
00479       ki = mct_aVect_indexRA(fractions_i,"ifrac")
00480       ko = mct_aVect_indexRA(fractions_i,"ofrac")
00481       kf = mct_aVect_indexRA(dom_i%data ,"frac" ,perrWith=subName)
00482       fractions_i%rAttr(ki,:) = fractions_i%rAttr(ki,:) * dom_i%data%rAttr(kf,:)
00483       fractions_i%rAttr(ko,:) = dom_i%data%rAttr(kf,:) - fractions_i%rAttr(ki,:)
00484       call seq_frac_check(fractions_i,'ice set')
00485 
00486       if (ocn_present) then
00487          call map_ice2ocn_mct(cdata_i, fractions_i, cdata_o, fractions_o, fluxlist='ofrac:ifrac')
00488          call seq_frac_check(fractions_o,'ocn set')
00489       endif
00490 
00491 !      if (atm_present) then
00492          call map_ice2atm_mct(cdata_i, fractions_i, cdata_a, fractions_a, fluxlist='ofrac:ifrac')
00493 !tcx---  fraction correction, this forces the fractions_a to sum to 1.0_r8.  
00494 !   ---  but it introduces a conservation error in mapping
00495          if (atm_frac_correct) then
00496             ki = mct_aVect_indexRA(fractions_a,"ifrac")
00497             ko = mct_aVect_indexRA(fractions_a,"ofrac")
00498             kl = mct_aVect_indexRA(fractions_a,"lfrac")
00499             lSize = mct_aVect_lSize(fractions_a)
00500             allocate(fcorr(lsize))
00501             do n = 1,lsize
00502                if ((fractions_a%rAttr(ki,n)+fractions_a%rAttr(ko,n)) > 0.0_r8) then
00503                   fcorr(n) = ((1.0_r8-fractions_a%rAttr(kl,n))/(fractions_a%rAttr(ki,n)+fractions_a%rAttr(ko,n)))
00504                else
00505                   fcorr(n) = 0.0_r8
00506                endif
00507             enddo
00508             fractions_a%rAttr(ki,:) = fractions_a%rAttr(ki,:) * fcorr(:)
00509             fractions_a%rAttr(ko,:) = fractions_a%rAttr(ko,:) * fcorr(:)
00510             deallocate(fcorr)
00511          endif
00512 
00513          call seq_frac_check(fractions_a,'atm set')
00514 !     endif
00515    end if
00516 
00517 end subroutine seq_frac_set
00518 
00519 !===============================================================================
00520 !BOP ===========================================================================
00521 !
00522 ! !IROUTINE: seq_frac_check
00523 !
00524 ! !DESCRIPTION:
00525 !    Check surface fractions
00526 !
00527 ! !REVISION HISTORY:
00528 !    2008-jun-11 - T. Craig - initial version
00529 !
00530 ! !INTERFACE: ------------------------------------------------------------------
00531 
00532 subroutine seq_frac_check(fractions,string)
00533 
00534 ! !INPUT/OUTPUT PARAMETERS:
00535 
00536    type(mct_aVect) , intent(in) :: fractions   ! Fractions datatype
00537    character(len=*), intent(in), optional :: string      ! character string
00538 
00539 !EOP
00540 
00541    !----- local -----
00542    integer  :: n, lsize
00543    integer  :: ncnt
00544    real(r8) :: sum,diff,maxerr
00545    real(r8) :: aminval,amaxval   ! used for lnd
00546    real(r8) :: lminval,lmaxval   ! used for lnd
00547    real(r8) :: ominval,omaxval   ! used for ocn
00548    real(r8) :: iminval,imaxval   ! used for ice
00549    real(r8) :: gminval,gmaxval   ! used for glc
00550    real(r8) :: kminval,kmaxval   ! used for lnd
00551    real(r8) :: sminval,smaxval   ! used for sum
00552    integer  :: ka,kl,ki,ko,kg,kk
00553    character(len=128) :: lstring
00554    logical :: error
00555 
00556    !----- formats -----
00557    character(*),parameter :: subName = '(seq_frac_check) '
00558    character(*),parameter :: F01 = "('(seq_frac_check) ',2a,i8,g26.18)"
00559    character(*),parameter :: F02 = "('(seq_frac_check) ',2a,2g26.18)"
00560 
00561 !-------------------------------------------------------------------------------
00562 !
00563 !-------------------------------------------------------------------------------
00564 
00565    if (present(string)) then
00566       lstring='['//trim(string)//']'
00567    else
00568       lstring=''
00569    endif
00570 
00571    ka = -1
00572    kl = -1
00573    ki = -1
00574    ko = -1
00575    kk = -1
00576    kg = -1
00577    aminval = 0.0_r8
00578    amaxval = 0.0_r8
00579    lminval = 0.0_r8
00580    lmaxval = 0.0_r8
00581    iminval = 0.0_r8
00582    imaxval = 0.0_r8
00583    ominval = 0.0_r8
00584    omaxval = 0.0_r8
00585    gminval = 0.0_r8
00586    gmaxval = 0.0_r8
00587    kminval = 0.0_r8
00588    kmaxval = 0.0_r8
00589    sminval = 0.0_r8
00590    smaxval = 0.0_r8
00591 
00592    lsize = mct_avect_lsize(fractions)
00593    ka = mct_aVect_indexRA(fractions,"afrac",perrWith='quiet')
00594    kl = mct_aVect_indexRA(fractions,"lfrac",perrWith='quiet')
00595    ki = mct_aVect_indexRA(fractions,"ifrac",perrWith='quiet')
00596    ko = mct_aVect_indexRA(fractions,"ofrac",perrWith='quiet')
00597    kg = mct_aVect_indexRA(fractions,"gfrac",perrWith='quiet')
00598    kk = mct_aVect_indexRA(fractions,"lfrin",perrWith='quiet')
00599 
00600    if (ka > 0) then
00601       aminval = minval(fractions%rAttr(ka,:))
00602       amaxval = maxval(fractions%rAttr(ka,:))
00603    endif
00604    if (kl > 0) then
00605       lminval = minval(fractions%rAttr(kl,:))
00606       lmaxval = maxval(fractions%rAttr(kl,:))
00607    endif
00608    if (ko > 0) then
00609       ominval = minval(fractions%rAttr(ko,:))
00610       omaxval = maxval(fractions%rAttr(ko,:))
00611    endif
00612    if (ki > 0) then
00613       iminval = minval(fractions%rAttr(ki,:))
00614       imaxval = maxval(fractions%rAttr(ki,:))
00615    endif
00616    if (kg > 0) then
00617       gminval = minval(fractions%rAttr(kg,:))
00618       gmaxval = maxval(fractions%rAttr(kg,:))
00619    endif
00620    if (kk > 0) then
00621       kminval = minval(fractions%rAttr(kk,:))
00622       kmaxval = maxval(fractions%rAttr(kk,:))
00623    endif
00624 
00625    ncnt = 0
00626    maxerr = 0.0_r8
00627    if (kl > 0 .and. ko > 0 .and. ki > 0) then
00628       sminval =  9999._r8
00629       smaxval = -9999._r8
00630       do n = 1,lsize
00631          sum = fractions%rAttr(ko,n) + fractions%rAttr(kl,n) + fractions%rAttr(ki,n)
00632          sminval = min(sum,sminval)
00633          smaxval = max(sum,smaxval)
00634          diff = abs(1.0_r8 - sum)
00635          if (diff > eps_fracsum) then
00636             ncnt = ncnt + 1
00637             maxerr = max(maxerr, diff)
00638 !tcx debug            write(logunit,*) trim(lstring),' err# ',ncnt,n,lsize,fractions%rAttr(ko,n),fractions%rAttr(kl,n),fractions%rAttr(ki,n),sum
00639          endif
00640       enddo
00641    endif
00642 
00643    error = .false.
00644    if (ncnt > 0) error = .true.
00645    if (aminval < 0.0_r8-eps_fracval .or. amaxval > 1.0_r8+eps_fracval) error = .true.
00646    if (lminval < 0.0_r8-eps_fracval .or. lmaxval > 1.0_r8+eps_fracval) error = .true.
00647    if (ominval < 0.0_r8-eps_fracval .or. omaxval > 1.0_r8+eps_fracval) error = .true.
00648    if (iminval < 0.0_r8-eps_fracval .or. imaxval > 1.0_r8+eps_fracval) error = .true.
00649    if (gminval < 0.0_r8-eps_fracval .or. gmaxval > 1.0_r8+eps_fracval) error = .true.
00650    if (kminval < 0.0_r8-eps_fracval .or. kmaxval > 1.0_r8+eps_fracval) error = .true.
00651 
00652    if (error .or. seq_frac_debug > 1) then
00653       if (ka > 0) &
00654          write(logunit,F02) trim(lstring),' afrac min/max   = ',aminval,amaxval
00655       if (kl > 0) &
00656          write(logunit,F02) trim(lstring),' lfrac min/max   = ',lminval,lmaxval
00657       if (kg > 0) &
00658          write(logunit,F02) trim(lstring),' gfrac min/max   = ',gminval,gmaxval
00659       if (ko > 0) &
00660          write(logunit,F02) trim(lstring),' ofrac min/max   = ',ominval,omaxval
00661       if (ki > 0) &
00662          write(logunit,F02) trim(lstring),' ifrac min/max   = ',iminval,imaxval
00663       if (kk > 0) &
00664          write(logunit,F02) trim(lstring),' lfrin min/max   = ',kminval,kmaxval
00665       if (kl > 0 .and. ko > 0 .and. ki > 0) then
00666          write(logunit,F02) trim(lstring),' sum min/max     = ',sminval,smaxval
00667          write(logunit,F01) trim(lstring),' sum ncnt/maxerr = ',ncnt,maxerr
00668       endif
00669       if (error .and. .not. seq_frac_dead .and. seq_frac_abort) then
00670          write(logunit,F02) trim(lstring),' ERROR aborting '
00671          call shr_sys_abort()
00672       elseif (error) then
00673          write(logunit,F02) trim(lstring),' ERROR but NOT aborting '
00674       endif
00675    endif
00676 
00677 end subroutine seq_frac_check
00678 
00679 end module seq_frac_mct