From 019661d8e57f6b84b86bc203cf196d1a7c58f428 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 21 Nov 2019 16:03:44 -0900 Subject: [PATCH] Adding TS_HSIMT option --- ROMS/Include/cppdefs.h | 1 + ROMS/Nonlinear/hsimt_tvd.F | 296 +++++++++++++++++++++++++++++++ ROMS/Nonlinear/pre_step3d.F | 39 +++- ROMS/Nonlinear/step3d_t.F | 185 ++++++++++++++++++- ROMS/Representer/rp_pre_step3d.F | 6 +- ROMS/Utility/inp_par.F | 2 +- 6 files changed, 517 insertions(+), 12 deletions(-) create mode 100755 ROMS/Nonlinear/hsimt_tvd.F diff --git a/ROMS/Include/cppdefs.h b/ROMS/Include/cppdefs.h index b56d9e65d..70d6b8090 100644 --- a/ROMS/Include/cppdefs.h +++ b/ROMS/Include/cppdefs.h @@ -74,6 +74,7 @@ ** TS_C2HADVECTION use if 2nd-order centered horizontal advection ** ** TS_C4HADVECTION use if 4th-order centered horizontal advection ** ** TS_MPDATA use if recursive MPDATA 3D advection ** +** TS_HSIMT use if scheme by Wu and Zhu (2010) ** ** TS_U3HADVECTION use if 3rd-order upstream horiz. advection ** ** TS_A4VADVECTION use if 4th-order Akima vertical advection ** ** TS_C2VADVECTION use if 2nd-order centered vertical advection ** diff --git a/ROMS/Nonlinear/hsimt_tvd.F b/ROMS/Nonlinear/hsimt_tvd.F new file mode 100755 index 000000000..e6fadcf01 --- /dev/null +++ b/ROMS/Nonlinear/hsimt_tvd.F @@ -0,0 +1,296 @@ +#include "cppdefs.h" +#define HSIMT_HOT + + MODULE hsimt_tvd_mod +#if defined NONLINEAR && defined TS_HSIMT && defined SOLVE3D +! +!================================================== Hernan G. Arango === +! Copyright (c) 2002-2014 The ROMS/TOMS Group John C. Warner ! +! Licensed under a MIT/X style license ! +! See License_ROMS.txt ! +!======================================================================= +! ! +!This routine computes anti-diffusive tracer flux based on HSIMT-TVD ! +!by Wu and Zhu (2010). This routine is for personal test only currently! +! ! +! On Output: FX, FE ! +! ! +! ! +! Reference: ! +! ! +! Hui Wu and Jianrong Zhu (2010), Advection scheme with 3rd ! +! high-order spatial interpolation at the middle temporal level ! +! and its application to saltwater intrusion in the Changjiang ! +! Estuary, Ocean Modelling 33, 33-51. ! +! Please contact Hui Wu (hwusklec@gmail.com) if have any questions ! +! ! +!======================================================================= +! + implicit none + + PUBLIC :: hsimt_tvd_tile + + CONTAINS +! +!*********************************************************************** + SUBROUTINE hsimt_tvd_tile (ng, tile, & + & LBi, UBi, LBj, UBj, & + & IminS, ImaxS, JminS, JmaxS, & +# ifdef MASKING + & rmask, umask, vmask, & +# endif +# ifdef WET_DRY + & rmask_wet, umask_wet, vmask_wet, & +# endif + & pm, pn, omn, om_u, on_v,u_k,v_k, & + & z_r, & + & Huon_k, Hvom_k, t_k, & + & FX,FE) +!*********************************************************************** +! + USE mod_param + USE mod_ncparam + USE mod_scalars +! +! Imported variable declarations. +! + integer, intent(in) :: ng, tile + integer, intent(in) :: LBi, UBi, LBj, UBj + integer, intent(in) :: IminS, ImaxS, JminS, JmaxS +! +# ifdef ASSUMED_SHAPE +# ifdef MASKING + real(r8), intent(in) :: rmask(LBi:,LBj:) + real(r8), intent(in) :: umask(LBi:,LBj:) + real(r8), intent(in) :: vmask(LBi:,LBj:) +# endif +# ifdef WET_DRY + real(r8), intent(in) :: rmask_wet(LBi:,LBj:) + real(r8), intent(in) :: umask_wet(LBi:,LBj:) + real(r8), intent(in) :: vmask_wet(LBi:,LBj:) +# endif + real(r8), intent(in) :: pm(LBi:,LBj:) + real(r8), intent(in) :: pn(LBi:,LBj:) + real(r8), intent(in) :: omn(LBi:,LBj:) + real(r8), intent(in) :: om_u(LBi:,LBj:) + real(r8), intent(in) :: on_v(LBi:,LBj:) + real(r8), intent(in) :: z_r(LBi:,LBj:,:) + real(r8), intent(in) :: Huon_k(LBi:,LBj:) + real(r8), intent(in) :: Hvom_k(LBi:,LBj:) + real(r8), intent(in) :: t_k(LBi:,LBj:) + real(r8), intent(in) :: u_k(LBi:,LBj:) + real(r8), intent(in) :: v_k(LBi:,LBj:) + real(r8), intent(out) :: FE(IminS:ImaxS,JminS:JmaxS) + real(r8), intent(out) :: FX(IminS:ImaxS,JminS:JmaxS) +# else +# ifdef MASKING + real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) +# endif +# ifdef WET_DRY + real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj) +# endif + real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) + real(r8), intent(in) :: Huon_k(LBi:,LBj:) + real(r8), intent(in) :: Hvom_k(LBi:,LBj:) + real(r8), intent(in) :: t_k(LBi:,LBj:) + real(r8), intent(in) :: u_k(LBi:,LBj:) + real(r8), intent(in) :: v_k(LBi:,LBj:) + real(r8), intent(out) :: FE(IminS:ImaxS,JminS:JmaxS) + real(r8), intent(out) :: FX(IminS:ImaxS,JminS:JmaxS) +# endif +! +! Local variable declarations. +! + integer :: i, is, j, k + real(r8) :: cc1,cc2,cc3 + real(r8) :: sl,rl,rkal,a1,b1,betal,sr,rtr,rkar,betar + real(r8) :: sd,rd,rkad,betad,su,ru,rkau,betau,epson + real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: kax + real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: kay + +# include "set_bounds.h" + + + DO j=JstrV-1,Jendp1 + DO i=IstrU-1,Iendp2 + kax(i,j)=1.0_r8-abs(u_k(i,j)*dt(ng)/om_u(i,j))+1.e-10_r8 + END DO + END DO + DO j=JstrV-1,Jendp2 + DO i=IstrU-1,Iendp1 + kay(i,j)=1.0_r8-abs(v_k(i,j)*dt(ng)/on_v(i,j))+1.e-10_r8 + END DO + END DO + + cc1=0.25_r8 + cc2=0.5_r8 + cc3=1.0_r8/12.0_r8 + epson=1.e-10_r8 + DO j=JstrV-1,Jendp1 + DO i=IstrU-1,Iendp2 +# ifdef MASKING + if (umask(i,j).eq.1.0_r8) then +# endif + + if (u_k(i,j).ge.0._r8) then + + if (i.le.3) then + sl=t_k(i-1,j) + else +# ifdef MASKING + if (rmask(i-2,j).eq.0.0_r8) then + sl=t_k(i-1,j) + else + if (rmask(i,j).eq.0.0_r8) then + sl=t_k(i-1,j) + else +# endif + + if (abs(t_k(i,j)-t_k(i-1,j)).le.epson) then + rl=0.0_r8 + rkal=0.0_r8 + else + rl=(t_k(i-1,j)-t_k(i-2,j))/(t_k(i,j)-t_k(i-1,j)) + rkal=kax(i-1,j)/kax(i,j) + endif + A1=cc1*kax(i,j)+cc2-cc3/kax(i,j) + B1=-cc1*kax(i,j)+cc2+cc3/kax(i,j) + betal=A1+B1*rl + sl=t_k(i-1,j)+ & + & 0.5_r8*max(0.0_r8,min(2.0_r8,2.0_r8*rl*rkal,betal))* & + & (t_k(i,j)-t_k(i-1,j))*kax(i,j) +# ifdef MASKING + endif + endif +# endif + endif + FX(i,j)=sl*huon_k(i,j) + else + if (i.ge.Lm(ng)-2) then + sr=t_k(i,j) + else +# ifdef MASKING + if (rmask(i+1,j).eq.0.0_r8) then + sr=t_k(i,j) + else + if (rmask(i-1,j).eq.0.0_r8) then + sr=t_k(i,j) + else +# endif + if (abs(t_k(i,j)-t_k(i-1,j)).le.epson) then + rtr=0.0_r8 + rkar=0.0_r8 + else + rtr=(t_k(i+1,j)-t_k(i,j))/(t_k(i,j)-t_k(i-1,j)) + rkar=kax(i+1,j)/kax(i,j) + endif + A1=cc1*kax(i,j)+cc2-cc3/kax(i,j) + B1=-cc1*kax(i,j)+cc2+cc3/kax(i,j) + betar=a1+b1*rtr + sr=t_k(i,j)- & + & 0.5_r8*max(0.0_r8,min(2.0_r8,2.0_r8*rtr*rkar,betar))* & + & (t_k(i,j)-t_k(i-1,j))*kax(i,j) +# ifdef MASKING + endif + endif +# endif + endif + FX(i,j)=sr*huon_k(i,j) + endif +# ifdef MASKING + else + FX(i,j)=0.0_r8 + endif +# endif + enddo + enddo + + DO j=JstrV-1,Jendp2 + DO i=IstrU-1,Iendp1 +# ifdef MASKING + if (vmask(i,j).eq.1.0_r8) then +# endif + if (v_k(i,j).ge.0.0_r8) then + if (j.le.3) then + sd=t_k(i,j-1) + else +# ifdef MASKING + if (rmask(i,j-2).eq.0.0_r8) then + sd=t_k(i,j-1) + else + if (rmask(i,j).eq.0.0_r8) then + sd=t_k(i,j-1) + else +# endif + if (abs(t_k(i,j)-t_k(i,j-1)).le.epson) then + rd=0.0_r8 + rkad=0.0_r8 + else + rd=(t_k(i,j-1)-t_k(i,j-2))/(t_k(i,j)-t_k(i,j-1)) + rkad=kay(i,j-1)/kay(i,j) + endif + a1=cc1*kay(i,j)+cc2-cc3/kay(i,j) + b1=-cc1*kay(i,j)+cc2+cc3/kay(i,j) + betad=a1+b1*rd + sd=t_k(i,j-1)+ & + & 0.5_r8*max(0.0_r8,min(2.0_r8,2.0_r8*rd*rkad,betad))* & + & (t_k(i,j)-t_k(i,j-1))*kay(i,j) +# ifdef MASKING + endif + endif +# endif + endif + FE(i,j)=sd*hvom_k(i,j) + else + if (j.ge.Mm(ng)-2) then + su=t_k(i,j) + else +# ifdef MASKING + if (rmask(i,j+1).eq.0.0_r8) then + su=t_k(i,j) + else + if (rmask(i,j-1).eq.0.0_r8) then + su=t_k(i,j) + else +# endif + if (abs(t_k(i,j)-t_k(i,j-1)).le.epson) then + ru=0.0_r8 + rkau=0.0_r8 + else + ru=(t_k(i,j+1)-t_k(i,j))/(t_k(i,j)-t_k(i,j-1)) + rkau=kay(i,j+1)/kay(i,j) + endif + a1=cc1*kay(i,j)+cc2-cc3/kay(i,j) + b1=-cc1*kay(i,j)+cc2+cc3/kay(i,j) + betau=a1+b1*ru + su=t_k(i,j)- & + & 0.5*max(0.0_r8,min(2.0_r8,2.0_r8*ru*rkau,betau))* & + & (t_k(i,j)-t_k(i,j-1))*kay(i,j) +# ifdef MASKING + endif + endif +# endif + endif + FE(i,j)=su*hvom_k(i,j) + endif +# ifdef MASKING + else + FE(i,j)=0.0_r8 + endif +# endif + enddo + enddo + + RETURN + END SUBROUTINE hsimt_tvd_tile +#endif + END MODULE hsimt_tvd_mod diff --git a/ROMS/Nonlinear/pre_step3d.F b/ROMS/Nonlinear/pre_step3d.F index 8abcc48d7..fbc0268b1 100644 --- a/ROMS/Nonlinear/pre_step3d.F +++ b/ROMS/Nonlinear/pre_step3d.F @@ -82,6 +82,14 @@ SUBROUTINE pre_step3d (ng, tile) & GRID(ng) % Hvom, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & +# ifdef TS_HSIMT + & GRID(ng) % omn, & + & GRID(ng) % om_u, & + & GRID(ng) % om_v, & + & GRID(ng) % on_u, & + & GRID(ng) % on_v, & + +# endif & FORCES(ng) % btflx, & & FORCES(ng) % bustr, & & FORCES(ng) % bvstr, & @@ -138,6 +146,9 @@ SUBROUTINE pre_step3d_tile (ng, tile, & & pm, pn, & & Hz, Huon, Hvom, & & z_r, z_w, & +# ifdef TS_HSIMT + & omn,om_u,om_v,on_u,on_v, & +# endif & btflx, bustr, bvstr, & & stflx, sustr, svstr, & # ifdef SOLAR_SOURCE @@ -187,6 +198,9 @@ SUBROUTINE pre_step3d_tile (ng, tile, & # ifdef T_PASSIVE USE pt3dbc_mod, ONLY : pt3dbc_tile # endif +# ifdef TS_HSIMT + USE hsimt_tvd_mod +# endif ! ! Imported variable declarations. ! @@ -211,6 +225,13 @@ SUBROUTINE pre_step3d_tile (ng, tile, & real(r8), intent(in) :: Hvom(LBi:,LBj:,:) real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(in) :: z_w(LBi:,LBj:,0:) +# ifdef TS_HSIMT + real(r8), intent(in) :: omn(LBi:,LBj:) + real(r8), intent(in) :: om_u(LBi:,LBj:) + real(r8), intent(in) :: om_v(LBi:,LBj:) + real(r8), intent(in) :: on_u(LBi:,LBj:) + real(r8), intent(in) :: on_v(LBi:,LBj:) +# endif real(r8), intent(in) :: btflx(LBi:,LBj:,:) real(r8), intent(in) :: bustr(LBi:,LBj:) real(r8), intent(in) :: bvstr(LBi:,LBj:) @@ -272,6 +293,13 @@ SUBROUTINE pre_step3d_tile (ng, tile, & real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng)) +# ifdef TS_HSIMT + real(r8), intent(in) :: omn(LBi:,LBj:) + real(r8), intent(in) :: om_u(LBi:,LBj:) + real(r8), intent(in) :: om_v(LBi:,LBj:) + real(r8), intent(in) :: on_u(LBi:,LBj:) + real(r8), intent(in) :: on_v(LBi:,LBj:) +# endif real(r8), intent(in) :: btflx(LBi:UBi,LBj:UBj,NT(ng)) real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj) @@ -320,7 +348,7 @@ SUBROUTINE pre_step3d_tile (ng, tile, & # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_UV integer :: idiag # endif -# ifdef TS_MPDATA +# if defined TS_MPDATA || defined TS_HSIMT real(r8), parameter :: Gamma = 0.5_r8 # else real(r8), parameter :: Gamma = 1.0_r8/6.0_r8 @@ -341,6 +369,11 @@ SUBROUTINE pre_step3d_tile (ng, tile, & real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curv real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad +# ifdef TS_HSIMT + real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: kaz + real(r8):: rd,rkad,a1,b1,betad,sw,rru,rkau,betau,epson,cc1,cc2,cc3 +# endif + # include "set_bounds.h" # ifndef TS_FIXED @@ -424,7 +457,7 @@ SUBROUTINE pre_step3d_tile (ng, tile, & END DO END DO -# elif defined TS_MPDATA +# elif defined TS_MPDATA || defined TS_HSIMT ! ! First-order, upstream differences horizontal advective fluxes. ! @@ -806,7 +839,7 @@ SUBROUTINE pre_step3d_tile (ng, tile, & FC(i,N(ng))=0.0_r8 END DO -# elif defined TS_MPDATA +# elif defined TS_MPDATA || defined TS_HSIMT ! ! First_order, upstream differences vertical advective flux. ! diff --git a/ROMS/Nonlinear/step3d_t.F b/ROMS/Nonlinear/step3d_t.F index 8583905b0..2afb4f45e 100644 --- a/ROMS/Nonlinear/step3d_t.F +++ b/ROMS/Nonlinear/step3d_t.F @@ -86,6 +86,21 @@ SUBROUTINE step3d_t (ng, tile) & GRID(ng) % om_v, & & GRID(ng) % on_u, & & GRID(ng) % on_v, & +# endif +# ifdef TS_HSIMT +# ifdef WET_DRY + & GRID(ng) % rmask_wet, & + & GRID(ng) % umask_wet, & + & GRID(ng) % vmask_wet, & +# endif + & GRID(ng) % omn, & + & GRID(ng) % om_u, & + & GRID(ng) % om_v, & + & GRID(ng) % on_u, & + & GRID(ng) % on_v, & + & OCEAN(ng) % u, & + & OCEAN(ng) % v, & + # endif & GRID(ng) % pm, & & GRID(ng) % pn, & @@ -138,6 +153,12 @@ SUBROUTINE step3d_t_tile (ng, tile, & & rmask_wet, umask_wet, vmask_wet, & # endif & omn, om_u, om_v, on_u, on_v, & +# endif +# ifdef TS_HSIMT +# ifdef WET_DRY + & rmask_wet, umask_wet, vmask_wet, & +# endif + & omn, om_u, om_v, on_u, on_v,u,v, & # endif & pm, pn, & & Hz, Huon, Hvom, & @@ -196,6 +217,9 @@ SUBROUTINE step3d_t_tile (ng, tile, & # ifdef TS_MPDATA USE mpdata_adiff_mod # endif +# ifdef TS_HSIMT + USE hsimt_tvd_mod +# endif # ifdef NESTING USE nesting_mod, ONLY : bry_fluxes # endif @@ -237,6 +261,21 @@ SUBROUTINE step3d_t_tile (ng, tile, & real(r8), intent(in) :: om_v(LBi:,LBj:) real(r8), intent(in) :: on_u(LBi:,LBj:) real(r8), intent(in) :: on_v(LBi:,LBj:) +# endif +# ifdef TS_HSIMT +# ifdef WET_DRY + real(r8), intent(in) :: rmask_wet(LBi:,LBj:) + real(r8), intent(in) :: umask_wet(LBi:,LBj:) + real(r8), intent(in) :: vmask_wet(LBi:,LBj:) +# endif + real(r8), intent(in) :: omn(LBi:,LBj:) + real(r8), intent(in) :: om_u(LBi:,LBj:) + real(r8), intent(in) :: om_v(LBi:,LBj:) + real(r8), intent(in) :: on_u(LBi:,LBj:) + real(r8), intent(in) :: on_v(LBi:,LBj:) + real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2) + real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2) + # endif real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) @@ -300,6 +339,19 @@ SUBROUTINE step3d_t_tile (ng, tile, & real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj) real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj) # endif +# ifdef TS_HSIMT +# ifdef WET_DRY + real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj) +# endif + real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj) + real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj) +# endif + real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj) real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) @@ -370,6 +422,14 @@ SUBROUTINE step3d_t_tile (ng, tile, & real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: Va real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: Wa # endif + +# ifdef TS_HSIMT + real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: kaz + real(r8):: rd,rkad,a1,b1,betad,sw,ru,rkau,betau,epson,cc1,cc2,cc3 + real(r8):: u_k(LBi:UBi,LBj:UBj) + real(r8):: v_k(LBi:UBi,LBj:UBj) +# endif + # ifdef AGE_DISTRIBUTION real(r8) :: fac # endif @@ -393,7 +453,7 @@ SUBROUTINE step3d_t_tile (ng, tile, & ! ! Compute inverse thickness. ! -# ifdef TS_MPDATA +# if defined TS_MPDATA || defined TS_HSIMT # define I_RANGE Istrm2,Iendp2 # define J_RANGE Jstrm2,Jendp2 # else @@ -438,6 +498,29 @@ SUBROUTINE step3d_t_tile (ng, tile, & # endif # endif +# ifdef TS_HSIMT +! +! The HSIMT algorithm requires a three-point footprint, so exchange +! boundary data on t(:,:,:,nnew,:) so other processes computed earlier +! (horizontal diffusion, biology, or sediment) are accounted. +! + IF (EWperiodic(ng).or.NSperiodic(ng)) THEN + DO itrc=1,NT(ng) + CALL exchange_r3d_tile (ng, tile, & + & LBi, UBi, LBj, UBj, 1, N(ng), & + & t(:,:,:,nnew,itrc)) + END DO + END IF + +# ifdef DISTRIBUTE + CALL mp_exchange4d (ng, tile, iNLM, 1, & + & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), & + & NghostPoints, & + & EWperiodic(ng), NSperiodic(ng), & + & t(:,:,:,nnew,:)) +# endif +# endif + #ifdef OFFLINE_BIOLOGY T_LOOP : DO ibt=1,NBT itrc=idbio(ibt) @@ -467,6 +550,27 @@ SUBROUTINE step3d_t_tile (ng, tile, & & t(i,j ,k,3,itrc)) END DO END DO + +# elif defined TS_HSIMT +! Third-order HSIMT-TVD horizontal advective fluxes. + + u_k(:,:)=0.5_r8*(u(:,:,k,1)+u(:,:,k,2)) + v_k(:,:)=0.5_r8*(v(:,:,k,1)+v(:,:,k,2)) + CALL hsimt_tvd_tile (ng, tile, & + & LBi, UBi, LBj, UBj, & + & IminS, ImaxS, JminS, JmaxS, & +# ifdef MASKING + & rmask, umask, vmask, & +# endif +# ifdef WET_DRY + & rmask_wet, umask_wet, vmask_wet, & +# endif + & pm, pn, omn, om_u, on_v, u_k, & + & v_k,z_r, & + & Huon(:,:,k), Hvom(:,:,k), & + & t(:,:,k,3,itrc), & + & FX,FE) + # elif defined TS_MPDATA ! @@ -630,7 +734,7 @@ SUBROUTINE step3d_t_tile (ng, tile, & Isrc=SOURCES(ng)%Isrc(is) Jsrc=SOURCES(ng)%Jsrc(is) IF (INT(SOURCES(ng)%Dsrc(is)).eq.0) THEN -# ifdef TS_MPDATA +# if defined TS_MPDATA || defined TS_HSIMT IF (((IstrUm2.le.Isrc).and.(Isrc.le.Iendp3)).and. & & ((JstrVm2.le.Jsrc).and.(Jsrc.le.Jendp2i))) THEN # else @@ -661,7 +765,7 @@ SUBROUTINE step3d_t_tile (ng, tile, & END IF END IF ELSE IF (INT(SOURCES(ng)%Dsrc(is)).eq.1) THEN -# ifdef TS_MPDATA +# if defined TS_MPDATA || defined TS_HSIMT IF (((IstrUm2.le.Isrc).and.(Isrc.le.Iendp2i)).and. & & ((JstrVm2.le.Jsrc).and.(Jsrc.le.Jendp3))) THEN # else @@ -804,6 +908,22 @@ SUBROUTINE step3d_t_tile (ng, tile, & ! Time-step vertical advection term. !----------------------------------------------------------------------- ! +# ifdef TS_HSIMT + do k=1,N(ng)-1 + do j=J_RANGE + do i=I_RANGE + kaz(i,j,k)=1.-abs(w(i,j,k)*pm(i,j)*pn(i,j)*dt(ng)/ & + & (z_r(i,j,k+1)-z_r(i,j,k))) + enddo + enddo + enddo + cc1=0.25_r8 + cc2=0.5_r8 + cc3=1.0_r8/12.0_r8 + epson=1.e-10_r8 + +# endif + DO j=J_RANGE #ifdef OFFLINE_BIOLOGY DO ibt=1,NBT @@ -935,6 +1055,61 @@ SUBROUTINE step3d_t_tile (ng, tile, & # endif FC(i,N(ng))=0.0_r8 END DO + +# elif defined TS_HSIMT + DO i=I_RANGE + + DO k=1,N(ng)-1 + if (k.eq.1.and.W(i,j,k).ge.0.0_r8) THEN + FC(i,k)=W(i,j,k)*t(i,j,k,3,itrc) + else if (k.eq.N(ng)-1.and.W(i,j,k).lt.0.0_r8) THEN + FC(i,k)=W(i,j,k)*t(i,j,k+1,3,itrc) + else + if (W(i,j,k).ge.0) THEN + if (abs(t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)).le.epson) THEN + rd=0.0_r8 + rkad=0.0_r8 + else + rd=(t(i,j,k,3,itrc)-t(i,j,k-1,3,itrc))/ & + & (t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)) + rkad=kaz(i,j,k-1)/(kaz(i,j,k)+epson) + endif + a1=cc1*kaz(i,j,k)+cc2-cc3/(kaz(i,j,k)+epson) + b1=-cc1*kaz(i,j,k)+cc2+cc3/(kaz(i,j,k)+epson) + betad=a1+b1*rd + sw=t(i,j,k,3,itrc)+ & + & 0.5_r8*max(0.0_r8,min(2.0_r8,2.0_r8*rd*rkad,betad))* & + & (t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc))*kaz(i,j,k) + else + if (abs(t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)).le.epson) THEN + ru=0.0_r8 + rkau=0.0_r8 + else + ru=(t(i,j,k+2,3,itrc)-t(i,j,k+1,3,itrc))/ & + & (t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)) + rkau=kaz(i,j,k+1)/(kaz(i,j,k)+epson) + endif + a1=cc1*kaz(i,j,k)+cc2-cc3/(kaz(i,j,k)+epson) + b1=-cc1*kaz(i,j,k)+cc2+cc3/(kaz(i,j,k)+epson) + betau=a1+b1*ru + sw=t(i,j,k+1,3,itrc)- & + & 0.5_r8*max(0.0_r8,min(2.0_r8,2.0_r8*ru*rkau,betau))* & + & (t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc))*kaz(i,j,k) + + endif + FC(i,k)=W(i,j,k)*sw + endif + + END DO +# ifdef SED_MORPH + FC(i,0)=W(i,j,0)*t(i,j,1,3,itrc) +# else + FC(i,0)=0.0_r8 +# endif + FC(i,N(ng))=0.0_r8 + END DO + + # else ! ! Fourth-order, central differences vertical advective flux. @@ -996,7 +1171,7 @@ SUBROUTINE step3d_t_tile (ng, tile, & Isrc=SOURCES(ng)%Isrc(is) Jsrc=SOURCES(ng)%Jsrc(is) IF (LtracerSrc(itrc,ng).and. & -# ifdef TS_MPDATA +# if defined TS_MPDATA || defined TS_HSIMT & ((IstrUm2.le.Isrc).and.(Isrc.le.Iendp2i)).and. & # else & ((IstrR.le.Isrc).and.(Isrc.le.IendR)).and. & @@ -1496,7 +1671,7 @@ SUBROUTINE step3d_t_tile (ng, tile, & END DO END DO # endif -# ifndef TS_MPDATA +# if !defined TS_MPDATA && !defined TS_HSIMT # define SALT_LIMIT # endif # ifdef SALT_LIMIT diff --git a/ROMS/Representer/rp_pre_step3d.F b/ROMS/Representer/rp_pre_step3d.F index 219a19a63..98eef9a04 100644 --- a/ROMS/Representer/rp_pre_step3d.F +++ b/ROMS/Representer/rp_pre_step3d.F @@ -382,7 +382,7 @@ SUBROUTINE rp_pre_step3d_tile (ng, tile, & # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_UV integer :: idiag # endif -# ifdef TS_MPDATA_NOT_YET +# if defined TS_MPDATA_NOT_YET || defined TS_HSIMT real(r8), parameter :: Gamma = 0.5_r8 # else real(r8), parameter :: Gamma = 1.0_r8/6.0_r8 @@ -505,7 +505,7 @@ SUBROUTINE rp_pre_step3d_tile (ng, tile, & END DO END DO -# elif defined TS_MPDATA_NOT_YET +# elif defined TS_MPDATA_NOT_YET || defined TS_HSIMT ! ! First-order, upstream differences horizontal advective fluxes. ! @@ -1183,7 +1183,7 @@ SUBROUTINE rp_pre_step3d_tile (ng, tile, & tl_FC(i,N(ng))=0.0_r8 END DO -# elif defined TS_MPDATA_NOT_YET +# elif defined TS_MPDATA_NOT_YET || defined TS_HSIMT ! ! First_order, upstream differences vertical advective flux. ! diff --git a/ROMS/Utility/inp_par.F b/ROMS/Utility/inp_par.F index aa8a36807..8b63d13ac 100644 --- a/ROMS/Utility/inp_par.F +++ b/ROMS/Utility/inp_par.F @@ -178,7 +178,7 @@ SUBROUTINE inp_par (model) ! ! Determine the number of ghost-points in the halo region. ! -#if defined TS_MPDATA || defined UV_VIS4 +#if defined TS_MPDATA || defined UV_VIS4 || defined TS_HSIMT NghostPoints=3 #else NghostPoints=2