Skip to content
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

Inline harmonic analysis #744

Merged
merged 3 commits into from
Dec 9, 2024
Merged

Inline harmonic analysis #744

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
83f44b2
Inline harmonic analysis
c2xu Oct 23, 2024
b101d1c
Inline harmonic analysis
c2xu Oct 30, 2024
1c32e40
Inline harmonic analysis
c2xu Nov 27, 2024
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
103 changes: 55 additions & 48 deletions src/diagnostics/MOM_harmonic_analysis.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,7 @@ end subroutine HA_register

!> This subroutine accumulates the temporal basis functions in FtF.
!! The tidal constituents are those used in MOM_tidal_forcing, plus the mean (of zero frequency).
!! Only the main diagonal and entries below it are calculated, which are needed for Cholesky decomposition.
subroutine HA_accum_FtF(Time, CS)
type(time_type), intent(in) :: Time !< The current model time
type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module
Expand All @@ -191,27 +192,31 @@ subroutine HA_accum_FtF(Time, CS)
nc = CS%nc
now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref)

! Accumulate FtF
CS%FtF(1,1) = CS%FtF(1,1) + 1.0 !< For the zero frequency
!< First entry, corresponding to the zero frequency constituent (mean)
CS%FtF(1,1) = CS%FtF(1,1) + 1.0

do c=1,nc
icos = 2*c
isin = 2*c+1
cosomegat = cos(CS%freq(c) * now + CS%phase0(c))
sinomegat = sin(CS%freq(c) * now + CS%phase0(c))

! First column, corresponding to the zero frequency constituent (mean)
CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat
CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat
CS%FtF(1,icos) = CS%FtF(icos,1)
CS%FtF(1,isin) = CS%FtF(isin,1)
do cc=c,nc

do cc=1,c
iccos = 2*cc
issin = 2*cc+1
ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc))
ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc))

! Interior of the matrix, corresponding to the products of cosine and sine terms
CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat
CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat
CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat
CS%FtF(isin,issin) = CS%FtF(isin,issin) + sinomegat * ssinomegat
enddo ! cc=c,nc
enddo ! cc=1,c
enddo ! c=1,nc

end subroutine HA_accum_FtF
Expand Down Expand Up @@ -276,14 +281,18 @@ subroutine HA_accum_FtSSH(key, data, Time, G, CS)

is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je

! Accumulate FtF and FtSSH
!< First entry, corresponding to the zero frequency constituent (mean)
do j=js,je ; do i=is,ie
ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j))
enddo ; enddo

!< The remaining entries
do c=1,nc
icos = 2*c
isin = 2*c+1
cosomegat = cos(CS%freq(c) * now + CS%phase0(c))
sinomegat = sin(CS%freq(c) * now + CS%phase0(c))
do j=js,je ; do i=is,ie
ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j))
ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat
ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat
enddo ; enddo
Expand Down Expand Up @@ -315,7 +324,7 @@ subroutine HA_write(ha1, Time, G, CS)
! Local variables
real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A]
integer :: year, month, day, hour, minute, second
integer :: nc, k, is, ie, js, je
integer :: nc, i, j, k, is, ie, js, je

character(len=255) :: filename !< Output file name
type(MOM_infra_file) :: cdf !< The file handle for output harmonic constants
Expand Down Expand Up @@ -348,6 +357,11 @@ subroutine HA_write(ha1, Time, G, CS)
call create_MOM_file(cdf, trim(filename), cdf_vars, &
2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G)

! Add the initial field back to the mean state
do j=js,je ; do i=is,ie
FtSSHw(i,j,1) = FtSSHw(i,j,1) + ha1%ref(i,j)
enddo ; enddo

! Write data
call MOM_write_field(cdf, cdf_fields(1), G%domain, FtSSHw(:,:,1), 0.0)
do k=1,nc
Expand All @@ -362,75 +376,68 @@ subroutine HA_write(ha1, Time, G, CS)

end subroutine HA_write

!> This subroutine computes the harmonic constants (stored in FtSSHw) using the dot products of the temporal
!> This subroutine computes the harmonic constants (stored in x) using the dot products of the temporal
!! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis
!! functions accumulated in FtSSH. The system is solved by Cholesky decomposition,
!!
!! FtF * FtSSHw = FtSSH, => FtFw * (FtFw' * FtSSHw) = FtSSH,
!! FtF * x = FtSSH, => L * (L' * x) = FtSSH, => L * y = FtSSH,
!!
!! where FtFw is a lower triangular matrix, and the prime denotes matrix transpose.
!! where L is the lower triangular matrix, y = L' * x, and x is the solution vector.
!!
subroutine HA_solver(ha1, nc, FtF, FtSSHw)
subroutine HA_solver(ha1, nc, FtF, x)
type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field
integer, intent(in) :: nc !< Number of harmonic constituents
real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim]
real, dimension(:,:,:), allocatable, intent(out) :: FtSSHw !< Work array for Cholesky decomposition [A]
real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1), &
intent(out) :: x !< Solution vector of harmonic constants [A]

! Local variables
real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim]
real, dimension(:), allocatable :: tmp1 !< Temporary variable for Cholesky decomposition [nondim]
real, dimension(:,:), allocatable :: tmp2 !< Temporary variable for Cholesky decomposition [A]
real, dimension(:,:), allocatable :: FtFw !< Lower triangular matrix for Cholesky decomposition [nondim]
integer :: k, m, n, is, ie, js, je

is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je

allocate(tmp1(1:2*nc+1), source=0.0)
allocate(tmp2(is:ie,js:je), source=0.0)
allocate(FtFw(1:2*nc+1,1:2*nc+1), source=0.0)
allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0)

! Construct FtFw
FtFw(:,:) = 0.0
real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim]
real, dimension(2*nc+1,2*nc+1) :: L !< Lower triangular matrix of Cholesky decomposition [nondim]
real, dimension(2*nc+1) :: tmp1 !< Inverse of the diagonal entries of L [nondim]
real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je) :: tmp2 !< 2D temporary array involving FtSSH [A]
real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1) :: y !< 3D temporary array, i.e., L' * x [A]
integer :: k, m, n

! Cholesky decomposition
do m=1,2*nc+1

! First, calculate the diagonal entries
tmp0 = 0.0
do k=1,m-1
tmp0 = tmp0 + FtFw(m,k) * FtFw(m,k)
do k=1,m-1 ! This loop operates along the m-th row
tmp0 = tmp0 + L(m,k) * L(m,k)
enddo
FtFw(m,m) = sqrt(FtF(m,m) - tmp0)
tmp1(m) = 1 / FtFw(m,m)
do k=m+1,2*nc+1
L(m,m) = sqrt(FtF(m,m) - tmp0) ! This is the m-th diagonal entry

! Now calculate the off-diagonal entries
tmp1(m) = 1 / L(m,m)
do k=m+1,2*nc+1 ! This loop operates along the column below the m-th diagonal entry
tmp0 = 0.0
do n=1,m-1
tmp0 = tmp0 + FtFw(k,n) * FtFw(m,n)
tmp0 = tmp0 + L(k,n) * L(m,n)
enddo
FtFw(k,m) = (FtF(k,m) - tmp0) * tmp1(m)
L(k,m) = (FtF(k,m) - tmp0) * tmp1(m) ! This is the k-th off-diagonal entry below the m-th diagonal entry
enddo
enddo

! Solve for (FtFw' * FtSSHw)
FtSSHw(:,:,:) = ha1%FtSSH(:,:,:)
! Solve for y from L * y = FtSSH
do k=1,2*nc+1
tmp2(:,:) = 0.0
do m=1,k-1
tmp2(:,:) = tmp2(:,:) + FtFw(k,m) * FtSSHw(:,:,m)
tmp2(:,:) = tmp2(:,:) + L(k,m) * y(:,:,m)
Hallberg-NOAA marked this conversation as resolved.
Show resolved Hide resolved
enddo
FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k)
y(:,:,k) = (ha1%FtSSH(:,:,k) - tmp2(:,:)) * tmp1(k)
enddo

! Solve for FtSSHw
! Solve for x from L' * x = y
do k=2*nc+1,1,-1
tmp2(:,:) = 0.0
do m=k+1,2*nc+1
tmp2(:,:) = tmp2(:,:) + FtSSHw(:,:,m) * FtFw(m,k)
tmp2(:,:) = tmp2(:,:) + L(m,k) * x(:,:,m)
enddo
FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k)
x(:,:,k) = (y(:,:,k) - tmp2(:,:)) * tmp1(k)
enddo

deallocate(tmp1)
deallocate(tmp2)
deallocate(FtFw)

end subroutine HA_solver

!> \namespace harmonic_analysis
Expand All @@ -441,7 +448,7 @@ end subroutine HA_solver
!! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent
!! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares,
!!
!! (F' * F) * x = F' * SSH_in,
!! (F' * F) * x = F' * SSH_in, => FtF * x = FtSSH,
!!
!! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by
!! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is
Expand Down
Loading