Skip to content

Remove discontinuities in the balance op #885

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

Merged
travissluka merged 9 commits into from
May 19, 2023
Merged

Remove discontinuities in the balance op #885

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
8484154
set to min/max instead of 0
guillaumevernieres May 8, 2023
dd0c5b9
...
guillaumevernieres May 8, 2023
2fce3cf
vanilla vert conv
guillaumevernieres May 10, 2023
a77d360
tidy sign
guillaumevernieres May 10, 2023
6ceb6d2
removed jac_mask
guillaumevernieres May 10, 2023
2ae9fda
tanh filt
guillaumevernieres May 11, 2023
0bfd935
updated yamls and test refs
guillaumevernieres May 11, 2023
79b18b7
merge develop and update answers after the develop history MESS i cre…
travissluka May 12, 2023
7ccd65d
code tidy/removed vanilla
guillaumevernieres May 18, 2023
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
92 changes: 41 additions & 51 deletions src/soca/LinearVariableChange/Balance/soca_balance_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,16 @@ module soca_balance_mod
! ------------------------------------------------------------------------------


function soca_tanh_filt(l, l0) result (coef)
real(kind=kind_real), intent(in) :: l
real(kind=kind_real), intent(in) :: l0

real(kind=kind_real) :: coef

coef = 0.5_kind_real*(tanh(l-l0)+1.0_kind_real)

end function soca_tanh_filt

! ------------------------------------------------------------------------------
!> Initialization of the balance operator and its trajectory.
!!
Expand All @@ -75,16 +85,12 @@ subroutine soca_balance_setup(self, f_conf, traj, geom)
integer :: isc, iec, jsc, jec
integer :: isd, ied, jsd, jed
integer :: i, j, k, nl
real(kind=kind_real), allocatable :: jac(:)
real(kind=kind_real), allocatable :: jac(:), coef_mld, coef_layers
type(soca_field), pointer :: tocn, socn, hocn, cicen, mld, layer_depth

! declarations related to the dynamic height Jacobians
character(len=:), allocatable :: filename, mask_name
real(kind=kind_real), allocatable :: jac_mask(:,:) !> mask for Jacobian
character(len=:), allocatable :: filename
real(kind=kind_real) :: threshold
integer :: nlayers !> dynamic height Jac=0 in nlayers upper layers
logical :: mask_detadt = .false. !> if true, set deta/dt to 0
logical :: mask_detads = .false. !> if true, set deta/ds to 0

! declarations related to the sea-ice Jacobian
character(len=:), allocatable :: kct_name
Expand All @@ -98,33 +104,6 @@ subroutine soca_balance_setup(self, f_conf, traj, geom)
isd=geom%isd; ied=geom%ied
jsd=geom%jsd; jed=geom%jed

! Setup mask for Jacobians related to the dynamic height balance
allocate(jac_mask(isd:ied,jsd:jed))
jac_mask = 1.0_kind_real
nlayers = 0
if ( f_conf%has("jac_mask") ) then
jac_mask = 0.0_kind_real
call f_conf%get_or_die("jac_mask.filename", filename)
call f_conf%get_or_die("jac_mask.name", mask_name)
call f_conf%get_or_die("jac_mask.threshold", threshold)
call f_conf%get_or_die("jac_mask.nlayers", nlayers)
call f_conf%get_or_die("jac_mask.detadt", mask_detadt)
call f_conf%get_or_die("jac_mask.detads", mask_detads)
call fms_io_init()
call read_data(filename, mask_name, jac_mask, domain=geom%Domain%mpp_domain)
call fms_io_exit()
where(jac_mask<threshold)
jac_mask = 0.0_kind_real
end where
end if

! Get configuration for Kst
call f_conf%get_or_die("dsdtmax", self%kst%dsdtmax)
call f_conf%get_or_die("dsdzmin", self%kst%dsdzmin)
call f_conf%get_or_die("dtdzmin", self%kst%dtdzmin)
call f_conf%get_or_die("nlayers", self%kst%nlayers) ! Set jac to 0 in the
! nlayers top layers

! Get required fields
call traj%get("tocn", tocn)
call traj%get("socn", socn)
Expand All @@ -139,27 +118,45 @@ subroutine soca_balance_setup(self, f_conf, traj, geom)
allocate(jac(nl))
self%kst%jacobian=0.0

! Get configuration for Kst
self%kst%dsdtmax = 0.1_kind_real
self%kst%dsdzmin = 3.0e-6_kind_real
self%kst%dtdzmin = 1.0e-6_kind_real
self%kst%nlayers = -999 ! input to the tanh filter
if ( f_conf%has("kst") ) then
call f_conf%get_or_die("kst.dsdtmax", self%kst%dsdtmax)
call f_conf%get_or_die("kst.dsdzmin", self%kst%dsdzmin)
call f_conf%get_or_die("kst.dtdzmin", self%kst%dtdzmin)
call f_conf%get_or_die("kst.nlayers", self%kst%nlayers)
end if

! Compute and store Jacobian of Kst
do i = isc, iec
do j = jsc, jec
jac=0.0
! do nothing if on land
if ( geom%mask2d(i, j) == 0 ) cycle

! compute dS(T)/dT
call soca_soft_jacobian(jac,&
&tocn%val(i,j,:),&
&socn%val(i,j,:),&
&hocn%val(i,j,:),&
&self%kst%dsdtmax, self%kst%dsdzmin, self%kst%dtdzmin)
! Set Jacobian to 0 above mixed layer

! filter out the Jacobian as specified in the configuration
do k=1,nl
if (layer_depth%val(i,j,k) < mld%val(i,j,1)) then
jac(k) = 0.0_kind_Real
end if
coef_mld = soca_tanh_filt(layer_depth%val(i,j,k),mld%val(i,j,1))

@travissluka travissluka May 12, 2023

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know you don't totally agree with how the nemovar balances were implemented, but I believe they set the jacobian to 0 at the surface and linearly interpolate down to the MLD, instead of such a sharp cutoff. Thoughts on trying that (someday?)

@guillaumevernieres guillaumevernieres May 12, 2023

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure.

coef_layers = soca_tanh_filt(real(k, kind=kind_real), real(self%kst%nlayers, kind=kind_real))
self%kst%jacobian(i,j,k) = jac(k)*coef_mld*coef_layers
end do
self%kst%jacobian(i,j,:) = jac(:)
self%kst%jacobian(i,j,1:self%kst%nlayers) = 0.0_kind_real
end do
end do
deallocate(jac)

! Get configuration for Ksshts
self%ksshts%nlayers = -999 ! input to the tanh filter
if ( f_conf%has("ksshts") ) call f_conf%get_or_die("ksshts.nlayers", self%ksshts%nlayers)

! Compute Jacobian of Ksshts
allocate(self%ksshts%kssht, mold=self%kst%jacobian)
allocate(self%ksshts%ksshs, mold=self%kst%jacobian)
Expand All @@ -176,21 +173,14 @@ subroutine soca_balance_setup(self, f_conf, traj, geom)
&hocn%val(i,j,k),&
&geom%lon(i,j),&
&geom%lat(i,j))
self%ksshts%kssht(i,j,k) = jac(1)*jac_mask(i,j)
self%ksshts%ksshs(i,j,k) = jac(2)*jac_mask(i,j)
end do
if (nlayers>0) then
self%ksshts%kssht(i,j,1:nlayers) = 0.0_kind_real
self%ksshts%ksshs(i,j,1:nlayers) = 0.0_kind_real
end if
coef_layers = soca_tanh_filt(real(k, kind=kind_real), real(self%ksshts%nlayers, kind=kind_real))
self%ksshts%kssht(i,j,k) = jac(1)*coef_layers
self%ksshts%ksshs(i,j,k) = jac(2)*coef_layers
end do
end do
end do
deallocate(jac)

! Zero-out Jacobians if required by configuration
if (mask_detadt) self%ksshts%kssht = 0.0_kind_real
if (mask_detads) self%ksshts%ksshs = 0.0_kind_real

! Compute Kct
if (traj%has("cicen")) then
! Setup dc/dT
Expand Down
1 change: 1 addition & 0 deletions src/soca/LinearVariableChange/Balance/soca_ksshts_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ module soca_ksshts_mod
type, public :: soca_ksshts
real(kind=kind_real), allocatable :: kssht(:,:,:) !< deta(i,j)/dT(i,j,k)
real(kind=kind_real), allocatable :: ksshs(:,:,:) !< deta(i,j)/dS(i,j,k)
integer :: nlayers
end type soca_ksshts


Expand Down
14 changes: 5 additions & 9 deletions src/soca/LinearVariableChange/Balance/soca_kst_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ subroutine soca_soft_jacobian (jac, t, s, h, dsdtmax, dsdzmin, dtdzmin)

real(kind=kind_real), allocatable :: dtdz(:), dsdz(:)
integer :: nl, z
real(kind=kind_real) :: j

! Allocate
nl = size(t,1)
Expand All @@ -59,20 +58,17 @@ subroutine soca_soft_jacobian (jac, t, s, h, dsdtmax, dsdzmin, dtdzmin)

jac = 0.0
do z=1,nl
jac(z) = 0.0

! Limit application of soft according to configuration
if ( abs(dtdz(z)) < dtdzmin ) cycle
if ( abs(dsdz(z)) < dsdzmin ) cycle
if ( abs(dtdz(z)) < dtdzmin ) dtdz(z) = sign(dtdzmin, dtdz(z))

if ( abs(dsdz(z)) < dsdzmin ) dsdz(z) = sign(dsdzmin, dsdz(z))

! Jacobian of soft
j=dsdz(z)/dtdz(z)
jac(z)=dsdz(z)/dtdz(z)

! Limit application of soft according to configuration
if ( abs(j) > dsdtmax ) cycle

! if we reach this point in the code, the jacobian is usable
jac(z) = j;
if ( abs(jac(z)) > dsdtmax ) jac(z) = sign(dsdtmax, jac(z))
end do

end subroutine soca_soft_jacobian
Expand Down
8 changes: 8 additions & 0 deletions src/soca/LinearVariableChange/VertConv/soca_vertconv_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ module soca_vertconv_mod
real(kind=kind_real) :: lz_min !> Vertical decorrelation minimum [m]
real(kind=kind_real) :: lz_mld !> if /= 0, Use MLD to calculate Lz
real(kind=kind_real) :: lz_mld_max !> if calculating Lz from MLD, max value to use
logical :: vanilla !> if true, only use the layer scaling
real(kind=kind_real) :: scale_layer_thick !> Set the minimum decorrelation scale
!> as a multiple of the layer thickness
type(soca_state), pointer :: bkg !> Background
Expand Down Expand Up @@ -67,6 +68,10 @@ subroutine soca_conv_setup (self, bkg, geom, f_conf)
if ( self%lz_mld /= 0) &
call f_conf%get_or_die("Lz_mld_max", self%lz_mld_max )
call f_conf%get_or_die("scale_layer_thick", self%scale_layer_thick )
self%vanilla = .false.
if ( f_conf%has("vanilla") ) then
call f_conf%get_or_die("vanilla", self%vanilla )
end if

! Associate background and geometry
self%bkg => bkg
Expand Down Expand Up @@ -100,6 +105,9 @@ subroutine soca_calc_lz(self, i, j, lz)
lz = self%lz_min
lz = max(lz, self%scale_layer_thick*abs(hocn%val(i,j,:)))

! if doing vanilla vertconv, we're done!
if (self%vanilla) return

@travissluka travissluka May 12, 2023

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

whats the point of adding this flag? It's the same effect as setting Lz_mld: 0 in the yaml

@guillaumevernieres guillaumevernieres May 12, 2023

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤦‍♂️
🍨 ?


! if the upper Lz should be calculated from the MLD
! interpolate values from top to bottom of ML
if ( self%lz_mld /= 0 ) then
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/3dhyb.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,10 +70,13 @@ cost function:
scale_layer_thick: 1.5

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/3dvar_godas.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,10 +72,13 @@ cost function:
scale_layer_thick: 1.5

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 999

@travissluka travissluka May 12, 2023

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@guillaumevernieres I see several different values of nlayers depending on the tests here. What values would you recommend for the real-world config?

@guillaumevernieres guillaumevernieres May 12, 2023

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@travissluka We need to turn it off for the global until we properly test it and add a high latitude limit. 999 would turn it off. Omitting the kst config in the yaml to turnt kst off might be a better idea.
10 maybe in the hat10?

ksshts:
nlayers: 10
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
14 changes: 7 additions & 7 deletions test/testinput/3dvarbump.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,13 @@ cost function:
hicen_max: 100.0

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 10
# dcdt:
# filename: data_static/72x35x25/dcdt.nc
# name: dcdt
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2

observations:
observers:
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/3dvarfgat.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,10 +85,13 @@ cost function:
scale_layer_thick: 1.5

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 10
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/3dvarfgat_pseudo.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,10 +88,13 @@ cost function:
scale_layer_thick: 1.5

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 10
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/3dvarlowres_soca.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,13 @@ cost function:
filter variables: *a_vars

- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2

observations:
observers:
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/convertincrement.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,13 @@ linear variable change:
output variables: [tocn, socn, ssh, hocn]
linear variable changes:
- linear variable change name: BalanceSOCA
dsdtmax: 0.0
dsdzmin: 0.0
dtdzmin: 0.0
nlayers: 2
kst:
dsdtmax: 0.0
dsdzmin: 0.0
dtdzmin: 0.0
nlayers: 2
ksshts:
nlayers: 2

increments:
- date: 2018-04-15T00:00:00Z
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/dirac_soca_cov.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,10 +67,13 @@ background error:
output variables: *soca_vars

- linear variable change name: BalanceSOCA
dsdtmax: 1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
11 changes: 7 additions & 4 deletions test/testinput/dirac_socahyb_cov.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,10 +72,13 @@ background error:
input variables: *soca_vars
output variables: *soca_vars
- linear variable change name: BalanceSOCA
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
kst:
dsdtmax: 0.1
dsdzmin: 3.0e-6
dtdzmin: 1.0e-6
nlayers: 2
ksshts:
nlayers: 2
dcdt:
filename: data_static/72x35x25/dcdt.nc
name: dcdt
Expand Down
Loading