Skip to content

Feature/wjiang/sphere pert #187

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
weiyuan-jiang wants to merge 5 commits into
base: develop
Choose a base branch
from
Draft

Feature/wjiang/sphere pert #187

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
2e77b4c
init check in for the sphere_random_fields
Jul 30, 2026
07259eb
just for test
Jul 30, 2026
fd1d159
use abstract class
Jul 30, 2026
fe04f03
merge the test version
Jul 30, 2026
57fd4cb
add sphere_random_fields to land_pert.F90
Jul 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 7 additions & 2 deletions GEOSlandpert_GridComp/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,21 @@ esma_set_this ()

set (SRCS
nr_ran2_gasdev.F90 nr_jacobi.F90 nr_fft.F90
random_fields.F90 land_pert.F90 force_and_cat_progn_pert_types.F90 LDAS_PertRoutines.F90
abstract_random_fields.F90 rectangle_random_fields.F90 sphere_random_fields.F90
land_pert.F90 force_and_cat_progn_pert_types.F90 LDAS_PertRoutines.F90
GEOS_LandPertGridComp.F90
)

esma_add_subdirectory(@smerfs)

esma_add_library (${this}
SRCS ${SRCS}
DEPENDENCIES GEOS_LdasShared GEOSensavg_GridComp GEOSland_GridComp MAPL ${MKL_LIBRARIES}
DEPENDENCIES GEOS_LdasShared GEOSensavg_GridComp GEOSland_GridComp MAPL smerfs_f ${MKL_LIBRARIES}
INCLUDES ${INC_ESMF} ${MKL_INCLUDE_DIRS})

# Enable MKL paths only when MKL was detected and the compiler supports it.
if (MKL_FOUND AND NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
target_compile_definitions(${this} PRIVATE MKL_AVAILABLE)
endif ()

esma_add_subdirectory(compare_pert)
82 changes: 70 additions & 12 deletions GEOSlandpert_GridComp/LDAS_PertRoutines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,8 @@ subroutine read_ens_prop_inputs( &
kw_xcorr_force_pert, &
kw_ycorr_force_pert, &
kw_tcorr_force_pert, &
kw_ccorr_force_pert &
kw_ccorr_force_pert, &
kw_use_sphere_pert &
)

! read ensemble propagation inputs from namelist file
Expand Down Expand Up @@ -243,6 +244,9 @@ subroutine read_ens_prop_inputs( &
type(force_pert_ccor_type), intent(out), optional :: &
kw_ccorr_force_pert

! Sphere-pert optional keyword output
logical, intent(out), optional :: kw_use_sphere_pert

! ------------------------

! locals
Expand Down Expand Up @@ -286,6 +290,9 @@ subroutine read_ens_prop_inputs( &

type(force_pert_ccor_type) :: ccorr_force_pert

! Sphere-pert local variable
logical :: use_sphere_pert

! Errlong variables
character(len=*), parameter :: Iam = 'read_ens_prop_inputs'

Expand Down Expand Up @@ -322,11 +329,15 @@ subroutine read_ens_prop_inputs( &
xcorr_force_pert, &
ycorr_force_pert, &
tcorr_force_pert, &
ccorr_force_pert
ccorr_force_pert, &
use_sphere_pert


root_proc = (myid ==0)

! Initialize sphere-pert default (safe: flat-FFT path is the default)
use_sphere_pert = .false.

! ---------------------------------------------------------------------
!
! initialize selected name list inputs
Expand Down Expand Up @@ -488,6 +499,9 @@ subroutine read_ens_prop_inputs( &

if (present(kw_ccorr_force_pert)) kw_ccorr_force_pert = ccorr_force_pert

! sphere-pert keyword
if (present(kw_use_sphere_pert)) kw_use_sphere_pert = use_sphere_pert

! ------------------------------------------------------------------
!
! save driver inputs into *ens_prop_inputs.nml file
Expand Down Expand Up @@ -804,13 +818,17 @@ subroutine get_progn_pert_param( pert_grid_l, N_progn_pert, &

integer, dimension(N_progn_pert_max) :: progn_pert_select

! Sphere-pert parameter for prognostic perturbations
logical :: use_sphere_pert

! ---------------------------------------------------------------

call get_progn_pert_inputs( pert_grid_l, &
descr_progn_pert, zeromean_progn_pert, coarsen_progn_pert, &
std_normal_max_progn_pert, std_progn_pert, &
xcorr_progn_pert, ycorr_progn_pert, &
tcorr_progn_pert, typ_progn_pert, ccorr_progn_pert )
tcorr_progn_pert, typ_progn_pert, ccorr_progn_pert, &
use_sphere_pert )

call get_pert_select( N_progn_pert_max, pert_grid_l, std_progn_pert, &
N_progn_pert, progn_pert_select )
Expand All @@ -827,7 +845,8 @@ subroutine get_progn_pert_param( pert_grid_l, N_progn_pert, &
std_normal_max_progn_pert, &
std_progn_pert, xcorr_progn_pert, ycorr_progn_pert, &
tcorr_progn_pert, typ_progn_pert, ccorr_progn_pert, &
progn_pert_select, progn_pert_param )
progn_pert_select, progn_pert_param, &
use_sphere_pert=use_sphere_pert )

end if

Expand Down Expand Up @@ -872,13 +891,17 @@ subroutine get_force_pert_param( pert_grid_l, N_force_pert, force_pert_param )

integer, dimension(N_force_pert_max) :: force_pert_select

! Shared sphere-pert parameter for forcing perturbations
logical :: use_sphere_pert

! ---------------------------------------------------------------

call get_force_pert_inputs( pert_grid_l, &
descr_force_pert, zeromean_force_pert, coarsen_force_pert, &
std_normal_max_force_pert, std_force_pert, &
xcorr_force_pert, ycorr_force_pert, &
tcorr_force_pert, typ_force_pert, ccorr_force_pert )
tcorr_force_pert, typ_force_pert, ccorr_force_pert, &
use_sphere_pert )

call get_pert_select( N_force_pert_max, pert_grid_l, std_force_pert, &
N_force_pert, force_pert_select )
Expand All @@ -894,7 +917,8 @@ subroutine get_force_pert_param( pert_grid_l, N_force_pert, force_pert_param )
std_normal_max_force_pert, &
std_force_pert, xcorr_force_pert, ycorr_force_pert, &
tcorr_force_pert, typ_force_pert, ccorr_force_pert, &
force_pert_select, force_pert_param )
force_pert_select, force_pert_param, &
use_sphere_pert=use_sphere_pert )

end if

Expand All @@ -906,7 +930,8 @@ subroutine get_force_pert_inputs( pert_grid_l, &
descr_force_pert, zeromean_force_pert, coarsen_force_pert, &
std_normal_max_force_pert, std_force_pert, &
xcorr_force_pert, ycorr_force_pert, &
tcorr_force_pert, typ_force_pert, ccorr_force_pert )
tcorr_force_pert, typ_force_pert, ccorr_force_pert, &
use_sphere_pert )

! get inputs for forcing perturbations for ALL forcing
! variables (including zero standard deviations) on a grid (typically,
Expand Down Expand Up @@ -959,7 +984,9 @@ subroutine get_force_pert_inputs( pert_grid_l, &
typ_force_pert

real, dimension(N_force_pert_max,N_force_pert_max), intent(out) :: &
ccorr_force_pert
ccorr_force_pert

logical, intent(out) :: use_sphere_pert

! ---------------------------------------------------------------
!
Expand Down Expand Up @@ -1004,6 +1031,13 @@ subroutine get_force_pert_inputs( pert_grid_l, &

root_proc = (myid==0)

! Read the shared sphere-pert option for forcing perturbations.
if (root_proc) then
call read_ens_prop_inputs(kw_echo=.false., &
kw_use_sphere_pert=use_sphere_pert)
end if
call MPI_Bcast(use_sphere_pert, 1, MPI_LOGICAL, 0, mpicomm, mpierr)

! ---------
!
! DESCR
Expand Down Expand Up @@ -1297,7 +1331,8 @@ subroutine get_progn_pert_inputs( pert_grid_l, &
descr_progn_pert, zeromean_progn_pert, coarsen_progn_pert, &
std_normal_max_progn_pert, std_progn_pert, &
xcorr_progn_pert, ycorr_progn_pert, &
tcorr_progn_pert, typ_progn_pert, ccorr_progn_pert )
tcorr_progn_pert, typ_progn_pert, ccorr_progn_pert, &
use_sphere_pert )

! get inputs for perturbations to prognostic variables for ALL prognostic
! variables (including zero standard deviations) on a grid (typically,
Expand Down Expand Up @@ -1348,7 +1383,9 @@ subroutine get_progn_pert_inputs( pert_grid_l, &
typ_progn_pert

real, dimension(N_progn_pert_max,N_progn_pert_max), intent(out) :: &
ccorr_progn_pert
ccorr_progn_pert

logical, intent(out) :: use_sphere_pert

! ---------------------------------------------------------------
!
Expand Down Expand Up @@ -1393,6 +1430,13 @@ subroutine get_progn_pert_inputs( pert_grid_l, &

root_proc = (myid==0)

! Read the shared sphere-pert option for prognostic perturbations.
if (root_proc) then
call read_ens_prop_inputs(kw_echo=.false., &
kw_use_sphere_pert=use_sphere_pert)
end if
call MPI_Bcast(use_sphere_pert, 1, MPI_LOGICAL, 0, mpicomm, mpierr)

! -------
!
! DESCR
Expand Down Expand Up @@ -1686,13 +1730,15 @@ subroutine assemble_pert_param( &
std_normal_max_pert, &
std_pert, xcorr_pert, ycorr_pert, tcorr_pert, typ_pert, &
ccorr_pert, &
pert_select, pert_param )
pert_select, pert_param, &
use_sphere_pert )

! assemble perturbation parameters for the prognostics or forcings
! that have been selected to be subject to perturbations
!
! reichle, 30 Nov 2001
! reichle, 27 May 2005 - redesign
! wjiang, 2026-07 - added optional sphere-pert arguments
!
! -----------------------------------------------------------------

Expand Down Expand Up @@ -1725,6 +1771,10 @@ subroutine assemble_pert_param( &

type(pert_param_type), dimension(:), pointer :: pert_param ! output

! Optional sphere-pert arguments:
! use_sphere_pert : .true. => use SMERFS for all selected fields in this call
logical, intent(in), optional :: use_sphere_pert

! --------------------------

! local variables
Expand All @@ -1738,8 +1788,14 @@ subroutine assemble_pert_param( &
character(len=400) :: err_msg
character(len=*), parameter :: Iam = 'assemble_pert_param'

logical :: do_sphere_pert

! -----------------------------------------------------------------

! Determine whether sphere-pert is requested for this call
do_sphere_pert = .false.
if (present(use_sphere_pert)) do_sphere_pert = use_sphere_pert

! extract info into pert_param of length N_pert
! (only those states listed in pert_select are affected by perturbations)

Expand All @@ -1765,8 +1821,10 @@ subroutine assemble_pert_param( &
pert_param(m)%tcorr = tcorr_pert(k)
pert_param(m)%typ = nint(typ_pert(k))

! Sphere-pert fields (safe defaults ensure flat-FFT when not requested)
pert_param(m)%use_sphere_pert = do_sphere_pert

n = 0
n = 0

do l=1,N_pert_max

Expand Down
68 changes: 68 additions & 0 deletions GEOSlandpert_GridComp/abstract_random_fields.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
module abstract_random_fieldsMod
use nr_ran2_gasdev, only: NRANDSEED
implicit none
private
public :: abstract_random_fields

type, abstract :: abstract_random_fields
character(len=:), allocatable :: ID_string
contains
procedure(generate_2d_random_field_interface), deferred :: generate_2d_Random_field
procedure(generate_white_field_interface), deferred :: generate_white_field
procedure(finalize_random_field_interface), deferred :: finalize
end type abstract_random_fields

abstract interface
subroutine generate_2d_random_field_interface(this, rseed, rfield, rfield2, lx, ly, dx, dy)
import :: abstract_random_fields, NRANDSEED
class(abstract_random_fields), intent(inout) :: this
integer, intent(inout) :: rseed(NRANDSEED)
real, intent(out) :: rfield(:,:), rfield2(:,:)
real, intent(in) :: lx, ly, dx, dy
end subroutine generate_2d_random_field_interface

subroutine generate_white_field_interface(this, rseed, rfield)
import :: abstract_random_fields, NRANDSEED
class(abstract_random_fields), intent(inout) :: this
integer, intent(inout) :: rseed(NRANDSEED)
real, target,intent(out) :: rfield(:,:)
end subroutine generate_white_field_interface

subroutine finalize_random_field_interface(this, rc)
import :: abstract_random_fields
class(abstract_random_fields), intent(inout) :: this
integer, optional, intent(out) :: rc
end subroutine finalize_random_field_interface
end interface
end module abstract_random_fieldsMod

module StringAbstractRandom_fieldsMapMod
use abstract_random_fieldsMod

#include "types/key_deferredLengthString.inc"
#define _value class (abstract_random_fields)
#define _value_allocatable
#define _value_equal_defined

! Work around Intel assignment handling for polymorphic map components.
#define _ASSIGN(dest,src) allocate(dest%key,source=src%key); if(allocated(src%value)) allocate(dest%value,source=src%value)
#define _MOVE(dest,src) call move_alloc(from=src%key,to=dest%key); if (allocated(src%value)) call move_alloc(from=src%value,to=dest%value)
#define _FREE(x) deallocate(x%key,x%value)

#define _map StringAbstractRandom_fieldsMap
#define _iterator StringAbstractRandom_fieldsMapIterator
#define _alt

#include "templates/map.inc"

#undef _alt
#undef _iterator
#undef _map
#undef _value
#undef _value_allocatable
#undef _key
#undef _value_equal_defined
#undef _ASSIGN
#undef _MOVE
#undef _FREE
end module StringAbstractRandom_fieldsMapMod
33 changes: 33 additions & 0 deletions GEOSlandpert_GridComp/compare_pert/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
# -----------------------------------------------------------------------
# compare_pert/CMakeLists.txt
#
# Builds the MPI comparison program for flat lat-lon FFT and sphere-native
# SMERFS perturbations on a 360x180 grid.
#
# Executable:
# compare_pert.x - runs both perturbation paths and writes binary output
#
# Python script (installed to ${CMAKE_INSTALL_PREFIX}/bin):
# compare_pert.py - reads both .bin outputs, produces comparison figure
#
# Both executables link against GEOSlandpert_GridComp, which transitively
# carries MAPL, ESMF, MKL, smerfs_f, MPI, netCDF, and HDF5.
# The MKL_AVAILABLE compile definition must be repeated here because the
# parent library applies it only as PRIVATE.
# -----------------------------------------------------------------------

# ---- Combined comparison executable ------------------------------------
ecbuild_add_executable(
TARGET compare_pert.x
SOURCES compare_pert.F90
LIBS GEOSlandpert_GridComp MAPL)

# ---- MKL compile definition (matches parent library treatment) ---------
if (MKL_FOUND AND NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
target_compile_definitions(compare_pert.x PRIVATE MKL_AVAILABLE)
endif ()

# ---- Install Python comparison script ----------------------------------
install(
PROGRAMS compare_pert.py
DESTINATION bin)
Loading
Loading