Revision f7b2f3fd43c140a005e5d215dfcbc131fe9c1c0f authored by Ryan Knox on 12 October 2021, 14:49:52 UTC, committed by Ryan Knox on 12 October 2021, 14:49:52 UTC
1 parent bf910e9
seq_map_mod.F90
module seq_map_mod
!---------------------------------------------------------------------
!
! Purpose:
!
! General mapping routines
! including self normalizing mapping routine with optional fraction
!
! Author: T. Craig, Jan-28-2011
!
!---------------------------------------------------------------------
use shr_kind_mod ,only: R8 => SHR_KIND_R8, IN=>SHR_KIND_IN
use shr_kind_mod ,only: CL => SHR_KIND_CL, CX => SHR_KIND_CX
use shr_sys_mod
use shr_const_mod
use shr_mct_mod, only: shr_mct_sMatPInitnc, shr_mct_queryConfigFile
use mct_mod
use seq_comm_mct
use component_type_mod
use seq_map_type_mod
implicit none
save
private ! except
!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------
public :: seq_map_init_rcfile ! cpl pes
public :: seq_map_init_rearrolap ! cpl pes
public :: seq_map_initvect ! cpl pes
public :: seq_map_map ! cpl pes
public :: seq_map_mapvect ! cpl pes
public :: seq_map_readdata ! cpl pes
interface seq_map_avNorm
module procedure seq_map_avNormArr
module procedure seq_map_avNormAvF
end interface seq_map_avNorm
!--------------------------------------------------------------------------
! Public data
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------
character(*),parameter :: seq_map_stroff = 'variable_unset'
character(*),parameter :: seq_map_stron = 'StrinG_is_ON'
real(R8),parameter,private :: deg2rad = shr_const_pi/180.0_R8 ! deg to rads
!=======================================================================
contains
!=======================================================================
subroutine seq_map_init_rcfile( mapper, comp_s, comp_d, &
maprcfile, maprcname, maprctype, samegrid, string, esmf_map)
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout),pointer :: mapper
type(component_type) ,intent(inout) :: comp_s
type(component_type) ,intent(inout) :: comp_d
character(len=*) ,intent(in) :: maprcfile
character(len=*) ,intent(in) :: maprcname
character(len=*) ,intent(in) :: maprctype
logical ,intent(in) :: samegrid
character(len=*) ,intent(in),optional :: string
logical ,intent(in),optional :: esmf_map
!
! Local Variables
!
type(mct_gsmap), pointer :: gsmap_s ! temporary pointers
type(mct_gsmap), pointer :: gsmap_d ! temporary pointers
integer(IN) :: mpicom
character(CX) :: mapfile
character(CL) :: maptype
integer(IN) :: mapid
character(len=*),parameter :: subname = "(seq_map_init_rcfile) "
!-----------------------------------------------------
if (seq_comm_iamroot(CPLID) .and. present(string)) then
write(logunit,'(A)') subname//' called for '//trim(string)
endif
call seq_comm_setptrs(CPLID, mpicom=mpicom)
gsmap_s => component_get_gsmap_cx(comp_s)
gsmap_d => component_get_gsmap_cx(comp_d)
if (mct_gsmap_Identical(gsmap_s,gsmap_d)) then
call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="copy")
if (mapid > 0) then
call seq_map_mappoint(mapid,mapper)
else
call seq_map_mapinit(mapper,mpicom)
mapper%copy_only = .true.
mapper%strategy = "copy"
mapper%gsmap_s => component_get_gsmap_cx(comp_s)
mapper%gsmap_d => component_get_gsmap_cx(comp_d)
endif
elseif (samegrid) then
call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="rearrange")
if (mapid > 0) then
call seq_map_mappoint(mapid,mapper)
else
! --- Initialize rearranger
call seq_map_mapinit(mapper,mpicom)
mapper%rearrange_only = .true.
mapper%strategy = "rearrange"
mapper%gsmap_s => component_get_gsmap_cx(comp_s)
mapper%gsmap_d => component_get_gsmap_cx(comp_d)
call seq_map_gsmapcheck(gsmap_s, gsmap_d)
call mct_rearr_init(gsmap_s, gsmap_d, mpicom, mapper%rearr)
endif
else
! --- Initialize Smatp
call shr_mct_queryConfigFile(mpicom,maprcfile,maprcname,mapfile,maprctype,maptype)
call seq_map_mapmatch(mapid,gsMap_s=gsMap_s,gsMap_d=gsMap_d,mapfile=mapfile,strategy=maptype)
if (mapid > 0) then
call seq_map_mappoint(mapid,mapper)
else
call seq_map_mapinit(mapper,mpicom)
mapper%mapfile = trim(mapfile)
mapper%strategy= trim(maptype)
mapper%gsmap_s => component_get_gsmap_cx(comp_s)
mapper%gsmap_d => component_get_gsmap_cx(comp_d)
call shr_mct_sMatPInitnc(mapper%sMatp, mapper%gsMap_s, mapper%gsMap_d, trim(mapfile),trim(maptype),mpicom)
if (present(esmf_map)) mapper%esmf_map = esmf_map
if (mapper%esmf_map) then
call shr_sys_abort(subname//' ERROR: esmf SMM not supported')
endif ! esmf_map
endif ! mapid >= 0
endif
if (seq_comm_iamroot(CPLID)) then
write(logunit,'(2A,I6,4A)') subname,' mapper counter, strategy, mapfile = ', &
mapper%counter,' ',trim(mapper%strategy),' ',trim(mapper%mapfile)
call shr_sys_flush(logunit)
endif
end subroutine seq_map_init_rcfile
!=======================================================================
subroutine seq_map_init_rearrolap(mapper, comp_s, comp_d, string)
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout),pointer :: mapper
type(component_type) ,intent(inout) :: comp_s
type(component_type) ,intent(inout) :: comp_d
character(len=*) ,intent(in),optional :: string
!
! Local Variables
!
integer(IN) :: mapid
type(mct_gsmap), pointer :: gsmap_s
type(mct_gsmap), pointer :: gsmap_d
integer(IN) :: mpicom
character(len=*),parameter :: subname = "(seq_map_init_rearrolap) "
!-----------------------------------------------------
if (seq_comm_iamroot(CPLID) .and. present(string)) then
write(logunit,'(A)') subname//' called for '//trim(string)
endif
call seq_comm_setptrs(CPLID, mpicom=mpicom)
gsmap_s => component_get_gsmap_cx(comp_s)
gsmap_d => component_get_gsmap_cx(comp_d)
if (mct_gsmap_Identical(gsmap_s,gsmap_d)) then
call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="copy")
if (mapid > 0) then
call seq_map_mappoint(mapid,mapper)
else
call seq_map_mapinit(mapper,mpicom)
mapper%copy_only = .true.
mapper%strategy = "copy"
mapper%gsmap_s => component_get_gsmap_cx(comp_s)
mapper%gsmap_d => component_get_gsmap_cx(comp_d)
endif
else
call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="rearrange")
if (mapid > 0) then
call seq_map_mappoint(mapid,mapper)
else
! --- Initialize rearranger
call seq_map_mapinit(mapper, mpicom)
mapper%rearrange_only = .true.
mapper%strategy = "rearrange"
mapper%gsmap_s => component_get_gsmap_cx(comp_s)
mapper%gsmap_d => component_get_gsmap_cx(comp_d)
call seq_map_gsmapcheck(gsmap_s, gsmap_d)
call mct_rearr_init(gsmap_s, gsmap_d, mpicom, mapper%rearr)
endif
endif
if (seq_comm_iamroot(CPLID)) then
write(logunit,'(2A,I6,4A)') subname,' mapper counter, strategy, mapfile = ', &
mapper%counter,' ',trim(mapper%strategy),' ',trim(mapper%mapfile)
call shr_sys_flush(logunit)
endif
end subroutine seq_map_init_rearrolap
!=======================================================================
subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, &
string, msgtag )
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout) :: mapper
type(mct_aVect) ,intent(in) :: av_s
type(mct_aVect) ,intent(inout) :: av_d
character(len=*),intent(in),optional :: fldlist
logical ,intent(in),optional :: norm
type(mct_aVect) ,intent(in),optional :: avwts_s
character(len=*),intent(in),optional :: avwtsfld_s
character(len=*),intent(in),optional :: string
integer(IN) ,intent(in),optional :: msgtag
!
! Local Variables
!
logical :: lnorm
integer(IN),save :: ltag ! message tag for rearrange
character(len=*),parameter :: subname = "(seq_map_map) "
!-----------------------------------------------------
if (seq_comm_iamroot(CPLID) .and. present(string)) then
write(logunit,'(A)') subname//' called for '//trim(string)
endif
lnorm = .true.
if (present(norm)) then
lnorm = norm
endif
if (present(msgtag)) then
ltag = msgtag
else
ltag = 2000
endif
if (present(avwts_s) .and. .not. present(avwtsfld_s)) then
write(logunit,*) subname,' ERROR: avwts_s present but avwtsfld_s not'
call shr_sys_abort(subname//' ERROR: avwts present')
endif
if (.not. present(avwts_s) .and. present(avwtsfld_s)) then
write(logunit,*) subname,' ERROR: avwtsfld_s present but avwts_s not'
call shr_sys_abort(subname//' ERROR: avwtsfld present')
endif
if (mapper%copy_only) then
!-------------------------------------------
! COPY data
!-------------------------------------------
if (present(fldlist)) then
call mct_aVect_copy(aVin=av_s,aVout=av_d,rList=fldlist,vector=mct_usevector)
else
call mct_aVect_copy(aVin=av_s,aVout=av_d,vector=mct_usevector)
endif
else if (mapper%rearrange_only) then
!-------------------------------------------
! REARRANGE data
!-------------------------------------------
if (present(fldlist)) then
call mct_rearr_rearrange_fldlist(av_s, av_d, mapper%rearr, tag=ltag, VECTOR=mct_usevector, &
ALLTOALL=mct_usealltoall, fldlist=fldlist)
else
call mct_rearr_rearrange(av_s, av_d, mapper%rearr, tag=ltag, VECTOR=mct_usevector, &
ALLTOALL=mct_usealltoall)
endif
else
!-------------------------------------------
! MAP data
!-------------------------------------------
if (present(avwts_s)) then
if (present(fldlist)) then
call seq_map_avNorm(mapper, av_s, av_d, avwts_s, trim(avwtsfld_s), &
rList=fldlist, norm=lnorm)
else
call seq_map_avNorm(mapper, av_s, av_d, avwts_s, trim(avwtsfld_s), &
norm=lnorm)
endif
else
if (present(fldlist)) then
call seq_map_avNorm(mapper, av_s, av_d, rList=fldlist, norm=lnorm)
else
call seq_map_avNorm(mapper, av_s, av_d, norm=lnorm)
endif
endif
end if
end subroutine seq_map_map
!=======================================================================
subroutine seq_map_initvect(mapper, type, comp_s, comp_d, string)
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout) :: mapper
character(len=*) ,intent(in) :: type
type(component_type) ,intent(inout) :: comp_s
type(component_type) ,intent(inout) :: comp_d
character(len=*) ,intent(in),optional :: string
!
! Local Variables
!
type(mct_gGrid), pointer :: dom_s
type(mct_gGrid), pointer :: dom_d
integer(IN) :: klon, klat, lsize, n
character(len=CL) :: lstring
character(len=*),parameter :: subname = "(seq_map_initvect) "
!-----------------------------------------------------
lstring = ' '
if (present(string)) then
if (seq_comm_iamroot(CPLID)) write(logunit,'(A)') subname//' called for '//trim(string)
lstring = trim(string)
endif
dom_s => component_get_dom_cx(comp_s)
dom_d => component_get_dom_cx(comp_d)
if (trim(type(1:6)) == 'cart3d') then
if (mapper%cart3d_init == trim(seq_map_stron)) return
!--- compute these up front for vector mapping ---
lsize = mct_aVect_lsize(dom_s%data)
allocate(mapper%slon_s(lsize),mapper%clon_s(lsize), &
mapper%slat_s(lsize),mapper%clat_s(lsize))
klon = mct_aVect_indexRa(dom_s%data, "lon" )
klat = mct_aVect_indexRa(dom_s%data, "lat" )
do n = 1,lsize
mapper%slon_s(n) = sin(dom_s%data%rAttr(klon,n)*deg2rad)
mapper%clon_s(n) = cos(dom_s%data%rAttr(klon,n)*deg2rad)
mapper%slat_s(n) = sin(dom_s%data%rAttr(klat,n)*deg2rad)
mapper%clat_s(n) = cos(dom_s%data%rAttr(klat,n)*deg2rad)
enddo
lsize = mct_aVect_lsize(dom_d%data)
allocate(mapper%slon_d(lsize),mapper%clon_d(lsize), &
mapper%slat_d(lsize),mapper%clat_d(lsize))
klon = mct_aVect_indexRa(dom_d%data, "lon" )
klat = mct_aVect_indexRa(dom_d%data, "lat" )
do n = 1,lsize
mapper%slon_d(n) = sin(dom_d%data%rAttr(klon,n)*deg2rad)
mapper%clon_d(n) = cos(dom_d%data%rAttr(klon,n)*deg2rad)
mapper%slat_d(n) = sin(dom_d%data%rAttr(klat,n)*deg2rad)
mapper%clat_d(n) = cos(dom_d%data%rAttr(klat,n)*deg2rad)
enddo
mapper%cart3d_init = trim(seq_map_stron)
endif
end subroutine seq_map_initvect
!=======================================================================
subroutine seq_map_mapvect( mapper, type, av_s, av_d, fldu, fldv, norm, string )
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout) :: mapper
character(len=*),intent(in) :: type
type(mct_aVect) ,intent(in) :: av_s
type(mct_aVect) ,intent(inout) :: av_d
character(len=*),intent(in) :: fldu
character(len=*),intent(in) :: fldv
logical ,intent(in),optional :: norm
character(len=*),intent(in),optional :: string
!
! Local Variables
!
logical :: lnorm
character(len=CL) :: lstring
character(len=*),parameter :: subname = "(seq_map_mapvect) "
!-----------------------------------------------------
lstring = ' '
if (present(string)) then
if (seq_comm_iamroot(CPLID)) write(logunit,'(A)') subname//' called for '//trim(string)
lstring = trim(string)
endif
if (mapper%copy_only .or. mapper%rearrange_only) then
return
endif
lnorm = .true.
if (present(norm)) then
lnorm = norm
endif
if (trim(type(1:6)) == 'cart3d') then
if (mapper%cart3d_init /= trim(seq_map_stron)) then
call shr_sys_abort(trim(subname)//' ERROR: cart3d not initialized '//trim(lstring))
endif
call seq_map_cart3d(mapper, type, av_s, av_d, fldu, fldv, norm=lnorm, string=string)
elseif (trim(type) == 'none') then
call seq_map_map(mapper, av_s, av_d, fldlist=trim(fldu)//':'//trim(fldv), norm=lnorm)
else
write(logunit,*) subname,' ERROR: type unsupported ',trim(type)
call shr_sys_abort(trim(subname)//' ERROR in type='//trim(type))
end if
end subroutine seq_map_mapvect
!=======================================================================
subroutine seq_map_cart3d( mapper, type, av_s, av_d, fldu, fldv, norm, string)
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) ,intent(inout) :: mapper
character(len=*),intent(in) :: type
type(mct_aVect) ,intent(in) :: av_s
type(mct_aVect) ,intent(inout) :: av_d
character(len=*),intent(in) :: fldu
character(len=*),intent(in) :: fldv
logical ,intent(in),optional :: norm
character(len=*),intent(in),optional :: string
!
! Local Variables
!
integer :: lsize
logical :: lnorm
integer :: ku,kv,kux,kuy,kuz,n
real(r8) :: ue,un,ur,ux,uy,uz,speed
real(r8) :: urmaxl,urmax,uravgl,uravg,spavgl,spavg
type(mct_aVect) :: av3_s, av3_d
integer(in) :: mpicom,my_task,ierr,urcnt,urcntl
character(len=*),parameter :: subname = "(seq_map_cart3d) "
lnorm = .true.
if (present(norm)) then
lnorm=norm
endif
mpicom = mapper%mpicom
ku = mct_aVect_indexRA(av_s, trim(fldu), perrwith='quiet')
kv = mct_aVect_indexRA(av_s, trim(fldv), perrwith='quiet')
if (ku /= 0 .and. kv /= 0) then
lsize = mct_aVect_lsize(av_s)
call mct_avect_init(av3_s,rList='ux:uy:uz',lsize=lsize)
lsize = mct_aVect_lsize(av_d)
call mct_avect_init(av3_d,rList='ux:uy:uz',lsize=lsize)
kux = mct_aVect_indexRA(av3_s,'ux')
kuy = mct_aVect_indexRA(av3_s,'uy')
kuz = mct_aVect_indexRA(av3_s,'uz')
lsize = mct_aVect_lsize(av_s)
do n = 1,lsize
ur = 0.0_r8
ue = av_s%rAttr(ku,n)
un = av_s%rAttr(kv,n)
ux = mapper%clon_s(n)*mapper%clat_s(n)*ur - &
mapper%clon_s(n)*mapper%slat_s(n)*un - &
mapper%slon_s(n)*ue
uy = mapper%slon_s(n)*mapper%clon_s(n)*ur - &
mapper%slon_s(n)*mapper%slat_s(n)*un + &
mapper%clon_s(n)*ue
uz = mapper%slat_s(n)*ur + &
mapper%clat_s(n)*un
av3_s%rAttr(kux,n) = ux
av3_s%rAttr(kuy,n) = uy
av3_s%rAttr(kuz,n) = uz
enddo
call seq_map_map(mapper, av3_s, av3_d, norm=lnorm)
kux = mct_aVect_indexRA(av3_d,'ux')
kuy = mct_aVect_indexRA(av3_d,'uy')
kuz = mct_aVect_indexRA(av3_d,'uz')
lsize = mct_aVect_lsize(av_d)
urmaxl = -1.0_r8
uravgl = 0.0_r8
urcntl = 0
spavgl = 0.0_r8
do n = 1,lsize
ux = av3_d%rAttr(kux,n)
uy = av3_d%rAttr(kuy,n)
uz = av3_d%rAttr(kuz,n)
ue = -mapper%slon_d(n) *ux + &
mapper%clon_d(n) *uy
un = -mapper%clon_d(n)*mapper%slat_d(n)*ux - &
mapper%slon_d(n)*mapper%slat_d(n)*uy + &
mapper%clat_d(n)*uz
ur = mapper%clon_d(n)*mapper%clat_d(n)*ux + &
mapper%slon_d(n)*mapper%clat_d(n)*uy - &
mapper%slat_d(n)*uz
speed = sqrt(ur*ur + ue*ue + un*un)
if (trim(type) == 'cart3d_diag' .or. trim(type) == 'cart3d_uvw_diag') then
if (speed /= 0.0_r8) then
urmaxl = max(urmaxl,abs(ur))
uravgl = uravgl + abs(ur)
spavgl = spavgl + speed
urcntl = urcntl + 1
endif
endif
if (type(1:10) == 'cart3d_uvw') then
!--- this adds ur to ue and un, while preserving u/v angle and total speed ---
if (un == 0.0_R8) then
!--- if ue is also 0.0 then just give speed to ue, this is arbitrary ---
av_d%rAttr(ku,n) = sign(speed,ue)
av_d%rAttr(kv,n) = 0.0_r8
else if (ue == 0.0_R8) then
av_d%rAttr(ku,n) = 0.0_r8
av_d%rAttr(kv,n) = sign(speed,un)
else
av_d%rAttr(ku,n) = sign(speed/sqrt(1.0_r8 + ((un*un)/(ue*ue))),ue)
av_d%rAttr(kv,n) = sign(speed/sqrt(1.0_r8 + ((ue*ue)/(un*un))),un)
endif
else
!--- this ignores ur ---
av_d%rAttr(ku,n) = ue
av_d%rAttr(kv,n) = un
endif
enddo
if (trim(type) == 'cart3d_diag' .or. trim(type) == 'cart3d_uvw_diag') then
call mpi_comm_rank(mpicom,my_task,ierr)
call shr_mpi_max(urmaxl,urmax,mpicom,'urmax')
call shr_mpi_sum(uravgl,uravg,mpicom,'uravg')
call shr_mpi_sum(spavgl,spavg,mpicom,'spavg')
call shr_mpi_sum(urcntl,urcnt,mpicom,'urcnt')
if (my_task == 0 .and. urcnt > 0) then
uravg = uravg / urcnt
spavg = spavg / urcnt
write(logunit,*) trim(subname),' cart3d uravg,urmax,spavg = ',uravg,urmax,spavg
endif
endif
call mct_avect_clean(av3_s)
call mct_avect_clean(av3_d)
endif ! ku,kv
end subroutine seq_map_cart3d
!=======================================================================
subroutine seq_map_readdata(maprcfile, maprcname, mpicom, ID, &
ni_s, nj_s, av_s, gsmap_s, avfld_s, filefld_s, &
ni_d, nj_d, av_d, gsmap_d, avfld_d, filefld_d, string)
!--- lifted from work by J Edwards, April 2011
use shr_pio_mod, only : shr_pio_getiosys, shr_pio_getiotype
use pio, only : pio_openfile, pio_closefile, pio_read_darray, pio_inq_dimid, &
pio_inq_dimlen, pio_inq_varid, file_desc_t, io_desc_t, iosystem_desc_t, &
var_desc_t, pio_int, pio_get_var, pio_double, pio_initdecomp, pio_freedecomp
implicit none
!-----------------------------------------------------
!
! Arguments
!
character(len=*),intent(in) :: maprcfile
character(len=*),intent(in) :: maprcname
integer(IN) ,intent(in) :: mpicom
integer(IN) ,intent(in) :: ID
integer(IN) ,intent(out) ,optional :: ni_s
integer(IN) ,intent(out) ,optional :: nj_s
type(mct_avect) ,intent(inout),optional :: av_s
type(mct_gsmap) ,intent(in) ,optional :: gsmap_s
character(len=*),intent(in) ,optional :: avfld_s
character(len=*),intent(in) ,optional :: filefld_s
integer(IN) ,intent(out) ,optional :: ni_d
integer(IN) ,intent(out) ,optional :: nj_d
type(mct_avect) ,intent(inout),optional :: av_d
type(mct_gsmap) ,intent(in) ,optional :: gsmap_d
character(len=*),intent(in) ,optional :: avfld_d
character(len=*),intent(in) ,optional :: filefld_d
character(len=*),intent(in) ,optional :: string
!
! Local Variables
!
type(iosystem_desc_t), pointer :: pio_subsystem
integer(IN) :: pio_iotype
type(file_desc_t) :: File ! PIO file pointer
type(io_desc_t) :: iodesc ! PIO parallel io descriptor
integer(IN) :: rcode ! pio routine return code
type(var_desc_t) :: vid ! pio variable ID
integer(IN) :: did ! pio dimension ID
integer(IN) :: na ! size of source domain
integer(IN) :: nb ! size of destination domain
integer(IN) :: i ! index
integer(IN) :: mytask ! my task
integer(IN), pointer :: dof(:) ! DOF pointers for parallel read
character(len=256):: fileName
character(len=64) :: lfld_s, lfld_d, lfile_s, lfile_d
character(*),parameter :: areaAV_field = 'aream'
character(*),parameter :: areafile_s = 'area_a'
character(*),parameter :: areafile_d = 'area_b'
character(len=*),parameter :: subname = "(seq_map_readdata) "
!-----------------------------------------------------
if (seq_comm_iamroot(CPLID) .and. present(string)) then
write(logunit,'(A)') subname//' called for '//trim(string)
call shr_sys_flush(logunit)
endif
call MPI_COMM_RANK(mpicom,mytask,rcode)
lfld_s = trim(areaAV_field)
if (present(avfld_s)) then
lfld_s = trim(avfld_s)
endif
lfld_d = trim(areaAV_field)
if (present(avfld_d)) then
lfld_s = trim(avfld_d)
endif
lfile_s = trim(areafile_s)
if (present(filefld_s)) then
lfile_s = trim(filefld_s)
endif
lfile_d = trim(areafile_d)
if (present(filefld_d)) then
lfile_d = trim(filefld_d)
endif
call I90_allLoadF(trim(maprcfile),0,mpicom,rcode)
if(rcode /= 0) then
write(logunit,*)"Cant find maprcfile file ",trim(maprcfile)
call shr_sys_abort(trim(subname)//"i90_allLoadF File Not Found")
endif
call i90_label(trim(maprcname),rcode)
if(rcode /= 0) then
write(logunit,*)"Cant find label ",maprcname
call shr_sys_abort(trim(subname)//"i90_label Not Found")
endif
call i90_gtoken(filename,rcode)
if(rcode /= 0) then
write(logunit,*)"Error reading token ",filename
call shr_sys_abort(trim(subname)//"i90_gtoken Error on filename read")
endif
pio_subsystem => shr_pio_getiosys(ID)
pio_iotype = shr_pio_getiotype(ID)
rcode = pio_openfile(pio_subsystem, File, pio_iotype, filename)
if (present(ni_s)) then
rcode = pio_inq_dimid (File, 'ni_a', did) ! number of lons in input grid
rcode = pio_inq_dimlen(File, did , ni_s)
end if
if(present(nj_s)) then
rcode = pio_inq_dimid (File, 'nj_a', did) ! number of lats in input grid
rcode = pio_inq_dimlen(File, did , nj_s)
end if
if(present(ni_d)) then
rcode = pio_inq_dimid (File, 'ni_b', did) ! number of lons in output grid
rcode = pio_inq_dimlen(File, did , ni_d)
end if
if(present(nj_d)) then
rcode = pio_inq_dimid (File, 'nj_b', did) ! number of lats in output grid
rcode = pio_inq_dimlen(File, did , nj_d)
endif
!--- read and load area_a ---
if (present(av_s)) then
if (.not.present(gsmap_s)) then
call shr_sys_abort(trim(subname)//' ERROR av_s must have gsmap_s')
endif
rcode = pio_inq_dimid (File, 'n_a', did) ! size of input vector
rcode = pio_inq_dimlen(File, did , na)
i = mct_avect_indexra(av_s, trim(lfld_s))
call mct_gsmap_OrderedPoints(gsMap_s, mytask, dof)
call pio_initdecomp(pio_subsystem, pio_double, (/na/), dof, iodesc)
deallocate(dof)
rcode = pio_inq_varid(File,trim(lfile_s),vid)
call pio_read_darray(File, vid, iodesc, av_s%rattr(i,:), rcode)
call pio_freedecomp(File,iodesc)
end if
!--- read and load area_b ---
if (present(av_d)) then
if (.not.present(gsmap_d)) then
call shr_sys_abort(trim(subname)//' ERROR av_d must have gsmap_d')
endif
rcode = pio_inq_dimid (File, 'n_b', did) ! size of output vector
rcode = pio_inq_dimlen(File, did , nb)
i = mct_avect_indexra(av_d, trim(lfld_d))
call mct_gsmap_OrderedPoints(gsMap_d, mytask, dof)
call pio_initdecomp(pio_subsystem, pio_double, (/nb/), dof, iodesc)
deallocate(dof)
rcode = pio_inq_varid(File,trim(lfile_d),vid)
call pio_read_darray(File, vid, iodesc, av_d%rattr(i,:), rcode)
call pio_freedecomp(File,iodesc)
endif
call pio_closefile(File)
end subroutine seq_map_readdata
!=======================================================================
subroutine seq_map_avNormAvF(mapper, av_i, av_o, avf_i, avfifld, rList, norm)
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) , intent(inout) :: mapper ! mapper
type(mct_aVect) , intent(in) :: av_i ! input
type(mct_aVect) , intent(inout) :: av_o ! output
type(mct_aVect) , intent(in) :: avf_i ! extra src "weight"
character(len=*), intent(in) :: avfifld ! field name in avf_i
character(len=*), intent(in),optional :: rList ! fields list
logical , intent(in),optional :: norm ! normalize at end
!
integer(IN) :: lsize_i, lsize_f, kf, j
real(r8),allocatable :: frac_i(:)
logical :: lnorm
character(*),parameter :: subName = '(seq_map_avNormAvF) '
!-----------------------------------------------------
lnorm = .true.
if (present(norm)) then
lnorm = norm
endif
lsize_i = mct_aVect_lsize(av_i)
lsize_f = mct_aVect_lsize(avf_i)
if (lsize_i /= lsize_f) then
write(logunit,*) subname,' ERROR: lsize_i ne lsize_f ',lsize_i,lsize_f
call shr_sys_abort(subname//' ERROR size_i ne lsize_f')
endif
!--- extract frac_i field from avf_i to pass to seq_map_avNormArr ---
allocate(frac_i(lsize_i))
do j = 1,lsize_i
kf = mct_aVect_indexRA(avf_i,trim(avfifld))
frac_i(j) = avf_i%rAttr(kf,j)
enddo
if (present(rList)) then
call seq_map_avNormArr(mapper, av_i, av_o, frac_i, rList=rList, norm=lnorm)
else
call seq_map_avNormArr(mapper, av_i, av_o, frac_i, norm=lnorm)
endif
deallocate(frac_i)
end subroutine seq_map_avNormAvF
!=======================================================================
subroutine seq_map_avNormArr(mapper, av_i, av_o, norm_i, rList, norm)
implicit none
!-----------------------------------------------------
!
! Arguments
!
type(seq_map) , intent(inout) :: mapper! mapper
type(mct_aVect) , intent(in) :: av_i ! input
type(mct_aVect) , intent(inout) :: av_o ! output
real(r8) , intent(in), optional :: norm_i(:) ! source "weight"
character(len=*), intent(in), optional :: rList ! fields list
logical , intent(in), optional :: norm ! normalize at end
!
! Local variables
!
type(mct_aVect) :: avp_i , avp_o
integer(IN) :: j,kf
integer(IN) :: lsize_i,lsize_o
real(r8) :: normval
character(CX) :: lrList,appnd
logical :: lnorm
character(*),parameter :: subName = '(seq_map_avNormArr) '
character(len=*),parameter :: ffld = 'norm8wt' ! want something unique
!-----------------------------------------------------
lsize_i = mct_aVect_lsize(av_i)
lsize_o = mct_aVect_lsize(av_o)
lnorm = .true.
if (present(norm)) then
lnorm = norm
endif
if (present(norm_i)) then
if (.not.lnorm) call shr_sys_abort(subname//' ERROR norm_i and norm = false')
if (size(norm_i) /= lsize_i) call shr_sys_abort(subname//' ERROR size(norm_i) ne lsize_i')
endif
!--- create temporary avs for mapping ---
if (lnorm .or. present(norm_i)) then
appnd = ':'//ffld
else
appnd = ''
endif
if (present(rList)) then
call mct_aVect_init(avp_i, rList=trim( rList)//trim(appnd), lsize=lsize_i)
call mct_aVect_init(avp_o, rList=trim( rList)//trim(appnd), lsize=lsize_o)
else
lrList = mct_aVect_exportRList2c(av_i)
call mct_aVect_init(avp_i, rList=trim(lrList)//trim(appnd), lsize=lsize_i)
lrList = mct_aVect_exportRList2c(av_o)
call mct_aVect_init(avp_o, rList=trim(lrList)//trim(appnd), lsize=lsize_o)
endif
!--- copy av_i to avp_i and set ffld value to 1.0
!--- then multiply all fields by norm_i if norm_i exists
!--- this will do the right thing for the norm_i normalization
call mct_aVect_copy(aVin=av_i, aVout=avp_i, VECTOR=mct_usevector)
if (lnorm .or. present(norm_i)) then
kf = mct_aVect_indexRA(avp_i,ffld)
do j = 1,lsize_i
avp_i%rAttr(kf,j) = 1.0_r8
enddo
if (present(norm_i)) then
!$omp simd
do j = 1,lsize_i
avp_i%rAttr(:,j) = avp_i%rAttr(:,j)*norm_i(j)
enddo
endif
endif
!--- map ---
if (mapper%esmf_map) then
call shr_sys_abort(subname//' ERROR: esmf SMM not supported')
else
! MCT based SMM
call mct_sMat_avMult(avp_i, mapper%sMatp, avp_o, VECTOR=mct_usevector)
endif
!--- renormalize avp_o by mapped norm_i ---
if (lnorm) then
kf = mct_aVect_indexRA(avp_o,ffld)
!$omp simd
do j = 1,lsize_o
normval = avp_o%rAttr(kf,j)
if (normval /= 0.0_r8) then
normval = 1.0_r8/normval
endif
avp_o%rAttr(:,j) = avp_o%rAttr(:,j)*normval
enddo
endif
!--- copy back into av_o and we are done ---
call mct_aVect_copy(aVin=avp_o, aVout=av_o, VECTOR=mct_usevector)
call mct_aVect_clean(avp_i)
call mct_aVect_clean(avp_o)
end subroutine seq_map_avNormArr
end module seq_map_mod
Computing file changes ...
