Skip to content

Modify gsibec workarounds to fix bug with linear variable change #527

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
SamuelDegelia-NOAA wants to merge 7 commits into
base: develop
Choose a base branch
from
+291 −5
Draft

Modify gsibec workarounds to fix bug with linear variable change #527

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fa66d58
add new changes
Jan 20, 2026
fd70859
Merge branch 'NOAA-EMC:develop' into feature/update_gsibec_lvc
SamuelDegelia-NOAA Jan 28, 2026
c58f82a
Merge branch 'NOAA-EMC:develop' into feature/update_gsibec_lvc
SamuelDegelia-NOAA Jan 28, 2026
c1a8afd
add back replacement codes for guess_grids.f90
Jan 29, 2026
6b67dc6
Merge branch 'feature/update_gsibec_lvc' of https://github.com/Samuel…
Jan 29, 2026
3ea1d69
Merge branch 'develop' into feature/update_gsibec_lvc
SamuelDegelia-NOAA Jan 29, 2026
2785741
fix white space
Jan 29, 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
6 changes: 5 additions & 1 deletion sorc/_workaround_/gsibec/compute_qvar3d.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ subroutine compute_qvar3d
use berror, only: dssv
use derivsmod, only: qsatg,qgues
use control_vectors, only: cvars3d
use gridmod, only: lat2,lon2,nsig
use gridmod, only: lat2,lon2,nsig,regional
use constants, only: zero,one,fv,r100,qmin
use guess_grids, only: fact_tv,ntguessig,nfldsig,ges_tsen,ges_prsl,ges_qsat
use mpeu_util, only: getindex
Expand Down Expand Up @@ -78,6 +78,7 @@ subroutine compute_qvar3d
real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL()
integer(i_kind):: maxvarq1

real(r_kind), parameter :: rmiss_th = -1.0e30

nrf3_q=getindex(cvars3d,'q')
nrf3_cw=getindex(cvars3d,'cw')
Expand Down Expand Up @@ -138,6 +139,9 @@ subroutine compute_qvar3d
do j=1,lon2
do i=1,lat2
rhgues(i,j,k)=qgues(i,j,k)/qsatg(i,j,k)
if(regional .and. ges_tsen(i,j,k,ntguessig) < rmiss_th) then
rhgues(i,j,k)=0.5
endif
end do
end do
end do
Expand Down
130 changes: 130 additions & 0 deletions sorc/_workaround_/gsibec/gsi_convert_cv_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
module gsi_convert_cv_mod
use m_kinds, only: r_kind
use constants, only: epsilon=>fv
use constants, only: zero,one
private

public :: gsi_t_to_tv_tl
public :: gsi_t_to_tv_ad
public :: gsi_tv_to_t_tl
public :: gsi_tv_to_t_ad

interface gsi_t_to_tv_tl
module procedure t_to_tv_tl_
end interface
interface gsi_t_to_tv_ad
module procedure t_to_tv_ad_
end interface
interface gsi_tv_to_t_tl
module procedure tv_to_t_tl_
end interface
interface gsi_tv_to_t_ad
module procedure tv_to_t_ad_
end interface

contains

subroutine t_to_tv_tl_(t,t_tl,q,q_tl,tv_tl)

implicit none
real(r_kind), intent(in ) :: t(:,:,:)
real(r_kind), intent(in ) :: t_tl(:,:,:)
real(r_kind), intent(in ) :: q(:,:,:)
real(r_kind), intent(in ) :: q_tl(:,:,:)
real(r_kind), intent(out) :: tv_tl(:,:,:)

tv_tl = t_tl*(one + epsilon*q) + t*epsilon*q_tl

end subroutine t_to_tv_tl_

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

subroutine t_to_tv_ad_(t,t_ad,q,q_ad,tv_ad)

implicit none
real(r_kind), intent(in ) :: t(:,:,:)
real(r_kind), intent(inout) :: t_ad(:,:,:)
real(r_kind), intent(in ) :: q(:,:,:)
real(r_kind), intent(inout) :: q_ad(:,:,:)
real(r_kind), intent(inout) :: tv_ad(:,:,:)

t_ad = t_ad + tv_ad * (one + epsilon*q)
q_ad = q_ad + tv_ad * epsilon*t
tv_ad= zero

end subroutine t_to_tv_ad_

subroutine tv_to_t_tl_(tv,tv_tl,q,q_tl,t_tl,t)

implicit none
real(r_kind), intent(in ) :: tv(:,:,:)
real(r_kind), intent(in ) :: tv_tl(:,:,:)
real(r_kind), intent(in ) :: q(:,:,:)
real(r_kind), intent(in ) :: q_tl(:,:,:)
real(r_kind), intent(inout) :: t_tl(:,:,:)
real(r_kind), intent(in ) :: t(:,:,:)

integer :: i,j,k
real(r_kind), parameter :: rmiss_th = -1.0e30

do k = 1, size(t_tl,3)
do j = 1, size(t_tl,2)
do i = 1, size(t_tl,1)

if(t(i,j,k) < rmiss_th) then
t_tl(i,j,k) = zero
cycle
endif

t_tl(i,j,k)= (tv_tl(i,j,k)*(one+epsilon*q(i,j,k))-tv(i,j,k)*epsilon*q_tl(i,j,k))/(one+epsilon*q(i,j,k))**2

enddo
enddo
enddo

end subroutine tv_to_t_tl_

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

subroutine tv_to_t_ad_(tv,tv_ad,q,q_ad,t_ad,t)

implicit none
real(r_kind), intent(in ) :: tv(:,:,:)
real(r_kind), intent(inout) :: tv_ad(:,:,:)
real(r_kind), intent(in ) :: q(:,:,:)
real(r_kind), intent(inout) :: q_ad(:,:,:)
real(r_kind), intent(inout) :: t_ad(:,:,:)
real(r_kind), intent(in ) :: t(:,:,:)

real(r_kind),allocatable :: temp(:,:,:)

integer :: i,j,k
real(r_kind), parameter :: rmiss_th = -1.0e30

allocate(temp(size(t_ad,1),size(t_ad,2),size(t_ad,3)))

do k = 1, size(t_ad,3)
do j = 1, size(t_ad,2)
do i = 1, size(t_ad,1)

if(t(i,j,k) < rmiss_th) then
tv_ad(i,j,k) = zero
q_ad(i,j,k) = zero
t_ad(i,j,k) = zero
cycle
endif

temp(i,j,k) = t_ad(i,j,k)/(epsilon*q(i,j,k)+one)

tv_ad(i,j,k) = tv_ad(i,j,k) + temp(i,j,k)
q_ad(i,j,k) = q_ad(i,j,k) - tv(i,j,k)*epsilon*temp(i,j,k)/(epsilon*q(i,j,k)+one)
t_ad(i,j,k) = zero
enddo
enddo
enddo

deallocate(temp)

end subroutine tv_to_t_ad_

end module gsi_convert_cv_mod
148 changes: 148 additions & 0 deletions sorc/_workaround_/gsibec/normal_rh_to_q.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
subroutine normal_rh_to_q(rhnorm,t,p,q)
!$$$ subprogram documentation block
! . . . .
! subprogram: normal_rh_to_q tlm for normalized RH to q
! prgmmr: wu org: np20 date: 2005-03-06
!
! abstract: get specific humidity q from normalized RH
!
! program history log:
! 2005-03-06 wu
! 2005-03-30 treadon - reformat code (cosmetic change only)
! 2005-11-21 kleist - use 3d pressure increment for coupling
! 2005-11-21 derber modify to make qoption =1 work same as =2
! 2006-01-09 derber move sigsum calculation to compute_derived and clean up
! 2006-07-31 kleist - analysis variable changed from ln(ps) to ps
! 2008-05-28 safford - rm unused uses
! 2015-10-27 mahajan - code clean-up
!
! input argument list:
! rhnorm - normalized RH
! t - virtual temperature
! p - psfc
!
! output argument list:
! q - specific humidity
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use m_kinds, only: r_kind,i_kind
use derivsmod, only: dqdrh,dqdp,dqdt
use jfunc, only: qoption
use gridmod, only: lat2,lon2,nsig,regional
use guess_grids, only: ges_tsen,ntguessig
use constants, only: zero

implicit none

real(r_kind),intent(in ) :: rhnorm(lat2,lon2,nsig)
real(r_kind),intent(in ) :: t(lat2,lon2,nsig)
real(r_kind),intent(in ) :: p(lat2,lon2,nsig+1)
real(r_kind),intent( out) :: q(lat2,lon2,nsig)

real(r_kind), parameter :: rmiss_th = -1.0e30
integer(i_kind) i,j,k

! Convert normalized rh to q
do k=1,nsig
do j=1,lon2
do i=1,lat2
if(regional .and. ges_tsen(i,j,k,ntguessig) < rmiss_th) then
q(i,j,k) = zero
cycle
endif
q(i,j,k) = dqdrh(i,j,k)*rhnorm(i,j,k)
if ( qoption == 2 ) then
q(i,j,k) = q(i,j,k) + &
dqdt(i,j,k)*t(i,j,k) - &
dqdp(i,j,k)*(p(i,j,k) + p(i,j,k+1))
endif
enddo
enddo
enddo

return

end subroutine normal_rh_to_q

subroutine normal_rh_to_q_ad(rhnorm,t,p,q)
!$$$ subprogram documentation block
! . . . .
! subprogram: normal_rh_to_q_ad adjoint of normal_rh_to_q
! prgmmr: wu org: np20 date: 2005-03-06
!
! abstract: adjoint of normal_rh_to_q
!
! program history log:
! 2005-03-06 wu
! 2005-03-30 treadon - reformat code (cosmetic change only)
! 2005-11-21 kleist - use 3d pressure increment for coupling
! 2005-11-21 derber modify to make qoption =1 work same as =2
! 2006-01-09 derber move sigsum calculation to compute_derived and clean up
! 2006-07-31 kleist - analysis variable changed from ln(ps) to ps
! 2006-08-16 parrish - correct adjoint error, which only has impact when
! using strong balance constraint.
! 2008-05-28 safford - rm unused uses
! 2015-10-27 mahajan - code clean-up
!
! input argument list:
! rhnorm - normalized RH
! t - virtual temperature
! p - psfc
!
! output argument list:
! q - specific humidity
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block

use m_kinds, only: r_kind,i_kind
use derivsmod, only: dqdrh,dqdp,dqdt
use jfunc, only: qoption
use gridmod, only: lat2,lon2,nsig,regional
use guess_grids, only: ges_tsen,ntguessig
use constants, only: zero
implicit none

real(r_kind),intent(inout) :: rhnorm(lat2,lon2,nsig)
real(r_kind),intent(inout) :: t(lat2,lon2,nsig)
real(r_kind),intent(inout) :: p(lat2,lon2,nsig+1)
real(r_kind),intent(inout) :: q(lat2,lon2,nsig)

real(r_kind), parameter :: rmiss_th = -1.0e30

! local variables:
integer(i_kind) i,j,k

! Adjoint of convert normalized rh to q
do k=1,nsig
do j=1,lon2
do i=1,lat2
if(regional .and. ges_tsen(i,j,k,ntguessig) < rmiss_th) then
rhnorm(i,j,k) = zero
t(i,j,k ) = zero
p(i,j,k ) = zero
p(i,j,k+1) = zero
q(i,j,k) = zero
cycle
endif
rhnorm(i,j,k) = rhnorm(i,j,k) + dqdrh(i,j,k)*q(i,j,k)
if ( qoption == 2 ) then
t(i,j,k ) = t(i,j,k ) + dqdt(i,j,k)*q(i,j,k)
p(i,j,k ) = p(i,j,k ) - dqdp(i,j,k)*q(i,j,k)
p(i,j,k+1) = p(i,j,k+1) - dqdp(i,j,k)*q(i,j,k)
endif
q(i,j,k) = zero
enddo
enddo
enddo

return

end subroutine normal_rh_to_q_ad
12 changes: 8 additions & 4 deletions sorc/_workaround_/saber/gsi_covariance_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -988,6 +988,7 @@ subroutine cvfix_(gsicv,jedicv,vflip,need,ntimes,which)
!
real(kind=kind_real), allocatable :: t_pt(:,:,:)
real(kind=kind_real), pointer :: tv(:,:,:)=>NULL()
real(kind=kind_real), pointer :: t(:,:,:)=>NULL()
real(kind=kind_real), pointer :: tv_pt(:,:,:)=>NULL()
real(kind=kind_real), pointer :: q(:,:,:)=>NULL()
real(kind=kind_real), pointer :: q_pt(:,:,:)=>NULL()
Expand All @@ -1002,6 +1003,7 @@ subroutine cvfix_(gsicv,jedicv,vflip,need,ntimes,which)
! from first guess ...
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'q' ,q ,ier)
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'tv',tv,ier)
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'tsen',t,ier)
! from GSI cv ...
call gsi_bundlegetpointer(gsicv%step(ii),'q' ,q_pt ,ier)
call gsi_bundlegetpointer(gsicv%step(ii),'tv',tv_pt,ier)
Expand All @@ -1020,7 +1022,7 @@ subroutine cvfix_(gsicv,jedicv,vflip,need,ntimes,which)
endif
! retrieve missing field
if(which=='tlm') then
call gsi_tv_to_t_tl(tv,tv_pt,q,q_pt,t_pt)
call gsi_tv_to_t_tl(tv,tv_pt,q,q_pt,t_pt,t)
! pass it back to JEDI ...
allocate(aux1(size(rank2,2)))
if (vflip) then
Expand All @@ -1040,7 +1042,7 @@ subroutine cvfix_(gsicv,jedicv,vflip,need,ntimes,which)
endwhere
endif
if(which=='adm') then
call gsi_tv_to_t_ad(tv,tv_pt,q,q_pt,t_pt)
call gsi_tv_to_t_ad(tv,tv_pt,q,q_pt,t_pt,t)
where(need=='tv')
need='filled-'//need
endwhere
Expand Down Expand Up @@ -1101,6 +1103,7 @@ subroutine svfix_(gsisv,jedicv,vflip,need,ntimes,which)
!
real(kind=kind_real), allocatable :: t_pt(:,:,:)
real(kind=kind_real), pointer :: tv(:,:,:)=>NULL()
real(kind=kind_real), pointer :: t(:,:,:)=>NULL()
real(kind=kind_real), pointer :: tv_pt(:,:,:)=>NULL()
real(kind=kind_real), pointer :: q(:,:,:)=>NULL()
real(kind=kind_real), pointer :: q_pt(:,:,:)=>NULL()
Expand Down Expand Up @@ -1129,6 +1132,7 @@ subroutine svfix_(gsisv,jedicv,vflip,need,ntimes,which)
! from first guess ...
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'q' ,q ,ier)
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'tv',tv,ier)
call gsi_bundlegetpointer(gsi_metguess_bundle(ii),'tsen',t,ier)
! from GSI cv ...
call gsi_bundlegetpointer(gsisv(ii),'q' ,q_pt ,ier)
call gsi_bundlegetpointer(gsisv(ii),'tv',tv_pt,ier)
Expand All @@ -1147,7 +1151,7 @@ subroutine svfix_(gsisv,jedicv,vflip,need,ntimes,which)
endif
! retrieve missing field
if(which=='tlm') then
call gsi_tv_to_t_tl(tv,tv_pt,q,q_pt,t_pt)
call gsi_tv_to_t_tl(tv,tv_pt,q,q_pt,t_pt,t)
where(need=='tv')
need='filled-'//need
endwhere
Expand All @@ -1167,7 +1171,7 @@ subroutine svfix_(gsisv,jedicv,vflip,need,ntimes,which)
deallocate(aux1)
endif
if(which=='adm') then
call gsi_tv_to_t_ad(tv,tv_pt,q,q_pt,t_pt)
call gsi_tv_to_t_ad(tv,tv_pt,q,q_pt,t_pt,t)
where(need=='tv')
need='filled-'//need
endwhere
Expand Down