Skip to content

Commit

Permalink
Merge pull request #2546 from GEOS-ESM/feature/ygyu/geostation
Browse files Browse the repository at this point in the history 
updates for geostationary output part 2
  • Loading branch information
@tclune
tclune authored Jan 24, 2024
2 parents 7baf233 + 4f53352 commit a16827d
Show file tree
Hide file tree
Showing 9 changed files with 465 additions and 61 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## [Unreleased]

### Added
- Convert from ABI Fixed Grid to lon/lat coordinates used in MAPL_XYGridFactory (supporting geostationary GOES-R series)
- Modify trajectory sampler for a collection with multiple platforms: P3B (air craft) + FIREX
- Modify swath sampler to handle two Epoch swath grids
- Handle regrid accumulate for time step (1 sec) during which no obs exists
Expand Down
50 changes: 50 additions & 0 deletions base/MAPL_Comms.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,8 @@ module MAPL_CommsMod

interface MAPL_BcastShared
module procedure MAPL_BcastShared_1DR4
module procedure MAPL_BcastShared_1DR8
module procedure MAPL_BcastShared_2DI4
module procedure MAPL_BcastShared_2DR4
module procedure MAPL_BcastShared_2DR8
end interface
Expand Down Expand Up @@ -1086,6 +1088,30 @@ subroutine MAPL_BcastShared_1DR4(VM, Data, N, Root, RootOnly, rc)

end subroutine MAPL_BcastShared_1DR4

subroutine MAPL_BcastShared_1DR8(VM, Data, N, Root, RootOnly, rc)
type(ESMF_VM) :: VM
real(kind=REAL64), pointer, intent(INOUT) :: Data(:)
integer, intent(IN ) :: N
integer, optional, intent(IN ) :: Root
logical, intent(IN ) :: RootOnly
integer, optional, intent( OUT) :: rc
integer :: status

if(.not.MAPL_ShmInitialized) then
if (RootOnly) then
_RETURN(ESMF_SUCCESS)
end if
call MAPL_CommsBcast(vm, DATA=Data, N=N, ROOT=Root, _RC)
else
call MAPL_SyncSharedMemory(_RC)
call MAPL_BroadcastToNodes(Data, N=N, ROOT=Root, _RC)
call MAPL_SyncSharedMemory(_RC)
endif

_RETURN(ESMF_SUCCESS)

end subroutine MAPL_BcastShared_1DR8

subroutine MAPL_BcastShared_2DR4(VM, Data, N, Root, RootOnly, rc)
type(ESMF_VM) :: VM
real, pointer, intent(INOUT) :: Data(:,:)
Expand Down Expand Up @@ -1142,6 +1168,30 @@ subroutine MAPL_BcastShared_2DR8(VM, Data, N, Root, RootOnly, rc)
_RETURN(ESMF_SUCCESS)

end subroutine MAPL_BcastShared_2DR8

subroutine MAPL_BcastShared_2DI4(VM, Data, N, Root, RootOnly, rc)
type(ESMF_VM) :: VM
integer, pointer, intent(INOUT) :: Data(:,:)
integer, intent(IN ) :: N
integer, optional, intent(IN ) :: Root
logical, intent(IN ) :: RootOnly
integer, optional, intent( OUT) :: rc
integer :: status

if(.not.MAPL_ShmInitialized) then
if (RootOnly) then
_RETURN(ESMF_SUCCESS)
end if
call MAPL_CommsBcast(vm, DATA=Data, N=N, ROOT=Root, _RC)
else
call MAPL_SyncSharedMemory(_RC)
call MAPL_BroadcastToNodes(Data, N=N, ROOT=Root, _RC)
call MAPL_SyncSharedMemory(_RC)
endif

_RETURN(ESMF_SUCCESS)

end subroutine MAPL_BcastShared_2DI4

! Rank 0
!---------------------------
Expand Down
25 changes: 9 additions & 16 deletions base/MAPL_EsmfRegridder.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1431,6 +1431,7 @@ subroutine create_route_handle(this, kind, rc)

integer :: srcTermProcessing
integer, pointer :: factorIndexList(:,:)
integer, allocatable :: dstMaskValues(:)
real(ESMF_KIND_R8), pointer :: factorList(:)
type(ESMF_RouteHandle) :: dummy_rh
type(ESMF_UnmappedAction_Flag) :: unmappedaction
Expand Down Expand Up @@ -1494,24 +1495,16 @@ subroutine create_route_handle(this, kind, rc)
call ESMF_AttributeGet(spec%grid_in, name='Global',value=global,rc=status)
if (.not.global) unmappedaction=ESMF_UNMAPPEDACTION_IGNORE
end if
if (has_mask) dstMaskValues = [MAPL_MASK_OUT] ! otherwise unallocated
select case (spec%regrid_method)
case (REGRID_METHOD_BILINEAR, REGRID_METHOD_BILINEAR_MONOTONIC)
if (has_mask) then
call ESMF_FieldRegridStore(src_field, dst_field, &
& dstMaskValues = [MAPL_MASK_OUT], &
& regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
& linetype=ESMF_LINETYPE_GREAT_CIRCLE, & ! closer to SJ Lin interpolation weights?
& srcTermProcessing = srcTermProcessing, &
& factorList=factorList, factorIndexList=factorIndexList, &
& routehandle=route_handle, unmappedaction=unmappedaction, _RC)
else
call ESMF_FieldRegridStore(src_field, dst_field, &
& regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
& linetype=ESMF_LINETYPE_GREAT_CIRCLE, & ! closer to SJ Lin interpolation weights?
& srcTermProcessing = srcTermProcessing, &
& factorList=factorList, factorIndexList=factorIndexList, &
& routehandle=route_handle, unmappedaction=unmappedaction, _RC)
end if
call ESMF_FieldRegridStore(src_field, dst_field, &
& dstMaskValues = dstMaskValues, &
& regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
& linetype=ESMF_LINETYPE_GREAT_CIRCLE, & ! closer to SJ Lin interpolation weights?
& srcTermProcessing = srcTermProcessing, &
& factorList=factorList, factorIndexList=factorIndexList, &
& routehandle=route_handle, unmappedaction=unmappedaction, _RC)
case (REGRID_METHOD_PATCH)

_ASSERT(.not.has_mask, "destination masking with this regrid type is unsupported")
Expand Down
125 changes: 122 additions & 3 deletions base/MAPL_ObsUtil.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,21 @@ module MAPL_ObsUtilMod
use ESMF
use Plain_netCDF_Time
use netCDF
use MAPL_BaseMod, only: MAPL_UNDEF
use MAPL_CommsMod, only : MAPL_AM_I_ROOT
use pFIO_FileMetadataMod, only : FileMetadata
use pFIO_NetCDF4_FileFormatterMod, only : NetCDF4_FileFormatter
use, intrinsic :: iso_fortran_env, only: REAL32, REAL64
implicit none
integer, parameter :: mx_ngeoval = 60
!! private
! GRS80 by Moritz
real(REAL64) :: r_eq=6378137.d0
real(REAL64) :: r_pol=6356752.31414d0
real(REAL64) :: H_sat=42164160.d0
! GOES-R
real(REAL64) :: lambda0_SatE=-1.308996939d0 ! -75 deg Satellite East
real(REAL64) :: lambda0_SatW=-2.39110107523d0 ! -137 deg Satellite West
real(REAL64) :: lambda0_SatT=-1.56206968053d0 ! -89.5 deg Satellite Test

public :: obs_unit
type :: obs_unit
Expand Down Expand Up @@ -293,7 +301,7 @@ subroutine read_M_files_4_swath ( filenames, Xdim, Ydim, &
character(len=ESMF_MAXSTR), optional, intent(in) :: var_name_time
real(ESMF_KIND_R8), allocatable, optional, intent(inout) :: lon(:,:)
real(ESMF_KIND_R8), allocatable, optional, intent(inout) :: lat(:,:)
real(ESMF_KIND_R8), allocatable, optional, intent(inout) :: time(:,:)
real(ESMF_KIND_R8), allocatable, optional, intent(inout) :: time(:,:)
logical, optional, intent(in) :: Tfilter
integer, optional, intent(out) :: rc

Expand Down Expand Up @@ -439,7 +447,7 @@ subroutine read_M_files_4_swath ( filenames, Xdim, Ydim, &
deallocate(lat)
allocate (lat(Xdim, Ydim))
end if

jx=0
do i = 1, M
filename = filenames(i)
Expand Down Expand Up @@ -723,4 +731,115 @@ function union_platform(a, b, rc)
end function union_platform


! From GOES-R SERIES PRODUCT DEFINITION AND USERS’ GUIDE
!
subroutine ABI_XY_2_lonlat (x, y, lambda0, lon, lat, mask)
implicit none
real(REAL64), intent(in) :: x, y
real(REAL64), intent(in) :: lambda0
real(REAL64), intent(out):: lon, lat
integer, optional, intent(out):: mask
real(REAL64) :: a0, b0, c0, rs, Sx, Sy, Sz, t
real(REAL64) :: a, b, H
real(REAL64) :: delta

a=r_eq; b=r_pol; H=H_sat

if (present(mask)) mask=0
a0 = sin(x)*sin(x) + cos(x)*cos(x)*( cos(y)*cos(y) + (a/b)*(a/b)*sin(y)*sin(y) )
b0 = -2.d0 * H * cos(x) * cos(y)
c0 = H*H - a*a
delta = b0*b0 - 4.d0*a0*c0
if (delta < 0.d0) then
lon = MAPL_UNDEF
lat = MAPL_UNDEF
if (present(mask)) mask=0
return
end if
rs = ( -b0 - sqrt(b0*b0 - 4.d0*a0*c0) ) / (2.d0*a0)
Sx = rs * cos(x) * cos(y)
Sy = -rs * sin(x)
Sz = rs * cos(x) * sin(y)
lon = lambda0 - atan (Sy/(H - Sx))
lat = atan ( (a/b)**2.d0 * Sz / sqrt ((H -Sx)**2.d0 + Sy*Sy) )

t = H*(H-Sx) - ( Sy*Sy + (a/b)**2.d0 *Sz*Sz )
if (t < 0) then
lon = MAPL_UNDEF
lat = MAPL_UNDEF
if (present(mask)) mask=0
else
if (present(mask)) mask=1
end if

end subroutine ABI_XY_2_lonlat


subroutine lonlat_2_ABI_XY (lon, lat, lambda0, x, y, mask)
implicit none
real(REAL64), intent(in) :: lon, lat
real(REAL64), intent(in) :: lambda0
real(REAL64), intent(out):: x, y
integer, intent(out):: mask
real(REAL64) :: theta_c
real(REAL64) :: e2, rc, Sx, Sy, Sz, t
real(REAL64) :: a, b, H
real*8 :: delta

a=r_eq; b=r_pol; H=H_sat

theta_c = atan( (b/a)**2.d0 * tan(lat) )
e2 = 1.d0 - (b/a)**2.d0 ! (a^2-b^2)/a^2
rc = b / sqrt( 1.d0 - e2 * cos(theta_c)**2.d0 )
Sx = H - rc * cos(theta_c) * cos( lon - lambda0 )
Sy = - rc * cos(theta_c) * sin( lon - lambda0 )
Sz = rc * sin(theta_c)
x = - asin ( Sy / sqrt (Sx*Sx + Sy*Sy + Sz*Sz) )
y = atan ( Sz / Sx )

t = H*(H-Sx) - ( Sy*Sy + (a/b)**2.d0 *Sz*Sz )
if (t < 0) then
mask = 1
else
mask = 0
end if

end subroutine lonlat_2_ABI_XY


subroutine test_conversion
implicit none
real*8 :: x0
real*8 :: y0
real*8 :: lam, the
real*8 :: lon, lat
integer :: mask
real*8 :: xnew, ynew

! two points mapping: (x0, y0) <--> (lam, the)
x0 = -0.024052d0
y0 = 0.095340d0
lam = -1.478135612d0
the = 0.590726971d0

call ABI_XY_2_lonlat (x0, y0, lambda0_SatE, lon, lat, mask)
write(6, 111) 'x,y 2 ll'
write(6, 111) 'x,y=', x0, y0
write(6, 111) 'lon,lat=', lon, lat
write(6, 121) 'mask=', mask
write(6, 111) 'errror lon,lat=', lon - lam, lat-the

call lonlat_2_ABI_XY (lam, the, lambda0_SatE, xnew, ynew, mask)
write(6, 111) 'll 2 xy'
write(6, 111) 'lon,lat=', lam, the
write(6, 111) 'x,y=', xnew, ynew
write(6, 121) 'mask=', mask
write(6, 111) 'errror lon,lat=', xnew -x0, ynew-y0

101 format (2x, a,10(2x,f15.8))
111 format (2x, a,20(2x,f25.11))
121 format (2x, a,10(2x,i8))

end subroutine test_conversion

end module MAPL_ObsUtilMod
7 changes: 5 additions & 2 deletions base/MAPL_SwathGridFactory.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,7 @@ function create_basic_grid(this, unusable, rc) result(grid)
if (this%lm /= MAPL_UNDEFINED_INTEGER) then
call ESMF_AttributeSet(grid, name='GRID_LM', value=this%lm, _RC)
end if
call ESMF_AttributeSet(grid, 'GridType', 'LatLon', _RC)
call ESMF_AttributeSet(grid, 'GridType', 'Swath', _RC)
call ESMF_AttributeSet(grid, 'Global', .false., _RC)

_RETURN(_SUCCESS)
Expand Down Expand Up @@ -246,7 +246,7 @@ subroutine add_horz_coordinates_from_file(this, grid, unusable, rc)

_UNUSED_DUMMY(unusable)

!! call ESMF_VMGetCurrent(vm,_RC)
call ESMF_VMGetCurrent(vm,_RC)
!! call ESMF_VMGet(vm, mpiCommunicator=mpic, localPet=mypet, petCount=petCount, _RC)

Xdim=this%im_world
Expand Down Expand Up @@ -480,6 +480,9 @@ subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc
call ESMF_ConfigGetAttribute(config, this%epoch, label=prefix//'Epoch:', default=300, _RC)
call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'Epoch_init:', default='2006', _RC)

write(6,'(2x,a,100i10)') 'nail 2, nx,ny,im,jm,lm',&
this%nx,this%ny,this%im_world,this%jm_world,this%lm

call lgr%debug(' %a %a', 'CurrTime =', trim(tmp))

if ( index(tmp, 'T') /= 0 .OR. index(tmp, '-') /= 0 ) then
Expand Down
Loading

0 comments on commit a16827d

Please sign in to comment.