Skip to content

Necessary changes for land implicit coupling #83

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
JosephMouallem wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Necessary changes for land implicit coupling #83

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2a57937
flip matrix
JosephMouallem Mar 27, 2026
0d845d4
populate statemid dkt
JosephMouallem Mar 27, 2026
079ffb4
add surf_diff related routines
JosephMouallem Mar 27, 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
1 change: 1 addition & 0 deletions GFS_layer/GFS_physics_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1573,6 +1573,7 @@ subroutine GFS_physics_driver_down &

endif
endif
statemid%dkt(:,:)=dkt(:,:)

! --- ... return updated smsoil and stsoil to global arrays
Sfcprop%smc(:,:) = smsoil(:,:)
Expand Down
3 changes: 3 additions & 0 deletions GFS_layer/GFS_typedefs.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1032,6 +1032,7 @@ module GFS_typedefs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

real(kind=kind_phys), allocatable :: stored_au_out(:,:), stored_f1_out(:,:), stored_f2_out(:,:), stored_diss_out(:,:)
real(kind=kind_phys), allocatable :: dkt(:,:)
integer, allocatable :: stored_kpbl(:)
real(kind=kind_phys), allocatable :: stored_flux_cg(:,:), stored_flux_en(:,:), stored_elm_pbl(:,:)
real(kind=kind_phys), allocatable :: stored_dudt(:,:), stored_dvdt(:,:), stored_dqdt(:,:,:)
Expand Down Expand Up @@ -4041,6 +4042,7 @@ subroutine statemid_create (Statemid, IM, Model)
allocate (Statemid%stored_dvsfc1(IM))
allocate (Statemid%stored_dtsfc1(IM))
allocate (Statemid%stored_dqsfc1(IM))
allocate (Statemid%dkt(IM,model%levs))

! Surface fluxes and state
allocate (Statemid%stored_hflx(IM), Statemid%stored_evap(IM), Statemid%stored_stress(IM), Statemid%stored_wind(IM))
Expand All @@ -4058,6 +4060,7 @@ subroutine statemid_create (Statemid, IM, Model)
Statemid%stored_snowc = clear_val
Statemid%stored_drain = clear_val
Statemid%stored_runof = clear_val
Statemid%dkt = clear_val
!Statemid%stored_ep1d = clear_val
!Statemid%stored_gflx = clear_val
Statemid%stored_kpbl = 1
Expand Down
225 changes: 225 additions & 0 deletions atmos_shared/surf_diff.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,225 @@
module surf_diff_mod

! Mostly taken from vert_diff in atmos_phys
! could be replaced by the original routines
! when merging SHiELD and AMx
! Scientific description can be found in the
! document vert_diff.tech.ps dated March 2001
! joseph.mouallem@noaa.gov
! March 2026


use constants_mod, only: GRAV, RDGAS, RVGAS, CP_AIR
use IPD_typedefs, only: IPD_data_type, kind_phys
use fv_iau_mod, only: IAU_external_data_type
use block_control_mod, only: block_control_type

implicit none

private

public :: surf_diff_type, compute_nu, compute_e, vert_diff_down_2

real(kind=kind_phys), parameter :: d608 = (RVGAS-RDGAS)/RDGAS
logical :: use_virtual_temp_vert_diff = .true.
logical :: do_mcm_plev = .false.



type surf_diff_type

real, pointer, dimension(:,:) :: dtmass => NULL(), &
dflux_t => NULL(), &
delta_t => NULL(), &
delta_u => NULL(), &
delta_v => NULL()
real, pointer, dimension(:,:,:) :: tdt_dyn => NULL(), &
qdt_dyn => NULL(), &
dgz_dyn => NULL(), &
ddp_dyn => NULL(), &
tdt_rad => NULL() !miz

real, pointer, dimension(:,:,:) :: dflux_tr => NULL(),& ! tracer flux tendency
delta_tr => NULL() ! tracer tendency
end type surf_diff_type
!
!
contains
!
!

! inspired by atmos_phys/atmos_param/vert_diff/vert_diff.F90
subroutine compute_nu (diff, p_half, p_full, z_full, t, q, nu)

!-----------------------------------------------------------------------
real, intent(in) , dimension(:,:,:) :: diff, p_half, p_full, &
z_full, t, q
real, intent(out) , dimension(:,:,:) :: nu

real, dimension(size(t,1),size(t,2),size(t,3)) :: rho_half, tt
integer :: nlev
!-----------------------------------------------------------------------

nlev = size(nu,3)

if ( use_virtual_temp_vert_diff ) then
tt = t * (1.0 + d608*q) ! virtual temperature
else
tt = t ! Take out virtual temperature effect here to mimic supersource
endif

rho_half(:,:,2:nlev) = & ! density at half levels
2.0*p_half(:,:,2:nlev)/(RDGAS*(tt(:,:,2:nlev)+tt(:,:,1:nlev-1)))

if(do_mcm_plev) then
nu(:,:,2:nlev) = GRAV*rho_half(:,:,2:nlev)*rho_half(:,:,2:nlev)*diff(:,:,2:nlev)/ &
(p_full(:,:,2:nlev)-p_full(:,:,1:nlev-1))
else
nu(:,:,2:nlev) = rho_half(:,:,2:nlev)*diff(:,:,2:nlev) / &
(z_full(:,:,1:nlev-1)-z_full(:,:,2:nlev))
endif
!-----------------------------------------------------------------------

end subroutine compute_nu

subroutine compute_e (delt, mu, nu, e, a, b, c, g)

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

real, intent(in) :: delt
real, intent(in) , dimension(:,:,:) :: mu, nu
real, intent(out) , dimension(:,:,:) :: e, a, b, c, g

integer :: k, nlev

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

nlev = size(mu,3)

a(:,:,1:nlev-1) = - mu(:,:,1:nlev-1)*nu(:,:,2:nlev)*delt
a(:,:,nlev ) = 0.0
c(:,:,2:nlev ) = - mu(:,:,2:nlev )*nu(:,:,2:nlev)*delt
c(:,:,1 ) = 0.0

b = 1.0 - a - c

e(:,:,1) = - a(:,:,1)/b(:,:,1)
do k= 2, nlev - 1
g(:,:,k) = 1.0/(b(:,:,k) + c(:,:,k)*e(:,:,k-1))
e(:,:,k) = - a(:,:,k)*g(:,:,k)
enddo

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

end subroutine compute_e

subroutine compute_f (dt_xi, b, c, g, f)

!-----------------------------------------------------------------------
real, intent(in) , dimension(:,:,:) :: dt_xi, b, c, g
real, intent(out) , dimension(:,:,:) :: f

integer :: k
!-----------------------------------------------------------------------

f(:,:,1) = dt_xi(:,:,1)/b(:,:,1)

do k = 2, size(b,3)-1
f(:,:,k) = (dt_xi(:,:,k) - c(:,:,k)*f(:,:,k-1))*g(:,:,k)
enddo

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

end subroutine compute_f

subroutine explicit_tend (mu, nu, xi, dt_xi)

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

real, intent(in) , dimension(:,:,:) :: mu, nu, xi
real, intent(inout) , dimension(:,:,:) :: dt_xi

real, dimension(size(xi,1),size(xi,2),size(xi,3)) :: fluxx

integer :: nlev

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

nlev = size(mu,3)

fluxx(:,:,1) = 0.0
fluxx(:,:,2:nlev) = nu(:,:,2:nlev)*(xi(:,:,2:nlev) - xi(:,:,1:nlev-1))

dt_xi(:,:,1:nlev-1) = dt_xi(:,:,1:nlev-1) + &
mu(:,:,1:nlev-1)*(fluxx(:,:,2:nlev) - fluxx(:,:,1:nlev-1))
dt_xi(:,:,nlev) = dt_xi(:,:,nlev) - mu(:,:,nlev)*fluxx(:,:,nlev)

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

end subroutine explicit_tend

subroutine vert_diff_down_2 &
(delt, mu, nu, xi_1, xi_2, dt_xi_1, dt_xi_2, e, f_1, f_2, &
mu_delt_n, nu_n, e_n1, f_1_delt_n1, f_2_delt_n1, &
delta_1_n, delta_2_n, kbot)

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

real, intent(in) :: delt
real, intent(in) , dimension(:,:,:) :: mu, nu, xi_1, xi_2
real, intent(in) , dimension(:,:,:) :: dt_xi_1, dt_xi_2
real, intent(out) , dimension(:,:,:) :: e, f_1, f_2
real, intent(out) , dimension(:,:) :: mu_delt_n, nu_n, e_n1, &
f_1_delt_n1, f_2_delt_n1, &
delta_1_n, delta_2_n

integer, intent(in), dimension(:,:), optional :: kbot

real, dimension(size(xi_1,1),size(xi_1,2),size(xi_1,3)) :: a, b, c, g, &
dt_xi_11, dt_xi_22

integer :: i, j, kb, nlev

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

! local copy of input
dt_xi_11 = dt_xi_1
dt_xi_22 = dt_xi_2

call explicit_tend (mu, nu, xi_1, dt_xi_11)
call explicit_tend (mu, nu, xi_2, dt_xi_22)

call compute_e (delt, mu, nu, e, a, b, c, g)

call compute_f (dt_xi_11, b, c, g, f_1)
call compute_f (dt_xi_22, b, c, g, f_2)

if (present(kbot)) then
do j=1,size(xi_1,2)
do i=1,size(xi_1,1)
kb = kbot(i,j)
mu_delt_n(i,j) = mu(i,j,kb )*delt
nu_n(i,j) = nu(i,j,kb )
e_n1(i,j) = e(i,j,kb-1)
f_1_delt_n1(i,j) = f_1(i,j,kb-1)*delt
f_2_delt_n1(i,j) = f_2(i,j,kb-1)*delt
delta_1_n(i,j) = dt_xi_11(i,j,kb )*delt
delta_2_n(i,j) = dt_xi_22(i,j,kb )*delt
enddo
enddo
else
nlev = size(mu,3)
mu_delt_n(:,:) = mu(:,:,nlev )*delt
nu_n(:,:) = nu(:,:,nlev )
e_n1(:,:) = e(:,:,nlev-1)
f_1_delt_n1(:,:) = f_1(:,:,nlev-1)*delt
f_2_delt_n1(:,:) = f_2(:,:,nlev-1)*delt
delta_1_n(:,:) = dt_xi_11(:,:,nlev )*delt
delta_2_n(:,:) = dt_xi_22(:,:,nlev )*delt
endif

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

end subroutine vert_diff_down_2

end module
52 changes: 46 additions & 6 deletions gsmphys/satmedmfvdifq_down_up.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,12 @@ subroutine satmedmfvdifq_down(ix,im,km,ntrac,ntcw,ntiw,ntke,
& ri(im,km-1), tkmnz(im,km-1),
& rdzt(im,km-1),rlmnz(im,km),
& al(im,km-1), ad(im,km), au(im,km-1),
& f1(im,km), f2(im,km*(ntrac-1))
& f1(im,km), f2(im,km*(ntrac-1)),
& al_new(im,km-1), ad_new (im,km),
& au_new(im,km-1),
& f1_new(im,km), f2_new(im,km*(ntrac-1))
!
!
!
real(kind=kind_phys) elm(im,km), ele(im,km),
& ckz(im,km), chz(im,km), frik(im),
Expand Down Expand Up @@ -1604,6 +1609,26 @@ subroutine satmedmfvdifq_down(ix,im,km,ntrac,ntcw,ntiw,ntke,
endif


! flip the matrix upside down
! since k=1 is surface
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do k=1,km-1
al_new(:, k) = au(:, km-1-k+1)
au_new(:, k) = al(:, km-1-k+1)
enddo

do k=1,km
ad_new(:, k) = ad(:, km-k+1)
f1_new(:, k) = f1(:, km-k+1)
enddo

do n = 1, ntrac-1
is = (n-1)*km
do k =1,km
f2_new(:, is+k) = f2(:, is+km-k+1)
enddo
end do

c
c solve tridiagonal problem for heat and moisture
Expand All @@ -1612,12 +1637,13 @@ subroutine satmedmfvdifq_down(ix,im,km,ntrac,ntcw,ntiw,ntke,
c call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)

c Modified: downward sweep only
call tridin_down(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)
call tridin_down(im,km,ntrac1,al_new,ad_new,
& au_new,f1_new,f2_new,au_new,f1_new,f2_new)

c Save upper diag and RHS for upward sweep
au_out(:,:)=au(:,:)
f1_out(:,:)=f1(:,:)
f2_out(:,:)=f2(:,:)
au_out(:,:)=au_new(:,:)
f1_out(:,:)=f1_new(:,:)
f2_out(:,:)=f2_new(:,:)
diss_out(:,:)=diss(:,:)

! to be rearranged in the upward sweep
Expand Down Expand Up @@ -1773,9 +1799,10 @@ subroutine satmedmfvdifq_up(ix,im,km,ntrac,ntke,
! Input arrays from part 1
real(kind=kind_phys) au_in(im,km-1), diss_in(im,km-1)
real(kind=kind_phys) f1_in(im,km), f2_in(im,km*(ntrac-1))
real(kind=kind_phys) f1_in_n(im,km), f2_in_n(im,km*(ntrac-1))
!
! Local variables
integer i, k, kk, is, ntrac1
integer i, n, k, kk, is, ntrac1
real(kind=kind_phys) rdt, ttend, qtend
real(kind=kind_phys) cont, conq
!
Expand All @@ -1789,6 +1816,19 @@ subroutine satmedmfvdifq_up(ix,im,km,ntrac,ntke,
! Perform upward sweep for heat, moisture and tracers
call tridin_up(im,km,ntrac1,au_in,f1_in,f2_in)

do k=1,km
f1_in_n(:,k)=f1_in(:,km-k+1)
enddo
f1_in=f1_in_n

do n = 1, ntrac-1
is = (n-1)*km
do k =1,km
f2_in_n(:, is+k) = f2_in(:, is+km-k+1)
enddo
end do
f2_in=f2_in_n

! Apply tendencies for heat and moisture
do k = 1,km
do i = 1,im
Expand Down