CESM_Driver: seq_domain_mct.F90 Source File

00001 module seq_domain_mct
00002 
00003   use shr_kind_mod, only: R8=>shr_kind_r8, IN=>shr_kind_in
00004   use shr_kind_mod, only: CL=>shr_kind_cl
00005   use shr_sys_mod,  only: shr_sys_flush, shr_sys_abort
00006   use shr_mpi_mod,  only: shr_mpi_min, shr_mpi_max
00007 
00008   use mct_mod
00009   use seq_cdata_mod
00010   use seq_comm_mct
00011   use seq_infodata_mod
00012 
00013   use map_atmlnd_mct
00014   use map_atmocn_mct
00015   use map_atmice_mct
00016   use map_iceocn_mct
00017   use map_snoglc_mct
00018 
00019   implicit none
00020   save
00021   private ! except
00022 
00023 !--------------------------------------------------------------------------
00024 ! Public interfaces
00025 !--------------------------------------------------------------------------
00026 
00027   public :: seq_domain_check_mct
00028   public :: domain_areafactinit_mct
00029 
00030 !--------------------------------------------------------------------------
00031 ! Public variables
00032 !--------------------------------------------------------------------------
00033 
00034   real(R8), parameter :: eps_tiny   = 1.0e-16_R8 ! roundoff eps
00035   real(R8), parameter :: eps_big    = 1.0e+02_R8 ! big eps
00036   real(R8), parameter :: eps_frac_samegrid = 1.0e-14_R8 ! epsilon for fractions for samegrid
00037 
00038 !--------------------------------------------------------------------------
00039 ! Private interfaces
00040 !--------------------------------------------------------------------------
00041 
00042 #ifdef CPP_VECTOR
00043   logical :: usevector = .true.
00044 #else
00045   logical :: usevector = .false.
00046 #endif
00047   
00048 #ifdef SYSUNICOS
00049   logical :: usealltoall = .true.
00050 #else
00051   logical :: usealltoall = .false.
00052 #endif
00053   
00054   private :: seq_domain_check_grid_mct
00055 
00056 #include <mpif.h>  
00057 
00058 !================================================================================
00059 contains
00060 !================================================================================
00061 
00062 !================================================================================
00063 
00064   subroutine seq_domain_check_mct( cdata_a, cdata_i, cdata_l, cdata_o, cdata_r, cdata_g, cdata_s)
00065 
00066     !-----------------------------------------------------------
00067     !
00068     ! Arguments
00069     !
00070     type(seq_cdata), intent(in) :: cdata_a
00071     type(seq_cdata), intent(in) :: cdata_i
00072     type(seq_cdata), intent(in) :: cdata_l
00073     type(seq_cdata), intent(in) :: cdata_o
00074     type(seq_cdata), intent(in) :: cdata_r
00075     type(seq_cdata), intent(in) :: cdata_g
00076     type(seq_cdata), intent(in) :: cdata_s
00077     !
00078     ! Local variables
00079     !
00080     type(mct_gGrid)  , pointer :: atmdom_a   ! atm domain
00081     type(mct_gGrid)  , pointer :: icedom_i   ! ice domain
00082     type(mct_gGrid)  , pointer :: lnddom_l   ! lnd domain
00083     type(mct_gGrid)  , pointer :: ocndom_o   ! ocn domain
00084     type(mct_gGrid)  , pointer :: glcdom_g   ! glc domain
00085     type(mct_gGrid)  , pointer :: snodom_s   ! sno domain
00086     !
00087     type(mct_gsMap)  , pointer :: gsMap_a    ! atm global seg map 
00088     type(mct_gsMap)  , pointer :: gsMap_i    ! ice global seg map 
00089     type(mct_gsMap)  , pointer :: gsMap_l    ! lnd global seg map 
00090     type(mct_gsMap)  , pointer :: gsMap_o    ! ocn global seg map 
00091     type(mct_gsMap)  , pointer :: gsMap_r    ! ocn global seg map 
00092     type(mct_gsMap)  , pointer :: gsMap_g    ! glc global seg map 
00093     type(mct_gsMap)  , pointer :: gsMap_s    ! sno global seg map 
00094     !
00095     type(mct_gGrid) :: lnddom_a              ! lnd domain info on atm decomp
00096     type(mct_gGrid) :: icedom_a              ! ice domain info on atm decomp (all grids same)
00097     type(mct_gGrid) :: ocndom_a              ! ocn domain info on atm decomp (all grids same)
00098     type(mct_gGrid) :: icedom_o              ! ocn domain info on ocn decomp (atm/ocn grid different)
00099     type(mct_gGrid) :: snodom_g              ! sno domain info on glc decomp
00100     !
00101     integer(IN) :: mpicom_a                  ! atm mpicom
00102     integer(IN) :: mpicom_i                  ! ice mpicom
00103     integer(IN) :: mpicom_l                  ! lnd mpicom
00104     integer(IN) :: mpicom_o                  ! ocn mpicom
00105     integer(IN) :: mpicom_g                  ! glc mpicom
00106     integer(IN) :: mpicom_s                  ! sno mpicom
00107     !
00108     type(seq_infodata_type), pointer :: infodata
00109     !
00110     real(R8), pointer :: fracl(:)            ! land fraction on atm decomp 
00111     real(R8), pointer :: fraco(:)            ! ocn  fraction on atm decomp 
00112     real(R8), pointer :: fraci(:)            ! ice  fraction on atm decomp 
00113     real(R8), pointer :: maskl(:)            ! land mask on atm decomp (all grids same)
00114     real(R8), pointer :: maski(:)            ! ice  mask on atm decomp (all grids same)
00115     real(R8), pointer :: masko(:)            ! ocn  mask on atm decomp (all grids same)
00116     !
00117     integer(IN) :: n, kl, ko, ki             ! indicies
00118     integer(IN) :: k1,k2,k3                  ! indicies
00119     !
00120     logical      :: lnd_present              ! lnd present flag
00121     logical      :: ocn_present              ! ocn present flag
00122     logical      :: ice_present              ! ice present flag
00123     logical      :: glc_present              ! glc present flag
00124     logical      :: sno_present              ! sno present flag
00125     logical      :: rof_present              ! rof present flag
00126     logical      :: ocnrof_prognostic        ! ocn rof prognostic flag
00127     logical      :: samegrid_ao              ! atm ocn grid same
00128     logical      :: samegrid_ro              ! rof ocn grid same
00129     logical      :: samegrid_al              ! atm lnd grid same
00130     integer(IN)  :: rcode                    ! error status
00131     integer(IN)  :: atmsize                  ! local  size of atm  grid
00132     integer(IN)  :: lndsize                  ! local  size of land grid
00133     integer(IN)  :: ocnsize                  ! local  size of ocn  grid
00134     integer(IN)  :: icesize                  ! local  size of ice  grid
00135     integer(IN)  :: glcsize                  ! local  size of glc  grid
00136     integer(IN)  :: snosize                  ! local  size of sno  grid
00137     integer(IN)  :: gatmsize                 ! global size of atm  grid
00138     integer(IN)  :: glndsize                 ! global size of land grid
00139     integer(IN)  :: gocnsize                 ! global size of ocn  grid
00140     integer(IN)  :: grofsize                 ! global size of ocn  grid
00141     integer(IN)  :: gicesize                 ! global size of ice  grid
00142     integer(IN)  :: gglcsize                 ! global size of glc  grid
00143     integer(IN)  :: gsnosize                 ! global size of sno  grid
00144     integer(IN)  :: npts                     ! local size temporary
00145     integer(IN)  :: ier                      ! error code
00146     real(R8)     :: diff,dmaxo,dmaxi         ! difference tracker
00147     logical      :: iamroot                  ! local masterproc
00148     real(R8)     :: eps_frac                 ! epsilon for fractions
00149     real(R8)     :: eps_axmask               ! epsilon for masks, atm/lnd
00150     real(R8)     :: eps_axgrid               ! epsilon for grid coords, atm/lnd
00151     real(R8)     :: eps_axarea               ! epsilon for areas, atm/lnd
00152     real(R8)     :: eps_oimask               ! epsilon for masks, ocn/ice
00153     real(R8)     :: eps_oigrid               ! epsilon for grid coords, ocn/ice
00154     real(R8)     :: eps_oiarea               ! epsilon for areas, ocn/ice
00155     real(R8)     :: my_eps_frac              ! local eps_frac value
00156     real(r8)     :: rmin1,rmax1,rmin,rmax    ! local min max computation
00157     !
00158     real(R8),allocatable :: mask (:)         ! temporary real vector, domain mask
00159     !
00160     character(*),parameter :: F00 = "('(domain_check_mct) ',4a)"
00161     character(*),parameter :: F01 = "('(domain_check_mct) ',a,i6,a)"
00162     character(*),parameter :: F02 = "('(domain_check_mct) ',a,g23.15)"
00163     character(*),parameter :: F0R = "('(domain_check_mct) ',2A,2g23.15,A )"
00164     character(*),parameter :: subName = '(domain_check_mct) '
00165     !-----------------------------------------------------------
00166 
00167     call seq_comm_setptrs(CPLID,iamroot=iamroot)
00168 
00169     call seq_cdata_setptrs(cdata_a, infodata=infodata)
00170     call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al, &
00171          samegrid_ro=samegrid_ro)
00172     call seq_infodata_GetData( infodata, lnd_present=lnd_present, ocn_present=ocn_present, &
00173          ice_present=ice_present, glc_present=glc_present, sno_present=sno_present, &
00174          rof_present=rof_present, ocnrof_prognostic=ocnrof_prognostic)
00175     call seq_infodata_GetData( infodata, eps_frac=eps_frac, &
00176          eps_amask=eps_axmask, eps_agrid=eps_axgrid, eps_aarea=eps_axarea, &
00177          eps_omask=eps_oimask, eps_ogrid=eps_oigrid, eps_oarea=eps_oiarea )
00178 
00179     ! Get info
00180 
00181     call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=atmdom_a, mpicom=mpicom_a)
00182     atmsize = mct_avect_lsize(atmdom_a%data)
00183     gatmsize = mct_gsMap_gsize(gsMap_a)
00184 
00185     if (lnd_present) then
00186        call seq_cdata_setptrs(cdata_l, gsMap=gsMap_l, dom=lnddom_l, mpicom=mpicom_l)
00187        lndsize = mct_avect_lsize(lnddom_l%data)
00188        glndsize = mct_gsMap_gsize(gsMap_l) 
00189        if (samegrid_al .and. gatmsize /= glndsize) then
00190           write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global lndsize= ',glndsize
00191           call shr_sys_flush(logunit)
00192           call mct_die(subname,' atm and lnd grid must have the same global size')
00193        end if
00194        call mct_gGrid_init(oGGrid=lnddom_a, iGGrid=lnddom_l, lsize=atmsize)
00195        call mct_aVect_zero(lnddom_a%data)
00196        call map_lnd2atm_mct(cdata_l, lnddom_l%data, cdata_a, lnddom_a%data)
00197        allocate(maskl(atmsize),stat=rcode)
00198        if(rcode /= 0) call mct_die(subName,'allocate maskl')
00199        allocate(fracl(atmsize),stat=rcode)
00200        if(rcode /= 0) call mct_die(subName,'allocate fracl')
00201        call mct_aVect_exportRAttr(lnddom_a%data, 'mask', maskl, atmsize)
00202        call mct_aVect_exportRAttr(lnddom_a%data, 'frac', fracl, atmsize)
00203     endif
00204 
00205     if (ocn_present) then
00206        call seq_cdata_setptrs(cdata_o, gsMap=gsMap_o, dom=ocndom_o, mpicom=mpicom_o)
00207        ocnsize = mct_avect_lsize(ocndom_o%data)
00208        gocnsize = mct_gsMap_gsize(gsMap_o)
00209        if (samegrid_ao .and. gatmsize /= gocnsize) then
00210           write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global ocnsize= ',gocnsize
00211           call shr_sys_flush(logunit)
00212           call mct_die(subname,' atm and ocn grid must have the same global size')
00213        end if
00214        call mct_gGrid_init(oGGrid=ocndom_a, iGGrid=ocndom_o, lsize=atmsize)
00215        call mct_aVect_zero(ocndom_a%data)
00216        call map_ocn2atm_mct(cdata_o, ocndom_o%data, cdata_a, ocndom_a%data)
00217        allocate(masko(atmsize),stat=rcode)
00218        if(rcode /= 0) call mct_die(subName,'allocate masko')
00219        allocate(fraco(atmsize),stat=rcode)
00220        if(rcode /= 0) call mct_die(subName,'allocate fraco')
00221        call mct_aVect_exportRAttr(ocndom_a%data, 'mask', masko, atmsize)
00222        if (samegrid_ao) then
00223           call mct_aVect_exportRattr(ocndom_a%data, 'frac', fraco, atmsize)
00224        else
00225           call mct_aVect_exportRattr(ocndom_a%data, 'mask', fraco, atmsize)
00226        endif
00227     endif
00228    
00229     if (ice_present) then
00230        call seq_cdata_setptrs(cdata_i, gsMap=gsMap_i, dom=icedom_i, mpicom=mpicom_i)
00231        icesize = mct_avect_lsize(icedom_i%data)
00232        gicesize = mct_gsMap_gsize(gsMap_i) 
00233        if (samegrid_ao .and. gatmsize /= gicesize) then
00234           write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global icesize= ',gicesize
00235           call shr_sys_flush(logunit)
00236           call mct_die(subname,' atm and ice grid must have the same global size')
00237        end if
00238        call mct_gGrid_init(oGGrid=icedom_a, iGGrid=icedom_i, lsize=atmsize)
00239        call mct_aVect_zero(icedom_a%data)
00240        call map_ice2atm_mct(cdata_i, icedom_i%data, cdata_a, icedom_a%data)
00241        allocate(maski(atmsize),stat=rcode)
00242        if(rcode /= 0) call mct_die(subName,'allocate maski')
00243        allocate(fraci(atmsize),stat=rcode)
00244        if(rcode /= 0) call mct_die(subName,'allocate fraci')
00245        call mct_aVect_exportRAttr(icedom_a%data, 'mask', maski, atmsize)
00246        if (samegrid_ao) then
00247           call mct_aVect_exportRattr(icedom_a%data, 'frac', fraci, atmsize)
00248        else
00249           call mct_aVect_exportRattr(icedom_a%data, 'mask', fraci, atmsize)
00250        endif
00251     endif
00252 
00253     if (sno_present .and. glc_present) then
00254        call seq_cdata_setptrs(cdata_g, gsMap=gsMap_g, dom=glcdom_g, mpicom=mpicom_g)
00255        glcsize = mct_avect_lsize(glcdom_g%data)
00256        gglcsize = mct_gsMap_gsize(gsMap_g) 
00257        call seq_cdata_setptrs(cdata_s, gsMap=gsMap_s, dom=snodom_s, mpicom=mpicom_s)
00258        snosize = mct_avect_lsize(snodom_s%data)
00259        gsnosize = mct_gsMap_gsize(gsMap_s) 
00260        if (gglcsize /= gsnosize) then
00261           write(logunit,*) subname,' error: global glcsize = ',gglcsize,' global snosize= ',gsnosize
00262           call shr_sys_flush(logunit)
00263           call mct_die(subname,' glc and sno grid must have the same global size')
00264        end if
00265        call mct_gGrid_init(oGGrid=snodom_g, iGGrid=snodom_s, lsize=glcsize)
00266        call mct_aVect_zero(snodom_g%data)
00267        call map_sno2glc_mct(cdata_s, snodom_s%data, cdata_g, snodom_g%data)
00268        if (iamroot) write(logunit,F00) ' --- checking glc/sno domains ---'
00269        npts = glcsize
00270        allocate(mask(npts),stat=rcode)
00271        if(rcode /= 0) call mct_die(subName,'allocate mask')
00272        call mct_aVect_getRAttr(snodom_g%data,"mask",mask,rcode)
00273        where (mask < eps_axmask) mask = 0.0_r8
00274        call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'mask', eps=eps_axmask, mpicom=mpicom_g, mask=mask)
00275        call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
00276        call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
00277        call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'area', eps=eps_axarea, mpicom=mpicom_g, mask=mask)
00278        deallocate(mask,stat=rcode)
00279        if(rcode /= 0) call mct_die(subName,'deallocate mask')
00280     endif
00281 
00282     if (ice_present .and. ocn_present) then
00283        if (gocnsize /= gicesize) then
00284           write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global icesize= ',gicesize
00285           call shr_sys_flush(logunit)
00286           call mct_die(subname,' ocean and ice grid must have the same global size')
00287        endif
00288        call mct_gGrid_init(oGGrid=icedom_o, iGGrid=icedom_i, lsize=ocnsize)
00289        call mct_aVect_zero(icedom_o%data)
00290        call map_ice2ocn_mct(cdata_i, icedom_i%data, cdata_o, icedom_o%data)
00291     end if
00292 
00293     if (rof_present .and. ocnrof_prognostic .and. samegrid_ro) then
00294        call seq_cdata_setptrs(cdata_r, gsMap=gsMap_r)
00295        grofsize = mct_gsMap_gsize(gsMap_r)
00296        if (gocnsize /= grofsize) then
00297           write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global rofsize= ',grofsize
00298           call shr_sys_flush(logunit)
00299           call mct_die(subname,' ocean and rof grid must have the same global size')
00300        endif
00301     end if
00302 
00303     !------------------------------------------------------------------------------
00304     ! Check ice/ocean grid consistency
00305     !------------------------------------------------------------------------------
00306 
00307      if (ocn_present .and. ice_present) then
00308 !    if (samegrid_oi) then       ! doesn't yet exist
00309 
00310        npts = ocnsize
00311        allocate(mask(npts),stat=rcode)
00312        if(rcode /= 0) call mct_die(subName,'allocate mask')
00313 
00314        if (iamroot) write(logunit,F00) ' --- checking ocn/ice domains ---'
00315        call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'mask', eps=eps_oigrid, mpicom=mpicom_o)
00316        call mct_aVect_getRAttr(ocndom_o%data,"mask",mask,rcode)
00317        where (mask < eps_oimask) mask = 0.0_r8
00318 
00319        call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lat' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
00320        call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lon' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
00321        call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'area', eps=eps_oiarea, mpicom=mpicom_o, mask=mask)
00322 
00323        deallocate(mask,stat=rcode)
00324        if(rcode /= 0) call mct_die(subName,'deallocate mask')
00325 
00326 !    endif
00327      endif
00328 
00329     !------------------------------------------------------------------------------
00330     ! Check atm/lnd grid consistency
00331     !------------------------------------------------------------------------------
00332 
00333     if (lnd_present .and. samegrid_al) then
00334        if (iamroot) write(logunit,F00) ' --- checking atm/land domains ---'
00335        call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
00336        call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
00337        call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maskl)
00338     endif
00339 
00340     !------------------------------------------------------------------------------
00341     ! Check atm/ocn and atm/ice grid consistency (if samegrid)
00342     !------------------------------------------------------------------------------
00343 
00344     if (ice_present .and. samegrid_ao) then
00345        if (iamroot) write(logunit,F00) ' --- checking atm/ice domains ---'
00346        call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
00347        call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
00348        call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maski)
00349     endif
00350 
00351     if (ocn_present .and. samegrid_ao) then
00352        if (iamroot) write(logunit,F00) ' --- checking atm/ocn domains ---'
00353        call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
00354        call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
00355        call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=masko)
00356     endif
00357 
00358     !------------------------------------------------------------------------------
00359     ! Check consistency of land fraction with ocean mask on grid
00360     !------------------------------------------------------------------------------
00361 
00362     my_eps_frac = eps_frac
00363     if (samegrid_ao) my_eps_frac = eps_frac_samegrid
00364     if (.not. samegrid_al) my_eps_frac = eps_big
00365 
00366     if (iamroot) write(logunit,F00) ' --- checking fractions in domains ---'
00367     dmaxi = 0.0_R8
00368     dmaxo = 0.0_R8
00369     do n = 1,atmsize
00370        if (lnd_present .and. ice_present) then
00371           diff = abs(1._r8 - fracl(n) - fraci(n))
00372           dmaxi = max(diff,dmaxi)
00373           if (diff > my_eps_frac) then
00374              write(logunit,*)'inconsistency between land fraction and ice land fraction'
00375              write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n),' sum= ',fracl(n)+fraci(n)
00376              call shr_sys_flush(logunit)
00377              call mct_die(subname,'inconsistency between land fraction and ice land fraction')
00378           end if
00379           if ((1._r8-fraci(n)) > eps_frac .and. fracl(n) < eps_tiny) then
00380              write(logunit,*)'inconsistency between land mask and ice land mask'
00381              write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n)
00382              call shr_sys_flush(logunit)
00383              call mct_die(subname,' inconsistency between land mask and ice land mask')
00384           end if
00385        endif
00386        if (lnd_present .and. ocn_present) then
00387           diff = abs(1._r8 - fracl(n) - fraco(n))
00388           dmaxo = max(diff,dmaxo)
00389           if (diff > my_eps_frac) then
00390              write(logunit,*)'inconsistency between land fraction and ocn land fraction'
00391              write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n),' sum= ',fracl(n)+fraco(n)
00392              call shr_sys_flush(logunit)
00393              call mct_die(subname,' inconsistency between land fraction and ocn land fraction')
00394           end if
00395           if ((1._r8-fraco(n)) > eps_frac .and. fracl(n) < eps_tiny) then
00396              write(logunit,*)'inconsistency between land mask and ocn land mask'
00397              write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n)
00398              call shr_sys_flush(logunit)
00399              call mct_die(subname,' inconsistency between land mask and ocn land mask')
00400           end if
00401        endif
00402     end do 
00403     if (iamroot) then
00404        write(logunit,F02) ' maximum           difference for ofrac sum ',dmaxo
00405        write(logunit,F02) ' maximum           difference for ifrac sum ',dmaxi
00406        write(logunit,F02) ' maximum allowable difference for  frac sum ',my_eps_frac
00407        write(logunit,F02) ' maximum allowable tolerance for valid frac ',eps_frac
00408        call shr_sys_flush(logunit)
00409     endif
00410 
00411     !------------------------------------------------------------------------------
00412     ! Set atm and lnd ascale
00413     !------------------------------------------------------------------------------
00414     if (lnd_present .and. ocn_present) then
00415        k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
00416        do k2 = 1,atmsize
00417           if (fracl(k2) /= 0.0_r8) then
00418              atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraco(k2))/(fracl(k2))
00419           else
00420              atmdom_a%data%rAttr(k1,k2) = 0.0
00421           endif
00422        enddo
00423        call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
00424     elseif (lnd_present .and. ice_present) then
00425        k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
00426        do k2 = 1,atmsize
00427           if (fracl(k2) /= 0.0_r8) then
00428              atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraci(k2))/(fracl(k2))
00429           else
00430              atmdom_a%data%rAttr(k1,k2) = 0.0
00431           endif
00432        enddo
00433        call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
00434     endif
00435 
00436     if (lnd_present .and. (ocn_present .or. ice_present)) then
00437        k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check atm ascale')
00438        rmin1 = minval(atmdom_a%data%rAttr(k1,:))
00439        rmax1 = maxval(atmdom_a%data%rAttr(k1,:))
00440        call shr_mpi_min(rmin1,rmin,mpicom_a)
00441        call shr_mpi_max(rmax1,rmax,mpicom_a)
00442        if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' atmdom_a'
00443 
00444        k1 = mct_aVect_indexRa(lnddom_l%data,"ascale",perrWith='domain_check lnd ascale')
00445        rmin1 = minval(lnddom_l%data%rAttr(k1,:))
00446        rmax1 = maxval(lnddom_l%data%rAttr(k1,:))
00447        call shr_mpi_min(rmin1,rmin,mpicom_l)
00448        call shr_mpi_max(rmax1,rmax,mpicom_l)
00449        if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' lnddom_l'
00450     endif
00451 
00452     !------------------------------------------------------------------------------
00453     ! Clean up allocated memory
00454     !------------------------------------------------------------------------------
00455 
00456     if (lnd_present) then
00457        deallocate(fracl,stat=rcode)
00458        if(rcode /= 0) call mct_die(subName,'deallocate fracl')
00459        deallocate(maskl,stat=rcode)
00460        if(rcode /= 0) call mct_die(subName,'deallocate maskl')
00461        call mct_gGrid_clean(lnddom_a, rcode)
00462        if(rcode /= 0) call mct_die(subName,'clean lnddom_a')
00463     endif
00464 
00465     if (ocn_present) then
00466        deallocate(fraco,stat=rcode)
00467        if(rcode /= 0) call mct_die(subName,'deallocate fraco')
00468        deallocate(masko,stat=rcode)
00469        if(rcode /= 0) call mct_die(subName,'deallocate masko')
00470        call mct_gGrid_clean(ocndom_a, rcode)
00471        if(rcode /= 0) call mct_die(subName,'clean ocndom_a')
00472     endif
00473 
00474     if (ice_present) then
00475        deallocate(fraci,stat=rcode)
00476        if(rcode /= 0) call mct_die(subName,'deallocate fraci')
00477        deallocate(maski,stat=rcode)
00478        if(rcode /= 0) call mct_die(subName,'deallocate maski')
00479        call mct_gGrid_clean(icedom_a, rcode)
00480        if(rcode /= 0) call mct_die(subName,'clean icedom_o')
00481     endif
00482 
00483     if (ocn_present .and. ice_present) then
00484        call mct_gGrid_clean(icedom_o, rcode)
00485        if(rcode /= 0) call mct_die(subName,'clean icedom_o')
00486     endif
00487 
00488   end subroutine seq_domain_check_mct
00489 
00490 !===============================================================================
00491   
00492   subroutine seq_domain_check_grid_mct(dom1, dom2, attr, eps, mpicom, mask)
00493    
00494     !-----------------------------------------------------------
00495 
00496     ! Arguments
00497 
00498     type(mct_aVect) , intent(in) :: dom1
00499     type(mct_aVect) , intent(in) :: dom2
00500     character(len=*), intent(in) :: attr   ! grid attribute to compare
00501     real(R8)        , intent(in) :: eps    ! error condition for compare
00502     integer(IN)     , intent(in) :: mpicom
00503     real(R8)        , intent(in), optional :: mask(:)
00504 
00505     ! Local variables
00506 
00507     integer(in)       :: n,ndiff            ! indices
00508     integer(in)       :: npts1,npts2,npts   ! counters
00509     integer(in)       :: rcode              ! error code
00510     real(R8)          :: diff,max_diff      ! temporaries
00511     real(R8)          :: tot_diff           ! maximum diff across all pes
00512     integer(IN)       :: ier                ! error code
00513     real(R8), pointer :: data1(:)           ! temporaries
00514     real(R8), pointer :: data2(:)           ! temporaries
00515     real(R8), pointer :: lmask(:)           ! temporaries
00516     logical           :: iamroot            ! local masterproc
00517 
00518     character(*),parameter :: F00 = "('(domain_check_grid_mct) ',4a)"
00519     character(*),parameter :: F01 = "('(domain_check_grid_mct) ',a,i12,a)"
00520     character(*),parameter :: F02 = "('(domain_check_grid_mct) ',2a,g23.15)"
00521     character(*),parameter :: F0R = "('(domain_check_grid_mct) ',2A,2g23.15,A )"
00522     character(*),parameter :: subName = '(domain_check_grid_mct) '
00523     !-----------------------------------------------------------
00524 
00525     call seq_comm_setptrs(CPLID,iamroot=iamroot)
00526 
00527    npts1 = mct_aVect_lsize(dom1)
00528    npts2 = mct_aVect_lsize(dom2)
00529    npts  = npts1
00530    
00531    if (npts1 == npts2) then
00532       if (iamroot) write(logunit,F01) " the domain size is = ", npts
00533    else
00534       write(logunit,*) trim(subname)," domain size #1 = ", npts1
00535       write(logunit,*) trim(subname)," domain size #2 = ", npts2
00536       write(logunit,*) trim(subname)," ERROR: domain size mis-match"
00537       call mct_die(subName,"ERROR: domain size mis-match")
00538    end if
00539    
00540    allocate(data1(npts),stat=rcode)
00541    if(rcode /= 0) call mct_die(subName,'allocate data1')
00542    allocate(data2(npts),stat=rcode)
00543    if(rcode /= 0) call mct_die(subName,'allocate data2')
00544    allocate(lmask(npts),stat=rcode)
00545    if(rcode /= 0) call mct_die(subName,'allocate lmask')
00546 
00547    call mct_aVect_exportRAttr(dom1, trim(attr), data1, npts)
00548    call mct_aVect_exportRAttr(dom2, trim(attr), data2, npts)
00549    lmask = 1.0_r8
00550    if (present(mask)) then
00551       if (size(mask) /= npts) then
00552          call mct_die(subName,"ERROR: mask size mis-match")
00553       endif
00554       lmask = mask
00555    endif
00556 
00557    ! --- adjust lons to address wraparound issues, we're assuming degree here! ---
00558 
00559    if (trim(attr) == "lon") then
00560       do n = 1,npts
00561          if (data2(n) > data1(n)) then
00562             do while ( (data1(n)+360.0_R8) < (data2(n)+180.0_R8) ) ! longitude is periodic
00563                data1(n) = data1(n) + 360.0_R8
00564             end do
00565          else
00566             do while ( (data2(n)+360.0_R8) < (data1(n)+180.0_R8) ) ! longitude is periodic
00567                data2(n) = data2(n) + 360.0_R8
00568             end do
00569          endif
00570       enddo
00571    endif
00572 
00573    ! Only check consistency where mask is greater than zero, if mask is present
00574 
00575    max_diff = 0.0_R8
00576    ndiff = 0
00577    do n=1,npts
00578       if (lmask(n) > eps_tiny) then
00579          diff = abs(data1(n)-data2(n))
00580          max_diff = max(max_diff,diff)
00581          if (diff > eps) then
00582 !debug            write(logunit,*)'n= ',n,' data1= ',data1(n),' data2= ',data2(n),' diff= ',diff, ' eps= ',eps
00583             ndiff = ndiff + 1
00584          endif
00585       end if
00586    end do
00587 
00588    call mpi_reduce(max_diff,tot_diff,1,MPI_REAL8,MPI_MAX,0,mpicom,ier)
00589    if (iamroot) then
00590       write(logunit,F02) " maximum           difference for ",trim(attr),tot_diff
00591       write(logunit,F02) " maximum allowable difference for ",trim(attr),eps
00592       call shr_sys_flush(logunit)
00593    endif
00594    call mpi_barrier(mpicom,ier)
00595 
00596    if (ndiff > 0) then
00597       write(logunit,*) trim(subname)," ERROR: incompatible domain grid coordinates"
00598       call shr_sys_flush(logunit)
00599       call mct_die(subName,"incompatible domain grid coordinates")
00600    endif
00601    
00602    deallocate(data1,stat=rcode)
00603    if(rcode /= 0) call mct_die(subName,'deallocate data1')
00604    deallocate(data2,stat=rcode)
00605    if(rcode /= 0) call mct_die(subName,'deallocate data2')
00606    deallocate(lmask,stat=rcode)
00607    if(rcode /= 0) call mct_die(subName,'deallocate lmask')
00608 
00609  end subroutine seq_domain_check_grid_mct
00610 
00611 !===============================================================================
00612 
00613  subroutine domain_areafactinit_mct( cdata, mdl2drv, drv2mdl, comment)
00614     !-----------------------------------------------------------
00615     !
00616     ! Arguments
00617     !
00618     type(seq_cdata) , intent(inout) :: cdata
00619     real(R8),pointer                :: mdl2drv(:)
00620     real(R8),pointer                :: drv2mdl(:)
00621     character(len=*),optional,intent(in) :: comment
00622     !
00623     ! Local variables
00624     !
00625     type(seq_infodata_type),pointer :: infodata
00626     type(mct_gGrid),pointer:: domain
00627     integer                :: ID
00628     integer                :: mpicom
00629     logical                :: iamroot
00630     logical                :: samegrid_ao, samegrid_al
00631     integer                :: j1,j2,m1,n,rcode
00632     integer                :: gridsize,m2dsize,d2msize
00633     real(r8)               :: rmin1,rmax1,rmin,rmax
00634     real(r8)               :: rmask,rarea,raream
00635     character(cl)          :: lcomment
00636     character(len=*),parameter :: subName = '(domain_areafactinit_mct) '
00637     character(len=*),parameter :: F0R = "(2A,2g23.15,A )"
00638     !
00639     !-----------------------------------------------------------
00640 
00641     lcomment = ''
00642     if (present(comment)) lcomment = comment
00643 
00644     call seq_cdata_setptrs(cdata, ID=ID, dom=domain, infodata=infodata)
00645     call seq_comm_setptrs(ID, mpicom=mpicom, iamroot=iamroot) 
00646     call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al)
00647 
00648     ! get sizes
00649 
00650     gridsize = mct_gGrid_lsize(domain)
00651     allocate(drv2mdl(gridsize),mdl2drv(gridsize),stat=rcode)
00652     if(rcode /= 0) call mct_die(subname,'allocate area correction factors')
00653 
00654     j1 = mct_gGrid_indexRA(domain,"area"    ,dieWith=subName)
00655     j2 = mct_gGrid_indexRA(domain,"aream"   ,dieWith=subName)
00656     m1 = mct_gGrid_indexRA(domain,"mask"    ,dieWith=subName)
00657 
00658     mdl2drv(:)=1.0_R8
00659     drv2mdl(:)=1.0_R8
00660 
00661     if (samegrid_ao .and. samegrid_al) then
00662         ! default 1.0
00663     else
00664        do n=1,gridsize
00665           rmask  = domain%data%rAttr(m1,n)
00666           rarea  = domain%data%rAttr(j1,n)
00667           raream = domain%data%rAttr(j2,n)
00668           if ( abs(rmask) >= 1.0e-06) then
00669              if (rarea * raream /= 0.0_r8) then
00670                 mdl2drv(n) = rarea/raream
00671                 drv2mdl(n) = 1.0_R8/mdl2drv(n)
00672                 !if (mdl2drv(n) > 10.0 .or. mdl2drv(n) < 0.1) then
00673                 !   write(logunit,*) trim(subname),' WARNING area,aream= ', &
00674                 !      domain%data%rAttr(j1,n),domain%data%rAttr(j2,n),' in ',n,gridsize
00675                 !endif
00676              else
00677                 write(logunit,*) trim(subname),' ERROR area,aream= ', &
00678                    rarea,raream,' in ',n,gridsize
00679                 call shr_sys_flush(logunit)
00680                 call shr_sys_abort()
00681              endif
00682           endif
00683        enddo
00684     end if
00685        
00686     rmin1 = minval(mdl2drv)
00687     rmax1 = maxval(mdl2drv)
00688     call shr_mpi_min(rmin1,rmin,mpicom)
00689     call shr_mpi_max(rmax1,rmax,mpicom)
00690     if (iamroot) write(logunit,F0R) trim(subname),' : min/max mdl2drv ',rmin,rmax,trim(lcomment)
00691 
00692     rmin1 = minval(drv2mdl)
00693     rmax1 = maxval(drv2mdl)
00694     call shr_mpi_min(rmin1,rmin,mpicom)
00695     call shr_mpi_max(rmax1,rmax,mpicom)
00696     if (iamroot) write(logunit,F0R) trim(subname),' : min/max drv2mdl ',rmin,rmax,trim(lcomment)
00697 
00698  end subroutine domain_areafactinit_mct
00699 
00700 !===============================================================================
00701 
00702 end module seq_domain_mct
00703 
00704 
00705