diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 97c8cc38..c382ff6c 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -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(:,:) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index ea6d8a37..d43152be 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -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(:,:,:) @@ -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)) @@ -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 diff --git a/atmos_shared/surf_diff.F90 b/atmos_shared/surf_diff.F90 new file mode 100644 index 00000000..4b2f62fd --- /dev/null +++ b/atmos_shared/surf_diff.F90 @@ -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 diff --git a/gsmphys/satmedmfvdifq_down_up.f b/gsmphys/satmedmfvdifq_down_up.f index 20cba263..2a35ae53 100644 --- a/gsmphys/satmedmfvdifq_down_up.f +++ b/gsmphys/satmedmfvdifq_down_up.f @@ -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), @@ -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 @@ -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 @@ -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 ! @@ -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