Skip to content
Snippets Groups Projects
Commit 5c5fa98c authored by Sibylle TECHENE's avatar Sibylle TECHENE
Browse files

merge branch 80 : save memory by reducing temporary arrays dimension for adv... 

merge branch 80 : save memory by reducing temporary arrays dimension for adv and zdf routines for tra and dyn and dynvor ( it passes sette with AtFC options but changes results )
parents 05ac2727 ee6ce8c5
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1970 additions and 573 deletions
src/OCE/DYN/dynadv.F90
+ 18
− 7
View file @ 5c5fa98c
......@@ -7,6 +7,7 @@ MODULE dynadv
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.0 ! 2017-07 (G. Madec) add a linear dynamics option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -50,7 +51,7 @@ MODULE dynadv
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
......@@ -64,7 +65,6 @@ CONTAINS
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum Eq.
!!----------------------------------------------------------------------
......@@ -73,12 +73,23 @@ CONTAINS
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 ) != vector form =!
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! horizontal gradient of kinetic energy
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) ! vertical advection
! !* horizontal gradient of kinetic energy
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_keg_hls1( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! lbc needed with Hollingsworth scheme
ELSE ! halo 2 case
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs )
ENDIF
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) !* vertical advection
!
CASE( np_FLX_c2 ) != flux form =!
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 2nd order centered scheme
CASE( np_FLX_ubs )
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw ) !* 2nd order centered scheme
!
CASE( np_FLX_ubs ) !* 3rd order UBS scheme (UP3)
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_adv_ubs_hls1( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ELSE ! halo 2 case
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ENDIF
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
+ 75
− 69
View file @ 5c5fa98c
......@@ -6,6 +6,7 @@ MODULE dynadv_cen2
!!======================================================================
!! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -35,7 +36,7 @@ MODULE dynadv_cen2
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
......@@ -51,15 +52,17 @@ CONTAINS
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kmm, Krhs ! ocean time-step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection computation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zzv ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
REAL(wp) :: zzu, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(1)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(1)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(1)) :: zfu_uw, zfv_vw, zfw
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -71,8 +74,9 @@ CONTAINS
ENDIF
!
IF( l_trddyn ) THEN ! trends: store the input trends
zfu_uw(:,:,:) = puu(:,:,:,Krhs)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
......@@ -89,84 +93,86 @@ CONTAINS
!
DO jk = 1, jpkm1 ! horizontal transport
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
zfu(ji,jj) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point)
zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ) = ( zfv(ji,jj) + zfv(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ) = ( zfu(ji,jj) + zfu(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
& + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
& + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
zfu_t(:,:,:) = puu(:,:,:,Krhs)
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( no_zad ) ) THEN !== No vertical advection ==! (except if linear free surface)
! ==
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,1,Kmm)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,1,Kmm)
END_2D
ENDIF
!
ELSE !== Vertical advection ==!
! !== Vertical advection ==!
!
! ! surface vertical fluxes
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO jk = 1, jpk-2 ! divergence of advective fluxes
!
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior advective fluxes
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
END_2D
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
END_2D
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
DO_2D( 0, 0, 0, 0 )
! ! vertical flux at level k+1
zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
! ! divergence of vertical momentum flux
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_kp1
zfv_vw(ji,jj) = zzfv_kp1
END_2D
END DO
!
jk = jpkm1
DO_2D( 0, 0, 0, 0 )
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_cen2
......
src/OCE/DYN/dynadv_ubs.F90
+ 276
− 58
View file @ 5c5fa98c
This diff is collapsed.
src/OCE/DYN/dynkeg.F90
+ 123
− 4
View file @ 5c5fa98c
......@@ -6,7 +6,8 @@ MODULE dynkeg
!! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
!! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
!! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -27,7 +28,8 @@ MODULE dynkeg
IMPLICIT NONE
PRIVATE
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg_hls1 ! routine called by step module
INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme)
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
......@@ -44,6 +46,123 @@ MODULE dynkeg
CONTAINS
SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
!! ** Purpose : Compute the now momentum trend due to the horizontal
!! gradient of the horizontal kinetic energy and add it to the
!! general momentum trend.
!!
!! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that
!! conserve kinetic energy. Compute the now horizontal kinetic energy
!! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
!! * kscheme = nkeg_HW : Hollingsworth correction following
!! Arakawa (2001). The now horizontal kinetic energy is given by:
!! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 )
!! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ]
!!
!! Take its horizontal gradient and add it to the general momentum
!! trend.
!! u(rhs) = u(rhs) - 1/e1u di[ zhke ]
!! v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
!!
!! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
!! - send this trends to trd_dyn (l_trddyn=T) for post-processing
!!
!! ** References : Arakawa, A., International Geophysics 2001.
!! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: kscheme ! =0/1 type of KEG scheme
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zu, zv ! local scalars
REAL(wp), DIMENSION(:,: ) , ALLOCATABLE :: zhke
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_keg')
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! Save the input trends
ALLOCATE( zu_trd(A2D(0),jpk), zv_trd(A2D(0),jpk) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
SELECT CASE ( kscheme )
!
CASE ( nkeg_C2 ) !== Standard scheme ==!
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 ) !* Horizontal kinetic energy at T-point
zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm)
zv = pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm)
zhke(ji,jj) = 0.25_wp * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
CASE ( nkeg_HW ) !* Hollingsworth scheme
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zu = ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) &
& + ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) * ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zv = ( pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) &
& + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zhke(ji,jj) = r1_48 * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
END SELECT
!
IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
SUBROUTINE dyn_keg_hls1( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
......@@ -86,7 +205,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) 'dyn_keg_hls1 : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
......@@ -147,7 +266,7 @@ CONTAINS
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
END SUBROUTINE dyn_keg_hls1
!!======================================================================
END MODULE dynkeg
src/OCE/DYN/dynvor.F90
+ 432
− 18
View file @ 5c5fa98c
This diff is collapsed.
src/OCE/DYN/dynzdf.F90
+ 312
− 297
View file @ 5c5fa98c
This diff is collapsed.
src/OCE/TRA/traadv.F90
+ 11
− 2
View file @ 5c5fa98c
......@@ -10,6 +10,7 @@ MODULE traadv
!! - ! 2014-12 (G. Madec) suppression of cross land advection option
!! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling
!! 4.5 ! 2021-04 (G. Madec, S. Techene) add advective velocities as optional arguments
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -178,7 +179,7 @@ CONTAINS
!
!!gm ???
! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(0),:) ) ! diagnose the effective MSF
IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF
!!gm ???
!
......@@ -191,11 +192,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
!
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 )
IF( nn_hls == 1 ) THEN
CALL tra_adv_cen_hls1( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ELSE
CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ENDIF
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 )
CASE ( np_MUS ) ! MUSCL
IF( nn_hls == 1 ) THEN
CALL tra_adv_mus_hls1( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
ELSE
CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
END IF
CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST
......
src/OCE/TRA/traadv_cen.F90
+ 333
− 5
View file @ 5c5fa98c
This diff is collapsed.
src/OCE/TRA/traadv_mus.F90
+ 217
− 3
View file @ 5c5fa98c
......@@ -9,6 +9,7 @@ MODULE traadv_mus
!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
!! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed
!! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -34,7 +35,9 @@ MODULE traadv_mus
IMPLICIT NONE
PRIVATE
PUBLIC tra_adv_mus ! routine called by traadv.F90
PUBLIC tra_adv_mus ! routine called by traadv.F90
PUBLIC tra_adv_mus_hls1 ! routine called by traadv.F90
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
! ! and in closed seas (orca 2 and 1 configurations)
......@@ -55,6 +58,217 @@ MODULE traadv_mus
CONTAINS
SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
!!
!! ** Purpose : Compute the now trend due to total advection of tracers
!! using a MUSCL scheme (Monotone Upstream-centered Scheme for
!! Conservation Laws) and add it to the general tracer trend.
!!
!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
!! ld_msc_ups=T :
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (ln_diaptr=T)
!!
!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, 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
LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
INTEGER :: ik
REAL(wp) :: zu, z0u, zzslpx, zzwx, zw , zalpha ! local scalars
REAL(wp) :: zv, z0v, zzslpy, zzwy, z0w ! - -
REAL(wp) :: zdzt_kp2, zslpz_kp1, zfW_kp1
REAL(wp), DIMENSION(A2D(2)) :: zdxt, zslpx, zwx ! 2D workspace
REAL(wp), DIMENSION(A2D(2)) :: zdyt, zslpy, zwy ! - -
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups
IF(lwp) WRITE(numout,*) '~~~~~~~'
IF(lwp) WRITE(numout,*)
!
! Upstream / MUSCL scheme indicator
!
ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed
IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked)
DO jk = 1, jpkm1
xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed
& - rnfmsk(:,:) * rnfmsk_z(jk) * tmask(:,:,jk) ! =>0 near runoff mouths (& closed sea outflows)
END DO
ENDIF
!
ENDIF
!
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( l_diaptr .AND. cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( l_iom .AND. cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
!
DO jn = 1, kjpt !== loop over the tracers ==!
!
DO jk = 1, jpkm1
! !* Horizontal advective fluxes
!
! !-- first guess of the slopes
DO_2D( 2, 1, 2, 1 )
zdxt(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * umask(ji,jj,jk)
zdyt(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * vmask(ji,jj,jk)
END_2D
! !-- Slopes of tracer
DO_2D( 1, 1, 1, 1 )
! ! 1/2 Slopes at T-point (set to 0 if adjectent slopes are of opposite sign)
zzslpx = ( zdxt(ji,jj) + zdxt(ji-1,jj ) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdxt(ji,jj) * zdxt(ji-1,jj ) ) )
zzslpy = ( zdyt(ji,jj) + zdyt(ji ,jj-1) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdyt(ji,jj) * zdyt(ji ,jj-1) ) )
! ! Slopes limitation
zslpx(ji,jj) = SIGN( 1.0_wp, zzslpx ) * MIN( ABS( zzslpx ), &
& 2._wp*ABS( zdxt (ji-1,jj) ), &
& 2._wp*ABS( zdxt (ji ,jj) ) )
zslpy(ji,jj) = SIGN( 1.0_wp, zzslpy ) * MIN( ABS( zzslpy ), &
& 2._wp*ABS( zdyt (ji,jj-1) ), &
& 2._wp*ABS( zdyt (ji,jj ) ) )
END_2D
!!gm + !!st ticket ? comparaison pommes et carrottes ABS(zzslpx) et 2._wp*ABS( zdxt (ji-1,jj) )
!
DO_2D( 1, 0, 1, 0 ) !-- MUSCL horizontal advective fluxes
z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
zalpha = 0.5_wp - z0u
zu = z0u - 0.5_wp * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj)
zzwy = pt(ji ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji ,jj)
zwx(ji,jj) = pU(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
!
z0v = SIGN( 0.5_wp, pV(ji,jj,jk) )
zalpha = 0.5_wp - z0v
zv = z0v - 0.5_wp * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1)
zzwy = pt(ji,jj ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj )
zwy(ji,jj) = pV(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj) - zwx(ji-1,jj ) &
& + zwy(ji,jj) - zwy(ji ,jj-1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!!
!! ! ! trend diagnostics
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_xad, zwx(:,:), pU, pt(:,:,:,jn,Kbb) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_yad, zwy(:,:), pV, pt(:,:,:,jn,Kbb) )
!! END IF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, jk, 'adv', zwy(:,:) )
!! ! ! heat transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, jk, 'adv', zwx(:,:), zwy(:,:) )
!
! !* Vertical advective fluxes
!
#define zdzt_kp1 zdxt
#define zslpz zslpx
#define zfW zwx
zfW (A2D(0)) = 0._wp ! anciennement zwx at jk = 1
! ! anciennement zwx at jk = 2
DO_2D( 0, 0, 0, 0 )
zdzt_kp1(ji,jj) = tmask(ji,jj,2) * ( pt(ji,jj,1,jn,Kbb) - pt(ji,jj,2,jn,Kbb) )
END_2D
zslpz (A2D(0)) = 0._wp ! anciennement zslpx at jk = 1
!
IF( ln_linssh ) THEN !-- linear ssh : non zero top values
DO_2D( 0, 0, 0, 0 ) ! at the ocean surface
zfW(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) ! surface flux
END_2D
IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
DO_2D( 0, 0, 0, 0 ) ! update pt(Krhs) under the ice-shelf
ik = mikt(ji,jj) ! the flux at ik-1 is zero ( inside ice-shelf )
IF( ik > 1 ) THEN
pt(ji,jj,ik,jn,Krhs) = pt(ji,jj,ik,jn,Krhs) - pW(ji,jj,ik) * pt(ji,jj,ik,jn,Kbb) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
ENDIF
END_2D
ENDIF
ENDIF
!
! wmask usage for computing zw and zwk is needed in isf case and linear ssh
!
!
DO jk = 1, jpkm1
IF( jk < jpkm1 ) THEN
DO_2D( 0, 0, 0, 0 )
! !-- Slopes of tracer
! ! masked vertical gradient at jk+2
zdzt_kp2 = ( pt(ji,jj,jk+1,jn,Kbb) - pt(ji,jj,jk+2,jn,Kbb) ) * tmask(ji,jj,jk+2) !!st wmask(ji,jj,jk+2)
! ! vertical slope at jk+1
zslpz_kp1 = ( zdzt_kp1(ji,jj) + zdzt_kp2 ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdzt_kp1(ji,jj) * zdzt_kp2 ) )
! ! slope limitation
zslpz_kp1 = SIGN( 1.0_wp, zslpz_kp1 ) * MIN( ABS( zslpz_kp1 ), &
& 2.*ABS( zdzt_kp2 ), &
& 2.*ABS( zdzt_kp1(ji,jj) ) )
! !-- vertical advective flux at jk+1
! ! caution: zfW_kp1 is masked for ice-shelf cavities
! ! since top fluxes already added to pt(Krhs) before the vertical loop
z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
zalpha = 0.5_wp + z0w
zw = z0w - 0.5_wp * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz_kp1
zzwy = pt(ji,jj,jk ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz(ji,jj)
zfW_kp1 = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)!!st * wmask(ji,jj,jk+1)
! !-- vertical advective trend at jk
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zfW(ji,jj) - zfW_kp1 ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
! ! updates for next level
zdzt_kp1(ji,jj) = zdzt_kp2
zslpz (ji,jj) = zslpz_kp1
zfW (ji,jj) = zfW_kp1
END_2D
ELSE
DO_2D( 0, 0, 0, 0 ) !-- vertical advective trend at jpkm1
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zfW(ji,jj) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
ENDIF
END DO ! end of jk loop
!
!!gm + !!st idem see above
!! ! ! send trends for diagnostic
!! IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
!
END DO ! end of tracer loop
!
END SUBROUTINE tra_adv_mus
SUBROUTINE tra_adv_mus_hls1( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
......@@ -178,7 +392,7 @@ CONTAINS
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
! ! trend diagnostics
......@@ -239,7 +453,7 @@ CONTAINS
!
END DO ! end of tracer loop
!
END SUBROUTINE tra_adv_mus
END SUBROUTINE tra_adv_mus_hls1
!!======================================================================
END MODULE traadv_mus
src/OCE/TRA/trazdf.F90
+ 156
− 106
View file @ 5c5fa98c
......@@ -6,6 +6,7 @@ MODULE trazdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA
!! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option
!! 4.5 ! 2022-06 (G. Madec) refactoring to reduce memory usage (j-k-i loops)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -22,7 +23,7 @@ MODULE trazdf
USE ldfslp ! lateral diffusion: iso-neutral slope
USE trd_oce ! trends: ocean variables
USE trdtra ! trends: tracer trend manager
USE eosbn2, ONLY: ln_SEOS, rn_b0
USE eosbn2 , ONLY: ln_SEOS, rn_b0
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
......@@ -77,7 +78,7 @@ CONTAINS
ENDIF
!
! !* compute lateral mixing trend and add it to the general trend
CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
CALL tra_zdf_imp( 'TRA', Kbb, Kmm, Krhs, pts, Kaa, jpts )
!!gm WHY here ! and I don't like that !
! DRAKKAR SSS control {
......@@ -116,7 +117,7 @@ CONTAINS
END SUBROUTINE tra_zdf
SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
SUBROUTINE tra_zdf_imp( cdtype, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_zdf_imp ***
!!
......@@ -136,128 +137,177 @@ CONTAINS
!!
!! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer
!!---------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! 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
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: zrhs, zzwi, zzws ! local scalars
REAL(wp), DIMENSION(A2D(0),jpk) :: zwi, zwt, zwd, zws
REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwt, zwd, zws
!!---------------------------------------------------------------------
!
! ! ============= !
DO jn = 1, kjpt ! tracer loop !
! ! ============= !
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
! ! ================= !
DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! ================= !
DO jn = 1, kjpt ! tracer loop !
! ! ================= !
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
DO_3D( 0, 0, 0, 0, 2, jpk )
zwt(ji,jj,jk) = avt(ji,jj,jk)
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 2, jpk )
zwt(ji,jj,jk) = avs(ji,jj,jk)
END_3D
ENDIF
zwt(:,:,1) = 0._wp
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
END_3D
ELSE ! standard or triad iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_3D
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
!
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ! use avt for temperature
!
IF( l_ldfslp ) THEN ! use avt + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE ! use avt only
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk)
END_2D
ENDIF
!
ELSE ! use avs for salinty or passive tracers
!
IF( l_ldfslp ) THEN ! use avs + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE !
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk)
END_2D
ENDIF
ENDIF
zwt(:,1) = 0._wp
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
& + rDt * ( MAX( wi(ji,jj,jk ) , 0._wp ) &
& - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jk) = zzwi + rDt * MIN( wi(ji,jj,jk ) , 0._wp )
zws(ji,jk) = zzws - rDt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_2D
ELSE
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zwi(ji,jk) = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jk) = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jk) - zws(ji,jk)
END_2D
ENDIF
!
!!gm BUG?? : if edmfm is equivalent to a w ==>>> just add +/- rDt * edmfm(ji,jj,jk+1/jk )
!! but edmfm is at t-point !!!! crazy??? why not keep it at w-point????
!
IF( ln_zdfmfc ) THEN ! add upward Mass Flux in the matrix
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zws(ji,jk) = zws(ji,jk) + e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = zwd(ji,jk) - e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
END_2D
ENDIF
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! edmfa(ji,jj,jk) = 0._wp
! edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
! edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
! END_3D
!!gm BUG : level jpk never used in the inversion
! DO_2D( 0, 0, 0, 0 )
! edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm)
! edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm)
! edmfc(ji,jj,jpk) = 0._wp
! END_2D
!!
!!gm BUG ??? below e3t_Kmm should be used ?
!! or even no multiplication by e3t unless there is a bug in wi calculation
!!
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!!gm edmfa = 0._wp except at jpk which is not used ==>> zdiagi update is useless !
! zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
! zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
! zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
! END_3D
!!gm CALL diag_mfc( zwi, zwd, zws, rDt, Kaa )
!!gm SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_1Di( 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,1) = zwd(ji,1)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
END_2D
!
ENDIF
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
& + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp )
zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
END_3D
IF( ln_zdfmfc ) THEN ! add Mass Flux to the RHS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + edmftra(ji,jj,jk,jn)
END_2D
!!gm CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
ENDIF
!
! Modification of diagonal to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
END IF
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,jj,1) = zwd(ji,jj,1)
DO_1Di( 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb ) &
& + rDt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb ) &
& + rDt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jk) / zwt(ji,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
END_3D
!
ENDIF
!
! Modification of rhs to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
END IF
!
DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) &
& + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) &
& + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
END_3D
DO_1Di( 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jpkm1) * tmask(ji,jj,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jk) * tmask(ji,jj,jk)
END_2D
! ! ================= !
END DO ! tracer loop !
! ! ================= !
END DO ! end tracer loop !
END_1D ! i-k slices loop !
! ! ================= !
END SUBROUTINE tra_zdf_imp
......
src/OCE/do_loop_substitute.h90
+ 6
− 2
View file @ 5c5fa98c
......@@ -58,7 +58,10 @@
!
#endif
#define DO_2D(L, R, B, T) DO jj = ntsj-(B), ntej+(T) ; DO ji = ntsi-(L), ntei+(R)
#define DO_1Di(L, R) DO ji = ntsi-(L), ntei+(R)
#define DO_1Dj(B, T) DO jj = ntsj-(B), ntej+(T)
#define DO_2Dik(L, R, ks, ke, ki) DO jk = ks, ke, ki ; DO_1Di(L, R)
#define DO_2D(L, R, B, T) DO_1Dj(B, T) ; DO_1Di(L, R)
#define DO_2D_OVR(L, R, B, T) DO_2D(L-(L+R)*nthl, R-(R+L)*nthr, B-(B+T)*nthb, T-(T+B)*ntht)
#define A1Di(H) ntsi-(H):ntei+(H)
#define A1Dj(H) ntsj-(H):ntej+(H)
......@@ -76,5 +79,6 @@
#define DO_3DS(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki ; DO_2D(L, R, B, T)
#define DO_3DS_OVR(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki ; DO_2D_OVR(L, R, B, T)
#define END_1D END DO
#define END_2D END DO ; END DO
#define END_3D END DO ; END DO ; END DO
#define END_3D END DO ; END DO ; END DO
\ No newline at end of file
src/TOP/TRP/trcadv.F90
+ 10
− 1
View file @ 5c5fa98c
......@@ -8,6 +8,7 @@ MODULE trcadv
!! 3.7 ! 2014-05 (G. Madec, C. Ethe) Add 2nd/4th order cases for CEN and FCT schemes
!! 4.0 ! 2017-09 (G. Madec) remove vertical time-splitting option
!! 4.5 ! 2021-08 (G. Madec, S. Techene) add advective velocities as optional arguments
!! - ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
#if defined key_top
!!----------------------------------------------------------------------
......@@ -156,11 +157,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
!
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 )
IF( nn_hls == 1 ) THEN
CALL tra_adv_cen_hls1( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
ELSE
CALL tra_adv_cen ( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
ENDIF
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 )
CASE ( np_MUS ) ! MUSCL
IF( nn_hls == 1 ) THEN
CALL tra_adv_mus_hls1( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups )
ELSE
CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups )
ENDIF
CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST
......
src/TOP/TRP/trczdf.F90
+ 1
− 1
View file @ 5c5fa98c
......@@ -53,7 +53,7 @@ CONTAINS
!
IF( l_trdtrc ) ztrtrd(:,:,:,:) = ptr(:,:,:,:,Krhs)
!
CALL tra_zdf_imp( kt, nittrc000, 'TRC', rDt_trc, Kbb, Kmm, Krhs, ptr, Kaa, jptra ) ! implicit scheme
CALL tra_zdf_imp( 'TRC', Kbb, Kmm, Krhs, ptr, Kaa, jptra ) ! implicit scheme
!
IF( l_trdtrc ) THEN ! save the vertical diffusive trends for further diagnostics
DO jn = 1, jptra
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment