module seq_domain_mct 1,13
use shr_kind_mod
, only: R8=>shr_kind_r8, IN=>shr_kind_in
use shr_kind_mod, only: CL=>shr_kind_cl
use shr_sys_mod, only: shr_sys_flush, shr_sys_abort
use shr_mpi_mod, only: shr_mpi_min, shr_mpi_max
use mct_mod
use seq_cdata_mod
use seq_comm_mct
use seq_infodata_mod
use map_atmlnd_mct
use map_atmocn_mct
use map_atmice_mct
use map_iceocn_mct
use map_snoglc_mct
implicit none
save
private ! except
!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------
public :: seq_domain_check_mct
public :: domain_areafactinit_mct
!--------------------------------------------------------------------------
! Public variables
!--------------------------------------------------------------------------
real(R8), parameter :: eps_tiny = 1.0e-16_R8 ! roundoff eps
real(R8), parameter :: eps_big = 1.0e+02_R8 ! big eps
real(R8), parameter :: eps_frac_samegrid = 1.0e-14_R8 ! epsilon for fractions for samegrid
!--------------------------------------------------------------------------
! Private interfaces
!--------------------------------------------------------------------------
#ifdef CPP_VECTOR
logical :: usevector = .true.
#else
logical :: usevector = .false.
#endif
#ifdef SYSUNICOS
logical :: usealltoall = .true.
#else
logical :: usealltoall = .false.
#endif
private :: seq_domain_check_grid_mct
#include <mpif.h>
!================================================================================
contains
!================================================================================
!================================================================================
subroutine seq_domain_check_mct( cdata_a, cdata_i, cdata_l, cdata_o, cdata_r, cdata_g, cdata_s) 1,49
!-----------------------------------------------------------
!
! Arguments
!
type(seq_cdata), intent(in) :: cdata_a
type(seq_cdata), intent(in) :: cdata_i
type(seq_cdata), intent(in) :: cdata_l
type(seq_cdata), intent(in) :: cdata_o
type(seq_cdata), intent(in) :: cdata_r
type(seq_cdata), intent(in) :: cdata_g
type(seq_cdata), intent(in) :: cdata_s
!
! Local variables
!
type(mct_gGrid) , pointer :: atmdom_a ! atm domain
type(mct_gGrid) , pointer :: icedom_i ! ice domain
type(mct_gGrid) , pointer :: lnddom_l ! lnd domain
type(mct_gGrid) , pointer :: ocndom_o ! ocn domain
type(mct_gGrid) , pointer :: glcdom_g ! glc domain
type(mct_gGrid) , pointer :: snodom_s ! sno domain
!
type(mct_gsMap) , pointer :: gsMap_a ! atm global seg map
type(mct_gsMap) , pointer :: gsMap_i ! ice global seg map
type(mct_gsMap) , pointer :: gsMap_l ! lnd global seg map
type(mct_gsMap) , pointer :: gsMap_o ! ocn global seg map
type(mct_gsMap) , pointer :: gsMap_r ! ocn global seg map
type(mct_gsMap) , pointer :: gsMap_g ! glc global seg map
type(mct_gsMap) , pointer :: gsMap_s ! sno global seg map
!
type(mct_gGrid) :: lnddom_a ! lnd domain info on atm decomp
type(mct_gGrid) :: icedom_a ! ice domain info on atm decomp (all grids same)
type(mct_gGrid) :: ocndom_a ! ocn domain info on atm decomp (all grids same)
type(mct_gGrid) :: icedom_o ! ocn domain info on ocn decomp (atm/ocn grid different)
type(mct_gGrid) :: snodom_g ! sno domain info on glc decomp
!
integer(IN) :: mpicom_a ! atm mpicom
integer(IN) :: mpicom_i ! ice mpicom
integer(IN) :: mpicom_l ! lnd mpicom
integer(IN) :: mpicom_o ! ocn mpicom
integer(IN) :: mpicom_g ! glc mpicom
integer(IN) :: mpicom_s ! sno mpicom
!
type(seq_infodata_type), pointer :: infodata
!
real(R8), pointer :: fracl(:) ! land fraction on atm decomp
real(R8), pointer :: fraco(:) ! ocn fraction on atm decomp
real(R8), pointer :: fraci(:) ! ice fraction on atm decomp
real(R8), pointer :: maskl(:) ! land mask on atm decomp (all grids same)
real(R8), pointer :: maski(:) ! ice mask on atm decomp (all grids same)
real(R8), pointer :: masko(:) ! ocn mask on atm decomp (all grids same)
!
integer(IN) :: n, kl, ko, ki ! indicies
integer(IN) :: k1,k2,k3 ! indicies
!
logical :: lnd_present ! lnd present flag
logical :: ocn_present ! ocn present flag
logical :: ice_present ! ice present flag
logical :: glc_present ! glc present flag
logical :: sno_present ! sno present flag
logical :: rof_present ! rof present flag
logical :: ocnrof_prognostic ! ocn rof prognostic flag
logical :: samegrid_ao ! atm ocn grid same
logical :: samegrid_ro ! rof ocn grid same
logical :: samegrid_al ! atm lnd grid same
integer(IN) :: rcode ! error status
integer(IN) :: atmsize ! local size of atm grid
integer(IN) :: lndsize ! local size of land grid
integer(IN) :: ocnsize ! local size of ocn grid
integer(IN) :: icesize ! local size of ice grid
integer(IN) :: glcsize ! local size of glc grid
integer(IN) :: snosize ! local size of sno grid
integer(IN) :: gatmsize ! global size of atm grid
integer(IN) :: glndsize ! global size of land grid
integer(IN) :: gocnsize ! global size of ocn grid
integer(IN) :: grofsize ! global size of ocn grid
integer(IN) :: gicesize ! global size of ice grid
integer(IN) :: gglcsize ! global size of glc grid
integer(IN) :: gsnosize ! global size of sno grid
integer(IN) :: npts ! local size temporary
integer(IN) :: ier ! error code
real(R8) :: diff,dmaxo,dmaxi ! difference tracker
logical :: iamroot ! local masterproc
real(R8) :: eps_frac ! epsilon for fractions
real(R8) :: eps_axmask ! epsilon for masks, atm/lnd
real(R8) :: eps_axgrid ! epsilon for grid coords, atm/lnd
real(R8) :: eps_axarea ! epsilon for areas, atm/lnd
real(R8) :: eps_oimask ! epsilon for masks, ocn/ice
real(R8) :: eps_oigrid ! epsilon for grid coords, ocn/ice
real(R8) :: eps_oiarea ! epsilon for areas, ocn/ice
real(R8) :: my_eps_frac ! local eps_frac value
real(r8) :: rmin1,rmax1,rmin,rmax ! local min max computation
!
real(R8),allocatable :: mask (:) ! temporary real vector, domain mask
!
character(*),parameter :: F00 = "('(domain_check_mct) ',4a)"
character(*),parameter :: F01 = "('(domain_check_mct) ',a,i6,a)"
character(*),parameter :: F02 = "('(domain_check_mct) ',a,g23.15)"
character(*),parameter :: F0R = "('(domain_check_mct) ',2A,2g23.15,A )"
character(*),parameter :: subName = '(domain_check_mct) '
!-----------------------------------------------------------
call seq_comm_setptrs(CPLID,iamroot=iamroot)
call seq_cdata_setptrs(cdata_a, infodata=infodata)
call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al, &
samegrid_ro=samegrid_ro)
call seq_infodata_GetData( infodata, lnd_present=lnd_present, ocn_present=ocn_present, &
ice_present=ice_present, glc_present=glc_present, sno_present=sno_present, &
rof_present=rof_present, ocnrof_prognostic=ocnrof_prognostic)
call seq_infodata_GetData( infodata, eps_frac=eps_frac, &
eps_amask=eps_axmask, eps_agrid=eps_axgrid, eps_aarea=eps_axarea, &
eps_omask=eps_oimask, eps_ogrid=eps_oigrid, eps_oarea=eps_oiarea )
! Get info
call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=atmdom_a, mpicom=mpicom_a)
atmsize = mct_avect_lsize(atmdom_a%data)
gatmsize = mct_gsMap_gsize(gsMap_a)
if (lnd_present) then
call seq_cdata_setptrs(cdata_l, gsMap=gsMap_l, dom=lnddom_l, mpicom=mpicom_l)
lndsize = mct_avect_lsize(lnddom_l%data)
glndsize = mct_gsMap_gsize(gsMap_l)
if (samegrid_al .and. gatmsize /= glndsize) then
write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global lndsize= ',glndsize
call shr_sys_flush(logunit)
call mct_die(subname,' atm and lnd grid must have the same global size')
end if
call mct_gGrid_init(oGGrid=lnddom_a, iGGrid=lnddom_l, lsize=atmsize)
call mct_aVect_zero(lnddom_a%data)
call map_lnd2atm_mct(cdata_l, lnddom_l%data, cdata_a, lnddom_a%data)
allocate(maskl(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate maskl')
allocate(fracl(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate fracl')
call mct_aVect_exportRAttr(lnddom_a%data, 'mask', maskl, atmsize)
call mct_aVect_exportRAttr(lnddom_a%data, 'frac', fracl, atmsize)
endif
if (ocn_present) then
call seq_cdata_setptrs(cdata_o, gsMap=gsMap_o, dom=ocndom_o, mpicom=mpicom_o)
ocnsize = mct_avect_lsize(ocndom_o%data)
gocnsize = mct_gsMap_gsize(gsMap_o)
if (samegrid_ao .and. gatmsize /= gocnsize) then
write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global ocnsize= ',gocnsize
call shr_sys_flush(logunit)
call mct_die(subname,' atm and ocn grid must have the same global size')
end if
call mct_gGrid_init(oGGrid=ocndom_a, iGGrid=ocndom_o, lsize=atmsize)
call mct_aVect_zero(ocndom_a%data)
call map_ocn2atm_mct(cdata_o, ocndom_o%data, cdata_a, ocndom_a%data)
allocate(masko(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate masko')
allocate(fraco(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate fraco')
call mct_aVect_exportRAttr(ocndom_a%data, 'mask', masko, atmsize)
if (samegrid_ao) then
call mct_aVect_exportRattr(ocndom_a%data, 'frac', fraco, atmsize)
else
call mct_aVect_exportRattr(ocndom_a%data, 'mask', fraco, atmsize)
endif
endif
if (ice_present) then
call seq_cdata_setptrs(cdata_i, gsMap=gsMap_i, dom=icedom_i, mpicom=mpicom_i)
icesize = mct_avect_lsize(icedom_i%data)
gicesize = mct_gsMap_gsize(gsMap_i)
if (samegrid_ao .and. gatmsize /= gicesize) then
write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global icesize= ',gicesize
call shr_sys_flush(logunit)
call mct_die(subname,' atm and ice grid must have the same global size')
end if
call mct_gGrid_init(oGGrid=icedom_a, iGGrid=icedom_i, lsize=atmsize)
call mct_aVect_zero(icedom_a%data)
call map_ice2atm_mct(cdata_i, icedom_i%data, cdata_a, icedom_a%data)
allocate(maski(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate maski')
allocate(fraci(atmsize),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate fraci')
call mct_aVect_exportRAttr(icedom_a%data, 'mask', maski, atmsize)
if (samegrid_ao) then
call mct_aVect_exportRattr(icedom_a%data, 'frac', fraci, atmsize)
else
call mct_aVect_exportRattr(icedom_a%data, 'mask', fraci, atmsize)
endif
endif
if (sno_present .and. glc_present) then
call seq_cdata_setptrs(cdata_g, gsMap=gsMap_g, dom=glcdom_g, mpicom=mpicom_g)
glcsize = mct_avect_lsize(glcdom_g%data)
gglcsize = mct_gsMap_gsize(gsMap_g)
call seq_cdata_setptrs(cdata_s, gsMap=gsMap_s, dom=snodom_s, mpicom=mpicom_s)
snosize = mct_avect_lsize(snodom_s%data)
gsnosize = mct_gsMap_gsize(gsMap_s)
if (gglcsize /= gsnosize) then
write(logunit,*) subname,' error: global glcsize = ',gglcsize,' global snosize= ',gsnosize
call shr_sys_flush(logunit)
call mct_die(subname,' glc and sno grid must have the same global size')
end if
call mct_gGrid_init(oGGrid=snodom_g, iGGrid=snodom_s, lsize=glcsize)
call mct_aVect_zero(snodom_g%data)
call map_sno2glc_mct(cdata_s, snodom_s%data, cdata_g, snodom_g%data)
if (iamroot) write(logunit,F00) ' --- checking glc/sno domains ---'
npts = glcsize
allocate(mask(npts),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate mask')
call mct_aVect_getRAttr(snodom_g%data,"mask",mask,rcode)
where (mask < eps_axmask) mask = 0.0_r8
call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'mask', eps=eps_axmask, mpicom=mpicom_g, mask=mask)
call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'area', eps=eps_axarea, mpicom=mpicom_g, mask=mask)
deallocate(mask,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate mask')
endif
if (ice_present .and. ocn_present) then
if (gocnsize /= gicesize) then
write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global icesize= ',gicesize
call shr_sys_flush(logunit)
call mct_die(subname,' ocean and ice grid must have the same global size')
endif
call mct_gGrid_init(oGGrid=icedom_o, iGGrid=icedom_i, lsize=ocnsize)
call mct_aVect_zero(icedom_o%data)
call map_ice2ocn_mct(cdata_i, icedom_i%data, cdata_o, icedom_o%data)
end if
if (rof_present .and. ocnrof_prognostic .and. samegrid_ro) then
call seq_cdata_setptrs(cdata_r, gsMap=gsMap_r)
grofsize = mct_gsMap_gsize(gsMap_r)
if (gocnsize /= grofsize) then
write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global rofsize= ',grofsize
call shr_sys_flush(logunit)
call mct_die(subname,' ocean and rof grid must have the same global size')
endif
end if
!------------------------------------------------------------------------------
! Check ice/ocean grid consistency
!------------------------------------------------------------------------------
if (ocn_present .and. ice_present) then
! if (samegrid_oi) then ! doesn't yet exist
npts = ocnsize
allocate(mask(npts),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate mask')
if (iamroot) write(logunit,F00) ' --- checking ocn/ice domains ---'
call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'mask', eps=eps_oigrid, mpicom=mpicom_o)
call mct_aVect_getRAttr(ocndom_o%data,"mask",mask,rcode)
where (mask < eps_oimask) mask = 0.0_r8
call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lat' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lon' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'area', eps=eps_oiarea, mpicom=mpicom_o, mask=mask)
deallocate(mask,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate mask')
! endif
endif
!------------------------------------------------------------------------------
! Check atm/lnd grid consistency
!------------------------------------------------------------------------------
if (lnd_present .and. samegrid_al) then
if (iamroot) write(logunit,F00) ' --- checking atm/land domains ---'
call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maskl)
endif
!------------------------------------------------------------------------------
! Check atm/ocn and atm/ice grid consistency (if samegrid)
!------------------------------------------------------------------------------
if (ice_present .and. samegrid_ao) then
if (iamroot) write(logunit,F00) ' --- checking atm/ice domains ---'
call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maski)
endif
if (ocn_present .and. samegrid_ao) then
if (iamroot) write(logunit,F00) ' --- checking atm/ocn domains ---'
call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=masko)
endif
!------------------------------------------------------------------------------
! Check consistency of land fraction with ocean mask on grid
!------------------------------------------------------------------------------
my_eps_frac = eps_frac
if (samegrid_ao) my_eps_frac = eps_frac_samegrid
if (.not. samegrid_al) my_eps_frac = eps_big
if (iamroot) write(logunit,F00) ' --- checking fractions in domains ---'
dmaxi = 0.0_R8
dmaxo = 0.0_R8
do n = 1,atmsize
if (lnd_present .and. ice_present) then
diff = abs(1._r8 - fracl(n) - fraci(n))
dmaxi = max(diff,dmaxi)
if (diff > my_eps_frac) then
write(logunit,*)'inconsistency between land fraction and ice land fraction'
write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n),' sum= ',fracl(n)+fraci(n)
call shr_sys_flush(logunit)
call mct_die(subname,'inconsistency between land fraction and ice land fraction')
end if
if ((1._r8-fraci(n)) > eps_frac .and. fracl(n) < eps_tiny) then
write(logunit,*)'inconsistency between land mask and ice land mask'
write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n)
call shr_sys_flush(logunit)
call mct_die(subname,' inconsistency between land mask and ice land mask')
end if
endif
if (lnd_present .and. ocn_present) then
diff = abs(1._r8 - fracl(n) - fraco(n))
dmaxo = max(diff,dmaxo)
if (diff > my_eps_frac) then
write(logunit,*)'inconsistency between land fraction and ocn land fraction'
write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n),' sum= ',fracl(n)+fraco(n)
call shr_sys_flush(logunit)
call mct_die(subname,' inconsistency between land fraction and ocn land fraction')
end if
if ((1._r8-fraco(n)) > eps_frac .and. fracl(n) < eps_tiny) then
write(logunit,*)'inconsistency between land mask and ocn land mask'
write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n)
call shr_sys_flush(logunit)
call mct_die(subname,' inconsistency between land mask and ocn land mask')
end if
endif
end do
if (iamroot) then
write(logunit,F02) ' maximum difference for ofrac sum ',dmaxo
write(logunit,F02) ' maximum difference for ifrac sum ',dmaxi
write(logunit,F02) ' maximum allowable difference for frac sum ',my_eps_frac
write(logunit,F02) ' maximum allowable tolerance for valid frac ',eps_frac
call shr_sys_flush(logunit)
endif
!------------------------------------------------------------------------------
! Set atm and lnd ascale
!------------------------------------------------------------------------------
if (lnd_present .and. ocn_present) then
k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
do k2 = 1,atmsize
if (fracl(k2) /= 0.0_r8) then
atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraco(k2))/(fracl(k2))
else
atmdom_a%data%rAttr(k1,k2) = 0.0
endif
enddo
call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
elseif (lnd_present .and. ice_present) then
k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
do k2 = 1,atmsize
if (fracl(k2) /= 0.0_r8) then
atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraci(k2))/(fracl(k2))
else
atmdom_a%data%rAttr(k1,k2) = 0.0
endif
enddo
call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
endif
if (lnd_present .and. (ocn_present .or. ice_present)) then
k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check atm ascale')
rmin1 = minval(atmdom_a%data%rAttr(k1,:))
rmax1 = maxval(atmdom_a%data%rAttr(k1,:))
call shr_mpi_min(rmin1,rmin,mpicom_a)
call shr_mpi_max(rmax1,rmax,mpicom_a)
if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' atmdom_a'
k1 = mct_aVect_indexRa(lnddom_l%data,"ascale",perrWith='domain_check lnd ascale')
rmin1 = minval(lnddom_l%data%rAttr(k1,:))
rmax1 = maxval(lnddom_l%data%rAttr(k1,:))
call shr_mpi_min(rmin1,rmin,mpicom_l)
call shr_mpi_max(rmax1,rmax,mpicom_l)
if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' lnddom_l'
endif
!------------------------------------------------------------------------------
! Clean up allocated memory
!------------------------------------------------------------------------------
if (lnd_present) then
deallocate(fracl,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate fracl')
deallocate(maskl,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate maskl')
call mct_gGrid_clean(lnddom_a, rcode)
if(rcode /= 0) call mct_die(subName,'clean lnddom_a')
endif
if (ocn_present) then
deallocate(fraco,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate fraco')
deallocate(masko,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate masko')
call mct_gGrid_clean(ocndom_a, rcode)
if(rcode /= 0) call mct_die(subName,'clean ocndom_a')
endif
if (ice_present) then
deallocate(fraci,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate fraci')
deallocate(maski,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate maski')
call mct_gGrid_clean(icedom_a, rcode)
if(rcode /= 0) call mct_die(subName,'clean icedom_o')
endif
if (ocn_present .and. ice_present) then
call mct_gGrid_clean(icedom_o, rcode)
if(rcode /= 0) call mct_die(subName,'clean icedom_o')
endif
end subroutine seq_domain_check_mct
!===============================================================================
subroutine seq_domain_check_grid_mct(dom1, dom2, attr, eps, mpicom, mask) 17,3
!-----------------------------------------------------------
! Arguments
type(mct_aVect) , intent(in) :: dom1
type(mct_aVect) , intent(in) :: dom2
character(len=*), intent(in) :: attr ! grid attribute to compare
real(R8) , intent(in) :: eps ! error condition for compare
integer(IN) , intent(in) :: mpicom
real(R8) , intent(in), optional :: mask(:)
! Local variables
integer(in) :: n,ndiff ! indices
integer(in) :: npts1,npts2,npts ! counters
integer(in) :: rcode ! error code
real(R8) :: diff,max_diff ! temporaries
real(R8) :: tot_diff ! maximum diff across all pes
integer(IN) :: ier ! error code
real(R8), pointer :: data1(:) ! temporaries
real(R8), pointer :: data2(:) ! temporaries
real(R8), pointer :: lmask(:) ! temporaries
logical :: iamroot ! local masterproc
character(*),parameter :: F00 = "('(domain_check_grid_mct) ',4a)"
character(*),parameter :: F01 = "('(domain_check_grid_mct) ',a,i12,a)"
character(*),parameter :: F02 = "('(domain_check_grid_mct) ',2a,g23.15)"
character(*),parameter :: F0R = "('(domain_check_grid_mct) ',2A,2g23.15,A )"
character(*),parameter :: subName = '(domain_check_grid_mct) '
!-----------------------------------------------------------
call seq_comm_setptrs(CPLID,iamroot=iamroot)
npts1 = mct_aVect_lsize(dom1)
npts2 = mct_aVect_lsize(dom2)
npts = npts1
if (npts1 == npts2) then
if (iamroot) write(logunit,F01) " the domain size is = ", npts
else
write(logunit,*) trim(subname)," domain size #1 = ", npts1
write(logunit,*) trim(subname)," domain size #2 = ", npts2
write(logunit,*) trim(subname)," ERROR: domain size mis-match"
call mct_die(subName,"ERROR: domain size mis-match")
end if
allocate(data1(npts),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate data1')
allocate(data2(npts),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate data2')
allocate(lmask(npts),stat=rcode)
if(rcode /= 0) call mct_die(subName,'allocate lmask')
call mct_aVect_exportRAttr(dom1, trim(attr), data1, npts)
call mct_aVect_exportRAttr(dom2, trim(attr), data2, npts)
lmask = 1.0_r8
if (present(mask)) then
if (size(mask) /= npts) then
call mct_die(subName,"ERROR: mask size mis-match")
endif
lmask = mask
endif
! --- adjust lons to address wraparound issues, we're assuming degree here! ---
if (trim(attr) == "lon") then
do n = 1,npts
if (data2(n) > data1(n)) then
do while ( (data1(n)+360.0_R8) < (data2(n)+180.0_R8) ) ! longitude is periodic
data1(n) = data1(n) + 360.0_R8
end do
else
do while ( (data2(n)+360.0_R8) < (data1(n)+180.0_R8) ) ! longitude is periodic
data2(n) = data2(n) + 360.0_R8
end do
endif
enddo
endif
! Only check consistency where mask is greater than zero, if mask is present
max_diff = 0.0_R8
ndiff = 0
do n=1,npts
if (lmask(n) > eps_tiny) then
diff = abs(data1(n)-data2(n))
max_diff = max(max_diff,diff)
if (diff > eps) then
!debug write(logunit,*)'n= ',n,' data1= ',data1(n),' data2= ',data2(n),' diff= ',diff, ' eps= ',eps
ndiff = ndiff + 1
endif
end if
end do
call mpi_reduce(max_diff,tot_diff,1,MPI_REAL8,MPI_MAX,0,mpicom,ier)
if (iamroot) then
write(logunit,F02) " maximum difference for ",trim(attr),tot_diff
write(logunit,F02) " maximum allowable difference for ",trim(attr),eps
call shr_sys_flush(logunit)
endif
call mpi_barrier(mpicom,ier)
if (ndiff > 0) then
write(logunit,*) trim(subname)," ERROR: incompatible domain grid coordinates"
call shr_sys_flush(logunit)
call mct_die(subName,"incompatible domain grid coordinates")
endif
deallocate(data1,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate data1')
deallocate(data2,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate data2')
deallocate(lmask,stat=rcode)
if(rcode /= 0) call mct_die(subName,'deallocate lmask')
end subroutine seq_domain_check_grid_mct
!===============================================================================
subroutine domain_areafactinit_mct( cdata, mdl2drv, drv2mdl, comment) 7,5
!-----------------------------------------------------------
!
! Arguments
!
type(seq_cdata) , intent(inout) :: cdata
real(R8),pointer :: mdl2drv(:)
real(R8),pointer :: drv2mdl(:)
character(len=*),optional,intent(in) :: comment
!
! Local variables
!
type(seq_infodata_type),pointer :: infodata
type(mct_gGrid),pointer:: domain
integer :: ID
integer :: mpicom
logical :: iamroot
logical :: samegrid_ao, samegrid_al
integer :: j1,j2,m1,n,rcode
integer :: gridsize,m2dsize,d2msize
real(r8) :: rmin1,rmax1,rmin,rmax
real(r8) :: rmask,rarea,raream
character(cl) :: lcomment
character(len=*),parameter :: subName = '(domain_areafactinit_mct) '
character(len=*),parameter :: F0R = "(2A,2g23.15,A )"
!
!-----------------------------------------------------------
lcomment = ''
if (present(comment)) lcomment = comment
call seq_cdata_setptrs(cdata, ID=ID, dom=domain, infodata=infodata)
call seq_comm_setptrs(ID, mpicom=mpicom, iamroot=iamroot)
call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al)
! get sizes
gridsize = mct_gGrid_lsize(domain)
allocate(drv2mdl(gridsize),mdl2drv(gridsize),stat=rcode)
if(rcode /= 0) call mct_die(subname,'allocate area correction factors')
j1 = mct_gGrid_indexRA(domain,"area" ,dieWith=subName)
j2 = mct_gGrid_indexRA(domain,"aream" ,dieWith=subName)
m1 = mct_gGrid_indexRA(domain,"mask" ,dieWith=subName)
mdl2drv(:)=1.0_R8
drv2mdl(:)=1.0_R8
if (samegrid_ao .and. samegrid_al) then
! default 1.0
else
do n=1,gridsize
rmask = domain%data%rAttr(m1,n)
rarea = domain%data%rAttr(j1,n)
raream = domain%data%rAttr(j2,n)
if ( abs(rmask) >= 1.0e-06) then
if (rarea * raream /= 0.0_r8) then
mdl2drv(n) = rarea/raream
drv2mdl(n) = 1.0_R8/mdl2drv(n)
!if (mdl2drv(n) > 10.0 .or. mdl2drv(n) < 0.1) then
! write(logunit,*) trim(subname),' WARNING area,aream= ', &
! domain%data%rAttr(j1,n),domain%data%rAttr(j2,n),' in ',n,gridsize
!endif
else
write(logunit,*) trim(subname),' ERROR area,aream= ', &
rarea,raream,' in ',n,gridsize
call shr_sys_flush(logunit)
call shr_sys_abort()
endif
endif
enddo
end if
rmin1 = minval(mdl2drv)
rmax1 = maxval(mdl2drv)
call shr_mpi_min(rmin1,rmin,mpicom)
call shr_mpi_max(rmax1,rmax,mpicom)
if (iamroot) write(logunit,F0R) trim(subname),' : min/max mdl2drv ',rmin,rmax,trim(lcomment)
rmin1 = minval(drv2mdl)
rmax1 = maxval(drv2mdl)
call shr_mpi_min(rmin1,rmin,mpicom)
call shr_mpi_max(rmax1,rmax,mpicom)
if (iamroot) write(logunit,F0R) trim(subname),' : min/max drv2mdl ',rmin,rmax,trim(lcomment)
end subroutine domain_areafactinit_mct
!===============================================================================
end module seq_domain_mct