From 7dc2bfb1e375d7454657cf9e94a2783139b60845 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me=20Chanut?= <jerome.chanut@mercator-ocean.fr> Date: Mon, 5 Dec 2022 15:52:11 +0000 Subject: [PATCH] Resolve "FCT is not monotonic with non-linear free surface" --- src/OCE/TRA/traadv.F90 | 12 +++++----- src/OCE/TRA/traadv_fct.F90 | 48 +++++++++++++++++++------------------- src/OCE/step.F90 | 2 +- src/OCE/stpmlf.F90 | 2 +- src/OCE/stprk3_stg.F90 | 6 ++--- src/TOP/TRP/trcadv.F90 | 6 ++--- src/TOP/TRP/trctrp.F90 | 2 +- 7 files changed, 39 insertions(+), 39 deletions(-) diff --git a/src/OCE/TRA/traadv.F90 b/src/OCE/TRA/traadv.F90 index a55ffff5..d8bd85cb 100644 --- a/src/OCE/TRA/traadv.F90 +++ b/src/OCE/TRA/traadv.F90 @@ -78,7 +78,7 @@ MODULE traadv !!---------------------------------------------------------------------- CONTAINS - SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs, pau, pav, paw ) + SUBROUTINE tra_adv( kt, Kbb, Kmm, Kaa, pts, Krhs, pau, pav, paw ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv *** !! @@ -86,10 +86,10 @@ CONTAINS !! !! ** Method : - Update ts(Krhs) with the advective trend following nadv !!---------------------------------------------------------------------- - INTEGER , INTENT(in ) :: kt ! ocean time-step index - INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices - REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity - REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt) , INTENT(inout) :: pts ! active tracers and RHS of tracer equation + INTEGER , INTENT(in ) :: kt ! ocean time-step index + INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa, Krhs ! time level indices + REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity + REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt) , INTENT(inout) :: pts ! active tracers and RHS of tracer equation ! INTEGER :: ji, jj, jk ! dummy loop index REAL(wp), DIMENSION(:,:,:), POINTER :: zptu, zptv, zptw @@ -193,7 +193,7 @@ CONTAINS CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order - CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) + CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, Kaa, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) CASE ( np_MUS ) ! MUSCL CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) CASE ( np_UBS ) ! UBS diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90 index e66936ad..541e61e0 100644 --- a/src/OCE/TRA/traadv_fct.F90 +++ b/src/OCE/TRA/traadv_fct.F90 @@ -54,7 +54,7 @@ MODULE traadv_fct CONTAINS SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW, & - & Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) + & Kbb, Kmm, Kaa, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_fct *** !! @@ -70,8 +70,8 @@ CONTAINS !! - send trends to trdtra module for further diagnostics (l_trdtra=T) !! - poleward advective heat and salt transport (ln_diaptr=T) !!---------------------------------------------------------------------- - INTEGER , INTENT(in ) :: kt ! ocean time-step index - INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices + INTEGER , INTENT(in ) :: kt ! ocean time-step index + INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa, Krhs ! ocean time level indices INTEGER , INTENT(in ) :: kit000 ! first time step index CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) INTEGER , INTENT(in ) :: kjpt ! number of tracers @@ -93,7 +93,7 @@ CONTAINS !!---------------------------------------------------------------------- ! #if defined key_loop_fusion - CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) + CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, Kaa, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) #else IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( kt == kit000 ) THEN @@ -143,9 +143,9 @@ CONTAINS ALLOCATE(zwdia(A2D(nn_hls),jpk), zwinf(A2D(nn_hls),jpk), zwsup(A2D(nn_hls),jpk)) DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & - & / e3t(ji,jj,jk,Krhs) - zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) - zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) + & / e3t(ji,jj,jk,Kaa) + zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Kaa) + zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Kaa) END_3D END IF ! @@ -189,7 +189,7 @@ CONTAINS pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & - & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + & / e3t(ji,jj,jk,Kaa) * tmask(ji,jj,jk) END_3D IF ( ll_zAimp ) THEN @@ -309,7 +309,7 @@ CONTAINS ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) - ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Kaa) * tmask(ji,jj,jk) END_3D ! CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) @@ -323,7 +323,7 @@ CONTAINS ! ! !== monotonicity algorithm ==! ! - CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt ) + CALL nonosc( Kaa, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt ) ! ! !== final trend with corrected fluxes ==! ! @@ -332,7 +332,7 @@ CONTAINS & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) - zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Kaa) * tmask(ji,jj,jk) END_3D ! IF ( ll_zAimp ) THEN @@ -384,7 +384,7 @@ CONTAINS END SUBROUTINE tra_adv_fct - SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) + SUBROUTINE nonosc( Kaa, pbef, paa, pbb, pcc, paft, p2dt ) !!--------------------------------------------------------------------- !! *** ROUTINE nonosc *** !! @@ -397,7 +397,7 @@ CONTAINS !! drange (1995) multi-dimensional forward-in-time and upstream- !! in-space based differencing for fluid !!---------------------------------------------------------------------- - INTEGER , INTENT(in ) :: Kmm ! time level index + INTEGER , INTENT(in ) :: Kaa ! time level index REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pbef ! before field REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(in ) :: paft ! after field @@ -451,7 +451,7 @@ CONTAINS & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) ! up & down beta terms - zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt + zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kaa) / p2dt zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt END_2D @@ -689,7 +689,7 @@ CONTAINS out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji,jj+1,jk,jn,Kbb) ) SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW, & - & Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) + & Kbb, Kmm, Kaa, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_fct *** !! @@ -705,8 +705,8 @@ CONTAINS !! - send trends to trdtra module for further diagnostics (l_trdtra=T) !! - poleward advective heat and salt transport (ln_diaptr=T) !!---------------------------------------------------------------------- - INTEGER , INTENT(in ) :: kt ! ocean time-step index - INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices + INTEGER , INTENT(in ) :: kt ! ocean time-step index + INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa, Krhs ! ocean time level indices INTEGER , INTENT(in ) :: kit000 ! first time step index CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) INTEGER , INTENT(in ) :: kjpt ! number of tracers @@ -766,9 +766,9 @@ CONTAINS ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) DO_3D( 1, 1, 1, 1, 1, jpkm1 ) zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & - & / e3t(ji,jj,jk,Krhs) - zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) - zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) + & / e3t(ji,jj,jk,Kaa) + zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Kaa) + zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Kaa) END_3D END IF ! @@ -813,7 +813,7 @@ CONTAINS pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & - & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + & / e3t(ji,jj,jk,Kaa ) * tmask(ji,jj,jk) END_2D END DO @@ -929,7 +929,7 @@ CONTAINS ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) & & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) - ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Kaa) * tmask(ji,jj,jk) END_3D ! CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) @@ -943,7 +943,7 @@ CONTAINS ! ! !== monotonicity algorithm ==! ! - CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt ) + CALL nonosc( Kaa, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt ) ! ! !== final trend with corrected fluxes ==! ! @@ -952,7 +952,7 @@ CONTAINS & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) - zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) + zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Kaa) * tmask(ji,jj,jk) END_3D ! IF ( ll_zAimp ) THEN diff --git a/src/OCE/step.F90 b/src/OCE/step.F90 index ff2759f8..f2f79061 100644 --- a/src/OCE/step.F90 +++ b/src/OCE/step.F90 @@ -336,7 +336,7 @@ CONTAINS DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) - CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS + CALL tra_adv ( kstp, Nbb, Nnn, Naa, ts, Nrhs ) ! hor. + vert. advection ==> RHS IF( ln_zdfmfc ) CALL tra_mfc ( kstp, Nbb, ts, Nrhs ) ! Mass Flux Convection IF( ln_zdfosm ) THEN CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90 index bae69670..32eb6c8d 100644 --- a/src/OCE/stpmlf.F90 +++ b/src/OCE/stpmlf.F90 @@ -351,7 +351,7 @@ CONTAINS DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) - CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS + CALL tra_adv ( kstp, Nbb, Nnn, Naa, ts, Nrhs ) ! hor. + vert. advection ==> RHS IF( ln_zdfmfc ) CALL tra_mfc ( kstp, Nbb, ts, Nrhs ) ! Mass Flux Convection IF( ln_zdfosm ) THEN CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS diff --git a/src/OCE/stprk3_stg.F90 b/src/OCE/stprk3_stg.F90 index b12a1beb..e413a253 100644 --- a/src/OCE/stprk3_stg.F90 +++ b/src/OCE/stprk3_stg.F90 @@ -271,7 +271,7 @@ CONTAINS ! CALL trc_sbc_RK3( kstp, Kmm, tr, Krhs, kstg ) ! surface boundary condition ! - CALL trc_adv ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection + CALL trc_adv ( kstp, Kbb, Kmm, Kaa, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection ! ! !== time integration ==! ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2) DO jn = 1, jptra @@ -299,7 +299,7 @@ CONTAINS ! CALL trc_sbc_RK3( kstp, Kmm, tr, Krhs, kstg ) ! surface boundary condition ! - CALL trc_adv ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection + CALL trc_adv ( kstp, Kbb, Kmm, Kaa, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection ! CALL trc_sms ( kstp, Kbb, Kbb, Krhs ) ! tracers: sinks and sources CALL trc_trp ( kstp, Kbb, Kmm, Krhs, Kaa ) ! transport of passive tracers (without advection) @@ -320,7 +320,7 @@ CONTAINS CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in potential density for tra_mle computation !===>>> CAUTION here may be without GM velocity but stokes drift required ! 0 barotropic divergence for GM != 0 barotropic divergence for SD !!st consistence 2D / 3D - flux de masse - CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww ) ! hor. + vert. advection ==> RHS + CALL tra_adv( kstp, Kbb, Kmm, Kaa, ts, Krhs, zaU, zaV, ww ) ! hor. + vert. advection ==> RHS !===>>>>>> stg1&2: Verify the necessity of these trends (we may need it as there are in the RHS of dynspg_ts ?) !!gm ====>>>> needed for heat and salt fluxes associated with mass/volume flux diff --git a/src/TOP/TRP/trcadv.F90 b/src/TOP/TRP/trcadv.F90 index 86776d87..92bee004 100644 --- a/src/TOP/TRP/trcadv.F90 +++ b/src/TOP/TRP/trcadv.F90 @@ -72,7 +72,7 @@ MODULE trcadv !!---------------------------------------------------------------------- CONTAINS - SUBROUTINE trc_adv( kt, Kbb, Kmm, ptr, Krhs, pau, pav, paw ) + SUBROUTINE trc_adv( kt, Kbb, Kmm, Kaa, ptr, Krhs, pau, pav, paw ) !!---------------------------------------------------------------------- !! *** ROUTINE trc_adv *** !! @@ -81,7 +81,7 @@ CONTAINS !! ** Method : - Update tr(Krhs) with the advective trend following nadv !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index - INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices + INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa, Krhs ! time level indices REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt) , INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation ! @@ -158,7 +158,7 @@ CONTAINS CASE ( np_CEN ) ! Centered : 2nd / 4th order CALL tra_adv_cen ( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v ) CASE ( np_FCT ) ! FCT : 2nd / 4th order - CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v ) + CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, Kaa, ptr, jptra, Krhs, nn_fct_h, nn_fct_v ) CASE ( np_MUS ) ! MUSCL CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups ) CASE ( np_UBS ) ! UBS diff --git a/src/TOP/TRP/trctrp.F90 b/src/TOP/TRP/trctrp.F90 index 135466dc..6996ccaa 100644 --- a/src/TOP/TRP/trctrp.F90 +++ b/src/TOP/TRP/trctrp.F90 @@ -89,7 +89,7 @@ CONTAINS #endif #if ! defined key_RK3 ! ! MLF only: add the advection trend to the RHS - CALL trc_adv ( kt, Kbb, Kmm, tr, Krhs ) ! horizontal & vertical advection + CALL trc_adv ( kt, Kbb, Kmm, Kaa, tr, Krhs ) ! horizontal & vertical advection #endif CALL trc_ldf ( kt, Kbb, Kmm, tr, Krhs ) ! lateral mixing CALL trc_zdf ( kt, Kbb, Kmm, Krhs, tr, Kaa ) ! vert. mixing & after tracer ==> after -- GitLab