Skip to content

Fix zero boundary and orbital velocity issue in the cap #1402

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

Open
wants to merge 3 commits into
base: dev/ufs-weather-model
Choose a base branch
from
Open

Fix zero boundary and orbital velocity issue in the cap #1402

Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2241e43
Fix zero boundary and orbital velocity issue in the cap
danishyo Apr 4, 2025
72b634a
Revert mapsta if condition within CESM block
danishyo Apr 14, 2025
444aa2b
Revert changes for unused fields
danishyo Apr 14, 2025
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
52 changes: 24 additions & 28 deletions model/src/wav_import_export.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -846,7 +846,7 @@ subroutine export_fields (gcomp, rc)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1)
iy = mapsf(isea,2)
if (mapsta(iy,ix) == 1) then
if (mapsta(iy,ix) /= 0) then
sw_tauox(jsea) = tauox(jsea)
sw_tauoy(jsea) = tauoy(jsea)
endif
Expand Down Expand Up @@ -883,7 +883,7 @@ subroutine export_fields (gcomp, rc)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1)
iy = mapsf(isea,2)
if (mapsta(iy,ix) == 1) then
if (mapsta(iy,ix) /= 0) then
sw_wnmean(jsea) = wnmean(jsea)
endif
enddo
Expand All @@ -902,7 +902,7 @@ subroutine export_fields (gcomp, rc)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1)
iy = mapsf(isea,2)
if (mapsta(iy,ix) == 1) then
if (mapsta(iy,ix) /= 0) then
sw_taubblx(jsea) = taubbl(jsea,1)
sw_taubbly(jsea) = taubbl(jsea,2)
endif
Expand Down Expand Up @@ -1183,8 +1183,8 @@ subroutine CalcRoughl ( wrln)
ix = mapsf(isea,1)
iy = mapsf(isea,2)
if ( firstCall ) then
if(( runtype == 'initial' .and. mapsta(iy,ix) == 1 ) .or. &
( runtype == 'continue' .and. abs(mapsta(iy,ix)) == 1 )) then
if(( runtype == 'initial' .and. mapsta(iy,ix) /= 0 ) .or. &
Copy link
Contributor

@DeniseWorthen DeniseWorthen Apr 14, 2025

Choose a reason for hiding this comment

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

@danishyo Did you intend to change this block? I don't remember that coastal app exports Z0 and this is w/in code that UFS global mode uses.

Copy link
Author

@danishyo danishyo Apr 14, 2025

Choose a reason for hiding this comment

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

No, I didn't, I will revert this one as well.

( runtype == 'continue' .and. abs(mapsta(iy,ix)) /= 0 )) then
charn(jsea) = zero
llws(:) = .true.
ustar = zero
Expand Down Expand Up @@ -1243,7 +1243,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
sxx1 = 0.0
syy1 = 0.0
sxy1 = 0.0
Expand Down Expand Up @@ -1271,7 +1271,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval)
syy1 = syy1 + fte * abyy/cg(nk,isea)
sxy1 = sxy1 + fte * abxy/cg(nk,isea)
end if
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
sxxn(jsea) = sxx1*dwat*grav
syyn(jsea) = syy1*dwat*grav
sxyn(jsea) = sxy1*dwat*grav
Expand Down Expand Up @@ -1323,7 +1323,7 @@ subroutine CalcEF (a, wave_elevation_spectrum)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
factor = dden(ik) / cg(ik,isea)
ebd = ab(jsea) * factor
ebd = ebd / dsii(ik)
Expand Down Expand Up @@ -1367,7 +1367,7 @@ subroutine CalcHS (a, hs, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
et = 0.0
do ik = 1,nk
factor = dden(ik) / cg(ik,isea)
Expand Down Expand Up @@ -1425,7 +1425,7 @@ subroutine CalcBHD (a, bhd, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
ebd = 0.0
bhd1 = 0.0
do ik = 1,nk
Expand Down Expand Up @@ -1481,7 +1481,7 @@ subroutine CalcStokes(a, us, vs, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
us1 = 0.0
vs1 = 0.0
do ik = 1,nk
Expand All @@ -1492,8 +1492,8 @@ subroutine CalcStokes(a, us, vs, fval)
abx = abx + a(ith,ik,jsea)*ecos(ith)
aby = aby + a(ith,ik,jsea)*esin(ith)
end do
kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
if (kd .lt. 6.0) then
kd = max ( 1.0e-7 , wn(ik,isea) * dw(isea) ) !0.001
Copy link
Collaborator

@JessicaMeixner-NOAA JessicaMeixner-NOAA Apr 25, 2025

Choose a reason for hiding this comment

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

Why are there changes to this calculation? This appears different than the mapsta issue updated elsewhere but I don't see anything about it in the PR description.

Copy link
Author

@danishyo danishyo May 1, 2025

Choose a reason for hiding this comment

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

I borrow routine coefficients from https://github.com/NOAA-EMC/WW3/blob/develop/model/src/wmesmfmd.F90#L7315
For mapsta, I am OK to change ">0" if this includes ocean and boundary points.

Copy link
Contributor

@DeniseWorthen DeniseWorthen May 1, 2025

Choose a reason for hiding this comment

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

The code for the ussco matches (or should) what is currently in w3iogomd.F90

Copy link
Collaborator

@JessicaMeixner-NOAA JessicaMeixner-NOAA May 2, 2025

Choose a reason for hiding this comment

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

@danishyo - Can you let us know how your testing goes for ">0" this should be sufficient.

I would think this calculation should match what is in w3iogo. I'd have to look at the differences between that and wmesmf as to why that would have been diverged. Let me know if you need help investigating this further @danishyo or if you will take the lead on that.

Copy link
Collaborator

@JessicaMeixner-NOAA JessicaMeixner-NOAA Jul 29, 2025

Choose a reason for hiding this comment

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

@danishyo were you ever able to run those tests?

Copy link
Author

@danishyo danishyo Aug 7, 2025

Choose a reason for hiding this comment

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

@danishyo were you ever able to run those tests?

I didn't do "mapsta>0" test, but base on info from manual, I believe this should fix the issue as well.

Copy link
Collaborator

@JessicaMeixner-NOAA JessicaMeixner-NOAA Aug 7, 2025

Choose a reason for hiding this comment

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

Let us know when you've been able to run those tests and confirm if that works or not.

We also have outstanding issues about this part of the code change and why it's diverged, etc.

This PR is waiting on those updates to move forward.

if (kd .lt. 18.0) then !6.0
fkd = factor / sinh(kd)**2
ussco = fkd*sig(ik)*wn(ik,isea)*cosh(2.0*kd)
else
Expand Down Expand Up @@ -1545,7 +1545,7 @@ subroutine CalcUVBed(a, ubrx, ubry, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
uba1 = 0.0
ubd1 = 0.0
ubr1 = 0.0
Expand All @@ -1559,23 +1559,19 @@ subroutine CalcUVBed(a, ubrx, ubry, fval)
abx = abx + a(ith,ik,jsea)*ecos(ith)
aby = aby + a(ith,ik,jsea)*esin(ith)
end do
kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
if (kd .lt. 6.0) then
!kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
kd = max ( 1.0e-7, min(18.,wn(ik,isea)) * dw(isea) )
!if (kd .lt. 6.0) then
fkd = factor / sinh(kd)**2
ubr1 = ubr1 + ab*sig(ik)**2 * fkd
uba1 = uba1 + abx*sig(ik)**2 * fkd
ubd1 = ubd1 + aby*sig(ik)**2 * fkd
end if
!end if
end do !ik
ubr1 = sqrt(2.0*max(0.0,ubr1))
if (ubr1 .ge. 1.0e-7) then
ubd1 = atan2(ubd1,uba1)
else
ubd1 = 0.0
end if
uba1 = ubr1
ubrx(jsea) = uba1*cos(ubd1)
ubry(jsea) = uba1*sin(ubd1)
ubr1 = sqrt(2.0*max(0.0,ubr1)) !make ubr1 >=0 already
ubd1 = atan2(ubd1,uba1)
ubrx(jsea) = ubr1*cos(ubd1)
ubry(jsea) = ubr1*sin(ubd1)
else
ubrx(jsea) = fval
ubry(jsea) = fval
Expand Down Expand Up @@ -1615,7 +1611,7 @@ subroutine CalcTHM (a, thm, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
etx = 0.0
ety = 0.0
do ik = 1,nk
Expand Down Expand Up @@ -1676,7 +1672,7 @@ subroutine CalcT0M1 (a, t0m1, fval)
call init_get_isea(isea, jsea)
ix = mapsf(isea,1) ! global ix
iy = mapsf(isea,2) ! global iy
if (mapsta(iy,ix) == 1) then ! active sea point
if (mapsta(iy,ix) /= 0) then ! active sea point
etr = 0.0
et = 0.0
do ik = 1,nk
Expand Down