Skip to content

Commit 751f07f

Browse files
kshedstromkshedstrom
committed
Merge remote-tracking branch 'gfdl/dev/gfdl' into fix_obc_maskingdepth
2 parents b93bfd2 + 88306bd commit 751f07f

File tree

6 files changed

+1757
-792
lines changed

6 files changed

+1757
-792
lines changed

‎src/core/MOM.F90‎

Lines changed: 277 additions & 163 deletions
Large diffs are not rendered by default.

‎src/core/MOM_PressureForce_FV.F90‎

Lines changed: 162 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,8 @@ module MOM_PressureForce_FV
6161
!! exceeds the vertical grid spacing, reset intxpa at the interface below
6262
!! a trusted interior cell. (This often applies in ice shelf cavities.)
6363
logical :: MassWghtInterpVanOnly !< If true, don't do mass weighting of T/S interpolation unless vanished
64+
logical :: reset_intxpa_flattest !< If true, use flattest interface rather than top for reset integral
65+
!! in cases where no best nonvanished interface
6466
real :: h_nonvanished !< A minimal layer thickness that indicates that a layer is thick enough
6567
!! to usefully reestimate the pressure integral across the interface
6668
!! below it [H ~> m or kg m-2]
@@ -214,8 +216,12 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
214216
logical, dimension(SZI_(G),SZJB_(G)) :: &
215217
seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate
216218
! of the curvature terms in the inty_pa.
217-
218-
219+
real, dimension(SZIB_(G),SZJ_(G)) :: &
220+
delta_p_x ! If using flattest interface for reset integral, store x interface
221+
! differences [R L2 T-2 ~> Pa]
222+
real, dimension(SZI_(G),SZJB_(G)) :: &
223+
delta_p_y ! If using flattest interface for reset integral, store y interface
224+
! differences [R L2 T-2 ~> Pa]
219225
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
220226
MassWt_u ! The fractional mass weighting at a u-point [nondim].
221227
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
@@ -581,6 +587,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
581587
intx_za_nonlin(:,:) = 0.0 ; intx_za_cor_ri(:,:) = 0.0 ; dp_int_x(:,:) = 0.0
582588
do j=js,je ; do I=Isq,Ieq
583589
seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.)
590+
delta_p_x(I,j) = 0.0
584591
enddo ; enddo
585592

586593
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
@@ -619,15 +626,40 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
619626
enddo
620627

621628
if (do_more_k) then
622-
! There are still points where a correction is needed, so use the top interface.
623-
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
624-
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
625-
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
626-
p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
627-
intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
628-
dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
629+
if (CS%reset_intxpa_flattest) then
630+
! There are still points where a correction is needed, so use flattest interface
631+
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
632+
! choose top layer first
633+
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
634+
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
635+
p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
636+
intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
637+
dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
638+
delta_p_x(I,j) = abs(p(i+1,j,1)-p(i,j,1))
639+
do k=1,nz
640+
if (abs(p(i+1,j,k+1)-p(i,j,k+1)) < delta_p_x(I,j)) then
641+
! bottom of layer is less sloped than top. Use this layer
642+
delta_p_x(I,j) = abs(p(i+1,j,k+1)-p(i,j,k+1))
643+
T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
644+
S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
645+
p_int_W(I,j) = p(i,j,K+1) ; p_int_E(I,j) = p(i+1,j,K+1)
646+
intx_za_nonlin(I,j) = intx_za(I,j,K+1) - 0.5*(za(i,j,K+1) + za(i+1,j,K+1))
647+
dp_int_x(I,j) = p(i+1,j,K+1)-p(i,j,K+1)
648+
endif
649+
enddo
629650
seek_x_cor(I,j) = .false.
630-
endif ; enddo ; enddo
651+
endif; enddo; enddo;
652+
else
653+
! There are still points where a correction is needed, so use the top interface.
654+
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
655+
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
656+
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
657+
p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
658+
intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
659+
dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
660+
seek_x_cor(I,j) = .false.
661+
endif ; enddo ; enddo
662+
endif
631663
endif
632664

633665
do j=js,je
@@ -658,6 +690,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
658690
inty_za_nonlin(:,:) = 0.0 ; inty_za_cor_ri(:,:) = 0.0 ; dp_int_y(:,:) = 0.0
659691
do J=Jsq,Jeq ; do i=is,ie
660692
seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.)
693+
delta_p_y(i,J) = 0.0
661694
enddo ; enddo
662695

663696
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
@@ -695,15 +728,40 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
695728
enddo
696729

697730
if (do_more_k) then
698-
! There are still points where a correction is needed, so use the top interface.
699-
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
700-
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
701-
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
702-
p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
703-
inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
704-
dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
705-
seek_y_cor(i,J) = .false.
706-
endif ; enddo ; enddo
731+
if (CS%reset_intxpa_flattest) then
732+
! There are still points where a correction is needed, so use flattest interface.
733+
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
734+
! choose top interface first
735+
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
736+
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
737+
p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
738+
inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
739+
dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
740+
delta_p_y(i,J) = abs(p(i,j+1,1)-p(i,j,1))
741+
do k=1,nz
742+
if (abs(p(i,j+1,k+1)-p(i,j,k+1)) < delta_p_y(i,J)) then
743+
! bottom of layer is less sloped than top. Use this layer
744+
delta_p_y(i,J) = abs(p(i,j+1,k+1)-p(i,j,k+1))
745+
T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
746+
S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
747+
p_int_S(i,J) = p(i,j,K+1) ; p_int_N(i,J) = p(i,j+1,K+1)
748+
inty_za_nonlin(i,J) = inty_za(i,J,K+1) - 0.5*(za(i,j,K+1) + za(i,j+1,K+1))
749+
dp_int_y(i,J) = p(i,j+1,K+1) - p(i,j,K+1)
750+
endif
751+
enddo
752+
seek_y_cor(i,J) = .false.
753+
endif ; enddo ; enddo
754+
else
755+
! There are still points where a correction is needed, so use the top interface.
756+
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
757+
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
758+
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
759+
p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
760+
inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
761+
dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
762+
seek_y_cor(i,J) = .false.
763+
endif ; enddo ; enddo
764+
endif
707765
endif
708766

709767
do J=Jsq,Jeq
@@ -946,6 +1004,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
9461004
logical, dimension(SZI_(G),SZJB_(G)) :: &
9471005
seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate
9481006
! of the curvature terms in the inty_pa.
1007+
real, dimension(SZIB_(G),SZJ_(G)) :: &
1008+
delta_z_x ! If using flattest interface for reset integral, store x interface differences [Z ~> m]
1009+
real, dimension(SZI_(G),SZJB_(G)) :: &
1010+
delta_z_y ! If using flattest interface for reset integral, store y interface differences [Z ~> m]
9491011

9501012
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
9511013
T_tmp, & ! Temporary array of temperatures where layers that are lighter
@@ -1472,6 +1534,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
14721534
intx_pa_nonlin(:,:) = 0.0 ; dgeo_x(:,:) = 0.0 ; intx_pa_cor_ri(:,:) = 0.0
14731535
do j=js,je ; do I=Isq,Ieq
14741536
seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.)
1537+
delta_z_x(I,j) = 0.0
14751538
enddo ; enddo
14761539

14771540
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
@@ -1514,16 +1577,43 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
15141577
enddo
15151578

15161579
if (do_more_k) then
1517-
! There are still points where a correction is needed, so use the top interface for lack of a better idea?
1518-
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
1519-
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
1520-
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
1521-
p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1522-
p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
1523-
intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
1524-
dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
1525-
seek_x_cor(I,j) = .false.
1526-
endif ; enddo ; enddo
1580+
if (CS%reset_intxpa_flattest) then
1581+
! There are still points where a correction is needed, so use flattest interface
1582+
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
1583+
! choose top layer first
1584+
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
1585+
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
1586+
p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1587+
p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
1588+
intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
1589+
dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
1590+
delta_z_x(I,j) = abs(e(i+1,j,1)-e(i,j,1))
1591+
do k=1,nz
1592+
if (abs(e(i+1,j,k+1)-e(i,j,k+1)) < delta_z_x(I,j)) then
1593+
! bottom of layer is less sloped than top. Use this layer
1594+
delta_z_x(I,j) = abs(e(i+1,j,k+1)-e(i,j,k+1))
1595+
T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
1596+
S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
1597+
p_int_W(I,j) = -GxRho0*(e(i,j,K+1) - Z_0p(i,j))
1598+
p_int_E(I,j) = -GxRho0*(e(i+1,j,K+1) - Z_0p(i,j))
1599+
intx_pa_nonlin(I,j) = intx_pa(I,j,K+1) - 0.5*(pa(i,j,K+1) + pa(i+1,j,K+1))
1600+
dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,K+1)-e(i,j,K+1))
1601+
endif
1602+
enddo
1603+
seek_x_cor(I,j) = .false.
1604+
endif ; enddo ; enddo
1605+
else
1606+
! There are still points where a correction is needed, so use the top interface for lack of a better idea?
1607+
do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
1608+
T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
1609+
S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
1610+
p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1611+
p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
1612+
intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
1613+
dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
1614+
seek_x_cor(I,j) = .false.
1615+
endif ; enddo ; enddo
1616+
endif
15271617
endif
15281618

15291619
do j=js,je
@@ -1554,6 +1644,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
15541644
inty_pa_nonlin(:,:) = 0.0 ; dgeo_y(:,:) = 0.0 ; inty_pa_cor_ri(:,:) = 0.0
15551645
do J=Jsq,Jeq ; do i=is,ie
15561646
seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.)
1647+
delta_z_y(i,J) = 0.0
15571648
enddo ; enddo
15581649

15591650
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
@@ -1595,16 +1686,43 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
15951686
enddo
15961687

15971688
if (do_more_k) then
1598-
! There are still points where a correction is needed, so use the top interface for lack of a better idea?
1599-
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
1600-
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
1601-
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
1602-
p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1603-
p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
1604-
inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
1605-
dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
1606-
seek_y_cor(i,J) = .false.
1607-
endif ; enddo ; enddo
1689+
if (CS%reset_intxpa_flattest) then
1690+
! There are still points where a correction is needed, so use flattest interface.
1691+
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
1692+
! choose top interface first
1693+
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
1694+
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
1695+
p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1696+
p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
1697+
inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
1698+
dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
1699+
delta_z_y(i,J) = abs(e(i,j+1,1)-e(i,j,1))
1700+
do k=1,nz
1701+
if (abs(e(i,j+1,k+1)-e(i,j,k+1)) < delta_z_y(i,J)) then
1702+
! bottom of layer is less sloped than top. Use this layer
1703+
delta_z_y(i,J) = abs(e(i,j+1,k+1)-e(i,j,k+1))
1704+
T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
1705+
S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
1706+
p_int_S(i,J) = -GxRho0*(e(i,j,k+1) - Z_0p(i,j))
1707+
p_int_N(i,J) = -GxRho0*(e(i,j+1,k+1) - Z_0p(i,j))
1708+
inty_pa_nonlin(i,J) = inty_pa(i,J,k+1) - 0.5*(pa(i,j,k+1) + pa(i,j+1,k+1))
1709+
dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,k+1)-e(i,j,k+1))
1710+
endif
1711+
enddo
1712+
seek_y_cor(i,J) = .false.
1713+
endif ; enddo ; enddo
1714+
else
1715+
! There are still points where a correction is needed, so use the top interface for lack of a better idea?
1716+
do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
1717+
T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
1718+
S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
1719+
p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
1720+
p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
1721+
inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
1722+
dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
1723+
seek_y_cor(i,J) = .false.
1724+
endif ; enddo ; enddo
1725+
endif
16081726
endif
16091727

16101728
do J=Jsq,Jeq
@@ -1940,9 +2058,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
19402058
call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL", CS%reset_intxpa_integral, &
19412059
"If true, reset INTXPA to match pressures at first nonvanished cell. "//&
19422060
"Includes pressure correction.", default=.false., do_not_log=.not.use_EOS)
2061+
call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL_FLATTEST", CS%reset_intxpa_flattest, &
2062+
"If true, use flattest interface as reference interface where there is no "//&
2063+
"better choice for RESET_INTXPA_INTEGRAL. Otherwise, use surface interface.", &
2064+
default=.false., do_not_log=.not.use_EOS)
19432065
if (.not.use_EOS) then ! These options do nothing without an equation of state.
19442066
CS%correction_intxpa = .false.
19452067
CS%reset_intxpa_integral = .false.
2068+
CS%reset_intxpa_flattest = .false.
19462069
endif
19472070
call get_param(param_file, mdl, "RESET_INTXPA_H_NONVANISHED", CS%h_nonvanished, &
19482071
"A minimal layer thickness that indicates that a layer is thick enough to usefully "//&

0 commit comments

Comments
 (0)