Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1205 additions and 597 deletions
src/OCE/DYN/dynadv.F90
View file @ 1398e060
......@@ -45,12 +45,12 @@ MODULE dynadv
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv.F90 14053 2020-12-03 13:48:38Z techene $
!! $Id: dynadv.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
......@@ -63,21 +63,22 @@ CONTAINS
!! it is the relative vorticity which is added to coriolis term
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, 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 , 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.
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start( 'dyn_adv' )
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 )
CALL dyn_keg ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs ) ! vector form : horizontal gradient of kinetic energy
CALL dyn_zad ( kt, Kmm, puu, pvv, Krhs ) ! vector form : vertical advection
CASE( np_FLX_c2 )
CALL dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs ) ! 2nd order centered scheme
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
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 ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
View file @ 1398e060
......@@ -30,29 +30,36 @@ MODULE dynadv_cen2
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv_cen2.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynadv_cen2.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! ** Purpose : Compute the momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! ** Method : Trend evaluated with a 2nd order centered scheme
!! using fields at Kmm time-level.
!! In RK3 time stepping case, the optional arguments (pau,pav,paw)
!! are present. They are used as advective velocity while
!! the advected velocity remains (puu,pvv).
!!
!! ** Action : (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the now vorticity term trend
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
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 , 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), 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, 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), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -68,12 +75,22 @@ CONTAINS
zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
zpt_u => pau(:,:,:)
zpt_v => pav(:,:,:)
zpt_w => paw(:,:,:)
ELSE ! MLF: advective velocity = (puu,pvv,ww)
zpt_u => puu(:,:,:,Kmm)
zpt_v => pvv(:,:,:,Kmm)
zpt_w => ww (:,:,: )
ENDIF
!
! !== Horizontal advection ==!
!
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) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
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)
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) )
......@@ -83,11 +100,11 @@ CONTAINS
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) &
& / e3u(ji,jj,jk,Kmm)
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * 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) &
& / e3v(ji,jj,jk,Kmm)
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
......@@ -99,42 +116,57 @@ CONTAINS
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
! !== Vertical advection ==!
!
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
END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(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) * ww(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) )
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 ==!
!
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
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) &
& / 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 )
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) &
& / 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' )
!
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
View file @ 1398e060
......@@ -36,12 +36,12 @@ MODULE dynadv_ubs
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv_ubs.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynadv_ubs.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
......@@ -64,20 +64,26 @@ CONTAINS
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
!! In RK3 time stepping case, the optional arguments
!! (pau,pav,paw) are present. They are used as advective velocity
!! while the advected velocity remains (puu,pvv).
!!
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, 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 , 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(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) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v, 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), DIMENSION(A2D(nn_hls),jpk,2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlv_vv, zlv_vu
REAL(wp), DIMENSION(:,:,:), POINTER :: zpt_u, zpt_v, zpt_w
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -102,13 +108,24 @@ CONTAINS
zfu_uw(:,:,:) = puu(:,:,:,Krhs)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
zpt_u => pau(:,:,:)
zpt_v => pav(:,:,:)
zpt_w => paw(:,:,:)
ELSE ! MLF: advective velocity = (puu,pvv,ww)
zpt_u => puu(:,:,:,Kmm)
zpt_v => pvv(:,:,:,Kmm)
zpt_w => ww (:,:,: )
ENDIF
!
! ! =========================== !
DO jk = 1, jpkm1 ! Laplacian of the velocity !
! ! =========================== !
! ! horizontal volume fluxes
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian
......@@ -157,8 +174,8 @@ CONTAINS
DO jk = 1, jpkm1 ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
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)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
......@@ -212,42 +229,62 @@ CONTAINS
! ! ==================== !
! ! Vertical advection !
! ! ==================== !
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
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior fluxes
DO_2D( 0, 1, 0, 1 )
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(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) &
& / 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 ! save the vertical advection trend 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 )
! ! ======================== !
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 !
! ! =================== !
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
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
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 fluxes
DO_2D( 0, 1, 0, 1 )
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) &
& / 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 ! save the vertical advection trend 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=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_ubs
......
src/OCE/DYN/dynhpg.F90
View file @ 1398e060
......@@ -83,7 +83,7 @@ MODULE dynhpg
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynhpg.F90 15529 2021-11-23 15:00:19Z techene $
!! $Id: dynhpg.F90 15157 2021-07-29 08:28:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......
src/OCE/DYN/dynspg.F90
View file @ 1398e060
......@@ -51,12 +51,12 @@ MODULE dynspg
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynspg.F90 14225 2020-12-19 14:58:39Z smasson $
!! $Id: dynspg.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_spg ***
!!
......@@ -78,10 +78,9 @@ CONTAINS
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: z2dt, zg_2, zintp, zgrho0r, zld ! local scalars
REAL(wp) :: zg_2, zintp, zgrho0r, zld ! local scalars
REAL(wp) , DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
......@@ -150,8 +149,8 @@ CONTAINS
!
IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head
DO_2D( 0, 0, 0, 0 )
zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e2v(ji,jj)
END_2D
ENDIF
!
......@@ -166,7 +165,7 @@ CONTAINS
!
SELECT CASE ( nspg ) !== surface pressure gradient computed and add to the general trend ==!
CASE ( np_EXP ) ; CALL dyn_spg_exp( kt, Kmm, puu, pvv, Krhs ) ! explicit
CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) ! time-splitting
CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
END SELECT
!
IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics
......
src/OCE/DYN/dynspg_ts.F90
View file @ 1398e060
This diff is collapsed.
src/OCE/DYN/dynvor.F90
View file @ 1398e060
......@@ -50,6 +50,10 @@ MODULE dynvor
IMPLICIT NONE
PRIVATE
INTERFACE dyn_vor
MODULE PROCEDURE dyn_vor_3D, dyn_vor_RK3
END INTERFACE
PUBLIC dyn_vor ! routine called by step.F90
PUBLIC dyn_vor_init ! routine called by nemogcm.F90
......@@ -73,8 +77,9 @@ MODULE dynvor
INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme
INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme
INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity
! ! associated indices:
! !: choice of calculated vorticity
INTEGER, PUBLIC :: ncor, nrvm, ntot ! Coriolis, relative vorticity, total vorticity
! ! associated indices:
INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary)
INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity
INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term
......@@ -98,12 +103,12 @@ MODULE dynvor
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynvor.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynvor.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!!
!! ** Purpose : compute the lateral ocean tracer physics.
......@@ -120,7 +125,7 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_vor')
IF( ln_timing ) CALL timing_start('dyn_vor_3D')
!
IF( l_trddyn ) THEN !== trend diagnostics case : split the added trend in two parts ==!
!
......@@ -208,9 +213,85 @@ CONTAINS
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_vor')
IF( ln_timing ) CALL timing_stop('dyn_vor_3D')
!
END SUBROUTINE dyn_vor_3D
SUBROUTINE dyn_vor_RK3( kt, Kmm, puu, pvv, Krhs, knoco )
!!----------------------------------------------------------------------
!!
!! ** Purpose : compute the lateral ocean tracer physics.
!!
!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
!! and planetary vorticity trends) and send them to trd_dyn
!! for futher diagnostics (l_trddyn=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
INTEGER , INTENT(in ) :: knoco ! specified vorticity trend (= np_MET or np_RVO)
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_vor_RK3')
!
! !== total vorticity trend added to the general trend ==!
!!st WARNING 22/02 !!!!!!!! stoke drift or not stoke drift ? => bar to do later !!!
!! stoke drift a garder pas vortex force a priori !!
!! ATTENTION déja appelé avec Kbb !!
!
SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==!
CASE( np_ENT ) !* energy conserving scheme (T-pts)
CALL vor_enT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
CALL vor_eeT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_ENS ) !* enstrophy conserving scheme
CALL vor_ens( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_MIX ) !* mixed ene-ens scheme
CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens)
IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force
CASE( np_EEN ) !* energy and enstrophy conserving scheme
CALL vor_een( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
END SELECT
!
! ! print sum trends (used for debugging)
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_vor_RK3')
!
END SUBROUTINE dyn_vor
END SUBROUTINE dyn_vor_RK3
SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......
src/OCE/DYN/dynzad.F90
View file @ 1398e060
......@@ -3,8 +3,9 @@ MODULE dynzad
!! *** MODULE dynzad ***
!! Ocean dynamics : vertical advection trend
!!======================================================================
!! History : OPA ! 1991-01 (G. Madec) Original code
!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
!! History : OPA ! 1991-01 (G. Madec) Original code
!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
!! 4.5 ! 2021-01 (S. Techene, G. Madec) memory optimization
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -32,7 +33,7 @@ MODULE dynzad
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynzad.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynzad.F90 14428 2021-02-10 18:12:36Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -53,15 +54,14 @@ CONTAINS
!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the vert. momentum adv. trends
!! - Send the trends to trddyn for diagnostics (l_trddyn=T)
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step inedx
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 , INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step & 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) :: zua, zva ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zww
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwuw, zwvw
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv
REAL(wp) :: zWf, zWfi, zzWfu, zzWdzU ! local scalars
REAL(wp) :: zWfj, zzWfv, zzWdzV ! - -
REAL(wp), DIMENSION(A2D(0)) :: zWdzU, zWdzV ! 2D inner workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_zad')
......@@ -72,44 +72,53 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
ztrdu(:,:,:) = puu(:,:,:,Krhs)
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical
IF( ln_vortex_force ) THEN ! vertical fluxes
DO_2D( 0, 1, 0, 1 )
zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
END_2D
ELSE
DO_2D( 0, 1, 0, 1 )
zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
END_2D
!
! !== vertical momentum advection ==! at u- and v-points
!
zWdzU(A2D(0)) = 0._wp ! set surface (jk=1) vertical advection to zero
zWdzV(A2D(0)) = 0._wp
!
DO_3D( 0, 0, 0, 0 , 1, jpk-2 ) != surface to jpk-2 vertical advection
! ! vertical transport at jk+1 uw/vw-level (x2): 2*mi/j[e1e2t*(We)]
IF( ln_vortex_force ) THEN ! We = ww+wsd
zWf = e1e2t(ji ,jj ) * ( ww(ji ,jj ,jk+1) + wsd(ji ,jj ,jk+1) )
zWfi = e1e2t(ji+1,jj ) * ( ww(ji+1,jj ,jk+1) + wsd(ji+1,jj ,jk+1) )
zWfj = e1e2t(ji ,jj+1) * ( ww(ji ,jj+1,jk+1) + wsd(ji ,jj+1,jk+1) )
ELSE ! We = ww
zWf = e1e2t(ji ,jj ) * ww(ji ,jj ,jk+1)
zWfi = e1e2t(ji+1,jj ) * ww(ji+1,jj ,jk+1)
zWfj = e1e2t(ji ,jj+1) * ww(ji ,jj+1,jk+1)
ENDIF
DO_2D( 0, 0, 0, 0 ) ! vertical momentum advection at w-point
zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
END_2D
END DO
zzWfu = zWfi + zWf ! averaging at uw- and vw-points (x2)
zzWfv = zWfj + zWf
! ! vertical advection at jk+1 uw-level (x4): zzWfu/v*dk+1[u/v]
zzWdzU = zzWfu * ( puu(ji,jj,jk,Kmm) - puu(ji,jj,jk+1,Kmm) )
zzWdzV = zzWfv * ( pvv(ji,jj,jk,Kmm) - pvv(ji,jj,jk+1,Kmm) )
!
! ! vertical advection at jk u/v-level
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * ( zWdzU(ji,jj) + zzWdzU )
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * ( zWdzV(ji,jj) + zzWdzV )
!
zWdzU(ji,jj) = zzWdzU ! save for next level computation
zWdzV(ji,jj) = zzWdzV
END_3D
!
! Surface and bottom advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zwuw(ji,jj, 1 ) = 0._wp
zwvw(ji,jj, 1 ) = 0._wp
zwuw(ji,jj,jpk) = 0._wp
zwvw(ji,jj,jpk) = 0._wp
jk = jpkm1
DO_2D( 0, 0, 0, 0 ) != bottom vertical advection at jpkm1
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * zWdzU(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * zWdzV(ji,jj)
END_2D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Vertical momentum advection at u- and v-points
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
......
src/OCE/DYN/dynzdf.F90
View file @ 1398e060
This diff is collapsed.
src/OCE/DYN/sshwzv.F90
View file @ 1398e060
......@@ -16,7 +16,9 @@ MODULE sshwzv
!!----------------------------------------------------------------------
!! ssh_nxt : after ssh
!! ssh_atf : time filter the ssh arrays
!! wzv : compute now vertical velocity
!! wzv : generic interface of vertical velocity calculation
!! wzv_MLF : MLF: compute NOW vertical velocity
!! wzv_RK3 : RK3: Compute a vertical velocity
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE isf_oce ! ice shelf
......@@ -43,6 +45,10 @@ MODULE sshwzv
IMPLICIT NONE
PRIVATE
! !! * Interface
INTERFACE wzv
MODULE PROCEDURE wzv_MLF, wzv_RK3
END INTERFACE
PUBLIC ssh_nxt ! called by step.F90
PUBLIC wzv ! called by step.F90
......@@ -54,7 +60,7 @@ MODULE sshwzv
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sshwzv.F90 15150 2021-07-27 10:38:24Z smasson $
!! $Id: sshwzv.F90 14618 2021-03-19 14:42:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -76,10 +82,11 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
!
!
INTEGER :: ji, jj, jk ! dummy loop index
REAL(wp) :: zcoef ! local scalar
REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('ssh_nxt')
......@@ -91,12 +98,11 @@ CONTAINS
ENDIF
!
zcoef = 0.5_wp * r1_rho0
! !------------------------------!
! ! After Sea Surface Height !
! !------------------------------!
IF(ln_wd_il) THEN
CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
ENDIF
CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
......@@ -110,7 +116,11 @@ CONTAINS
! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
!
DO_2D_OVR( 1, nn_hls, 1, nn_hls ) ! Loop bounds limited by hdiv definition in div_hor
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#if defined key_RK3
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( r1_rho0 * emp(ji,jj) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#else
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#endif
END_2D
! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
......@@ -123,23 +133,21 @@ CONTAINS
IF ( .NOT.ln_dynspg_ts ) THEN
IF( ln_bdy ) THEN
IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
ENDIF
ENDIF
! !------------------------------!
! ! outputs !
! !------------------------------!
!
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
!
IF( ln_timing ) CALL timing_stop('ssh_nxt')
!
END SUBROUTINE ssh_nxt
SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
!!----------------------------------------------------------------------
!! *** ROUTINE wzv ***
!! *** ROUTINE wzv_MLF ***
!!
!! ** Purpose : compute the now vertical velocity
!!
......@@ -160,12 +168,12 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('wzv')
IF( ln_timing ) CALL timing_start('wzv_MLF')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~~~'
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
ENDIF
......@@ -198,7 +206,7 @@ CONTAINS
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
DEALLOCATE( zhdiv )
DEALLOCATE( zhdiv )
! !=================================!
ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
! !=================================!
......@@ -209,9 +217,146 @@ CONTAINS
ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
! !==========================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
#if defined key_qco
!!gm slightly faster
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
& + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
#else
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
#endif
END_3D
ENDIF
IF( ln_bdy ) THEN
DO jk = 1, jpkm1
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
END_2D
END DO
ENDIF
!
#if defined key_agrif
IF( .NOT. AGRIF_Root() ) THEN
!
! Mask vertical velocity at first/last columns/row
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
!
END DO
!
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('wzv_MLF')
!
END SUBROUTINE wzv_MLF
SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
!!----------------------------------------------------------------------
!! *** ROUTINE wzv_RK3 ***
!!
!! ** Purpose : compute the now vertical velocity
!!
!! ** Method : - Using the incompressibility hypothesis, the vertical
!! velocity is computed by integrating the horizontal divergence
!! from the bottom to the surface minus the scale factor evolution.
!! The boundary conditions are w=0 at the bottom (no flux) and.
!!
!! ** action : pww : now vertical velocity
!!
!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity at Kmm
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
REAL(wp) , DIMENSION(jpi,jpj,jpk) :: ze3div
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('wzv_RK3')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
ENDIF
!
CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
! !------------------------------!
! ! Now Vertical Velocity !
! !------------------------------!
!
! !===============================!
IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
! !===============================!
ALLOCATE( zhdiv(jpi,jpj,jpk) )
!
DO jk = 1, jpkm1
! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
END_2D
END DO
IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
! ! Is it problematic to have a wrong vertical velocity in boundary cells?
! ! Same question holds for hdiv. Perhaps just for security
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( ze3div(ji,jj,jk) + zhdiv(ji,jj,jk) &
& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
!
DEALLOCATE( zhdiv )
! !=================================!
ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
! !=================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ze3div(ji,jj,jk)
END_3D
! !==========================================!
ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
! !==========================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( ze3div(ji,jj,jk) &
& + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
ENDIF
......@@ -263,9 +408,9 @@ CONTAINS
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('wzv')
IF( ln_timing ) CALL timing_stop('wzv_RK3')
!
END SUBROUTINE wzv
END SUBROUTINE wzv_RK3
SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
......
src/OCE/IOM/iom.F90
View file @ 1398e060
......@@ -55,6 +55,7 @@ MODULE iom
#else
LOGICAL, PUBLIC, PARAMETER :: lk_iomput = .FALSE. !: iom_put flag
#endif
LOGICAL, PUBLIC :: l_iom = .TRUE. !: RK3 iom flag prevent writing at stage 1&2
PUBLIC iom_init, iom_init_closedef, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_get_var
PUBLIC iom_chkatt, iom_getatt, iom_putatt, iom_getszuld, iom_rstput, iom_delay_rst, iom_put
PUBLIC iom_use, iom_context_finalize, iom_update_file_name, iom_miss_val
......@@ -98,7 +99,7 @@ MODULE iom
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: iom.F90 15033 2021-06-21 10:24:45Z smasson $
!! $Id: iom.F90 15512 2021-11-15 17:22:03Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......
src/OCE/IOM/iom_nf90.F90
View file @ 1398e060
......@@ -556,6 +556,7 @@ CONTAINS
LOGICAL :: lchunk ! logical switch to activate chunking and compression
! ! when appropriate (currently chunking is applied to 4d fields only)
INTEGER :: idlv ! local variable
CHARACTER(LEN=256) :: ccname ! local variable
!---------------------------------------------------------------------
!
clinfo = ' iom_nf90_rp0123d, file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(cdvar)
......@@ -571,7 +572,10 @@ CONTAINS
! define the dimension variables if it is not already done
DO jd = 1, 2
CALL iom_nf90_check(NF90_INQUIRE_DIMENSION(if90id,jd,iom_file(kiomid)%cn_var(jd),iom_file(kiomid)%dimsz(jd,jd)),clinfo)
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(iom_file(kiomid)%cn_var(jd)), NF90_FLOAT , (/ 1, 2 /), &
ccname = TRIM(iom_file(kiomid)%cn_var(jd))
IF ( ccname == 'x') ccname = 'nav_lon'
IF ( ccname == 'y') ccname = 'nav_lat'
CALL iom_nf90_check(NF90_DEF_VAR( if90id, ccname, NF90_FLOAT , (/ 1, 2 /), &
& iom_file(kiomid)%nvid(jd) ), clinfo)
END DO
iom_file(kiomid)%dimsz(2,1) = iom_file(kiomid)%dimsz(2,2) ! second dim of first variable
......
src/OCE/IOM/restart.F90
View file @ 1398e060
......@@ -52,7 +52,7 @@ MODULE restart
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: restart.F90 15141 2021-07-23 14:20:12Z smasson $
!! $Id: restart.F90 15027 2021-06-19 08:14:22Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -143,7 +143,7 @@ CONTAINS
END SUBROUTINE rst_opn
SUBROUTINE rst_write( kt, Kbb, Kmm )
SUBROUTINE rst_write( kt, Kbb, Kmm, Kaa )
!!---------------------------------------------------------------------
!! *** ROUTINE rstwrite ***
!!
......@@ -157,6 +157,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER, OPTIONAL, INTENT(in) :: Kaa ! ocean time level index required for RK3
!!----------------------------------------------------------------------
!
CALL iom_rstput( kt, nitrst, numrow, 'rdt' , rn_Dt ) ! dynamics time step
......@@ -164,19 +165,27 @@ CONTAINS
IF( .NOT.lwxios ) CALL iom_delay_rst( 'WRITE', 'OCE', numrow ) ! save only ocean delayed global communication variables
!
IF( .NOT.ln_diurnal_only ) THEN
!
#if defined key_RK3
CALL iom_rstput( kt, nitrst, numrow, 'sshn', ssh(:,: ,Kbb) ) ! before fields
CALL iom_rstput( kt, nitrst, numrow, 'un' , uu(:,:,: ,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'vn' , vv(:,:,: ,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'tn' , ts(:,:,:,jp_tem,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'sn' , ts(:,:,:,jp_sal,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'uu_n' , uu_b(:,: ,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'vv_n' , vv_b(:,: ,Kbb) )
#else
CALL iom_rstput( kt, nitrst, numrow, 'sshb', ssh(:,: ,Kbb) ) ! before fields
CALL iom_rstput( kt, nitrst, numrow, 'ub' , uu(:,:,: ,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'vb' , vv(:,:,: ,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'tb' , ts(:,:,:,jp_tem,Kbb) )
CALL iom_rstput( kt, nitrst, numrow, 'sb' , ts(:,:,:,jp_sal,Kbb) )
!
#if ! defined key_RK3
CALL iom_rstput( kt, nitrst, numrow, 'sshn', ssh(:,: ,Kmm) ) ! now fields
CALL iom_rstput( kt, nitrst, numrow, 'sshn', ssh(:,: ,Kmm) ) ! now fields
CALL iom_rstput( kt, nitrst, numrow, 'un' , uu(:,:,: ,Kmm) )
CALL iom_rstput( kt, nitrst, numrow, 'vn' , vv(:,:,: ,Kmm) )
CALL iom_rstput( kt, nitrst, numrow, 'tn' , ts(:,:,:,jp_tem,Kmm) )
CALL iom_rstput( kt, nitrst, numrow, 'sn' , ts(:,:,:,jp_sal,Kmm) )
IF( .NOT.lk_SWE ) CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
#endif
ENDIF
......@@ -264,6 +273,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER :: jk
INTEGER :: id1
REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept ! 3D workspace for QCO
!!----------------------------------------------------------------------
......@@ -284,10 +294,30 @@ CONTAINS
#if defined key_RK3
! !* Read Kbb fields (NB: in RK3 Kmm = Kbb = Nbb)
IF(lwp) WRITE(numout,*) ' Kbb u, v and T-S fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'ub', uu(:,:,: ,Kbb), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vb', vv(:,:,: ,Kbb), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
CALL iom_get( numror, jpdom_auto, 'un' , uu(:,:,: ,Kbb), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vn' , vv(:,:,: ,Kbb), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'tn' , ts(:,:,:,jp_tem,Kbb) )
CALL iom_get( numror, jpdom_auto, 'sn' , ts(:,:,:,jp_sal,Kbb) )
id1 = iom_varid( numror, 'uu_n', ldstop = .FALSE. ) !* check presence
IF( id1 > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'uu_n' , uu_b(:,:,Kbb), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vv_n' , vv_b(:,:,Kbb), cd_type = 'V', psgn = -1._wp )
ELSE
uu_b(:,:,Kbb) = uu(:,:,1,Kbb)*e3u_0(:,:,1)*umask(:,:,1)
vv_b(:,:,Kbb) = vv(:,:,1,Kbb)*e3v_0(:,:,1)*vmask(:,:,1)
DO jk = 2, jpkm1
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) + uu(:,:,jk,Kbb)*e3u_0(:,:,jk)*umask(:,:,jk)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) + vv(:,:,jk,Kbb)*e3v_0(:,:,jk)*vmask(:,:,jk)
END DO
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_0(:,:)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_0(:,:)
ENDIF
uu(:,:,: ,Kmm) = uu(:,:,: ,Kbb) ! Kmm values set to Kbb for initialisation (sbc_ssm_init)
vv(:,:,: ,Kmm) = vv(:,:,: ,Kbb)
ts(:,:,:,:,Kmm) = ts(:,:,:,:,Kbb)
!
uu_b(:,:,Kmm) = uu_b(:,:,Kbb) ! Kmm value set to Kbb for initialisation in Agrif_Regrid
vv_b(:,:,Kmm) = vv_b(:,:,Kbb)
#else
! !* Read Kmm fields (MLF only)
IF(lwp) WRITE(numout,*) ' Kmm u, v and T-S fields read in the restart file'
......@@ -311,23 +341,6 @@ CONTAINS
ENDIF
#endif
!
IF( .NOT.lk_SWE ) THEN
IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'rhop' , rhop ) ! now potential density
ELSE
#if defined key_qco
ALLOCATE( zgdept(jpi,jpj,jpk) )
DO jk = 1, jpk
zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
END DO
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, zgdept )
DEALLOCATE( zgdept )
#else
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept(:,:,:,Kmm) )
#endif
ENDIF
ENDIF
!
END SUBROUTINE rst_read
......@@ -367,10 +380,11 @@ CONTAINS
! !* RK3: Read ssh at Kbb
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' Kbb sea surface height read in the restart file'
CALL iom_get( numror, jpdom_auto, 'sshb' , ssh(:,:,Kbb) )
CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kbb) )
!
! !* RK3: Set ssh at Kmm for AGRIF
ssh(:,:,Kmm) = 0._wp
ssh(:,:,Kmm) = ssh(:,:,Kbb)
!
#else
! !* MLF: Read ssh at Kmm
IF(lwp) WRITE(numout,*)
......@@ -420,14 +434,11 @@ CONTAINS
!
ENDIF
#if defined key_agrif
! Set ghosts points from parent
IF (.NOT.Agrif_Root()) CALL Agrif_istate_ssh( Kbb, Kmm, Kaa, .true. )
#endif
#if defined key_RK3
ssh(:,:,Kmm) = 0._wp !* RK3: set Kmm to 0 for AGRIF
#else
ssh(:,:,Kmm) = ssh(:,:,Kbb) !* MLF: set now values from to before ones
! Set ghosts points from parent
IF (.NOT.Agrif_Root()) CALL Agrif_istate_ssh( Kbb, Kmm, Kaa, .true. )
#endif
!
ssh(:,:,Kmm) = ssh(:,:,Kbb) !* set now values from to before ones
ENDIF
!
! !==========================!
......
src/OCE/ISF/isfhdiv.F90
View file @ 1398e060
......@@ -71,6 +71,7 @@ CONTAINS
!
END SUBROUTINE isf_hdiv
SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
!!----------------------------------------------------------------------
!! *** SUBROUTINE sbc_isf_div ***
......@@ -97,7 +98,11 @@ CONTAINS
!
! compute integrated divergence correction
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
#if defined key_RK3
zhdiv(ji,jj) = pfwf(ji,jj) * r1_rho0 / phtbl(ji,jj)
#else
zhdiv(ji,jj) = 0.5_wp * ( pfwf(ji,jj) + pfwf_b(ji,jj) ) * r1_rho0 / phtbl(ji,jj)
#endif
END_2D
!
! update divergence at each level affected by ice shelf top boundary layer
......
src/OCE/LDF/ldftra.F90
View file @ 1398e060
......@@ -67,6 +67,8 @@ MODULE ldftra
! != Use/diagnose eiv =!
LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag
LOGICAL , PUBLIC :: ln_ldfeiv_dia !: diagnose & output eiv streamfunction and velocity (IOM)
LOGICAL , PUBLIC :: l_ldfeiv_dia !: RK3: modified w.r.t. kstg diagnose & output eiv streamfunction and velocity flag
! != Coefficients =!
INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff.
REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s]
......@@ -97,7 +99,7 @@ MODULE ldftra
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: ldftra.F90 15475 2021-11-05 14:14:45Z cdllod $
!! $Id: ldftra.F90 15512 2021-11-15 17:22:03Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -768,7 +770,11 @@ CONTAINS
END_3D
!
! ! diagnose the eddy induced velocity and associated heat transport
#if defined key_RK3
IF( l_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
#else
IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
#endif
!
END SUBROUTINE ldf_eiv_trp
......@@ -867,7 +873,7 @@ CONTAINS
CALL iom_put( "veiv_heattr" , zztmp * zw2d ) ! heat transport in j-direction
CALL iom_put( "veiv_heattr3d", zztmp * zw3d ) ! heat transport in j-direction
!
IF( iom_use( 'sophteiv' ) ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
IF( iom_use( 'sophteiv' ) .AND. l_diaptr ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
!
zztmp = 0.5_wp * 0.5
IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
......@@ -891,7 +897,7 @@ CONTAINS
CALL iom_put( "veiv_salttr" , zztmp * zw2d ) ! salt transport in j-direction
CALL iom_put( "veiv_salttr3d", zztmp * zw3d ) ! salt transport in j-direction
!
IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
IF( iom_use( 'sopsteiv' ) .AND. l_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
!
!
END SUBROUTINE ldf_eiv_dia
......
src/OCE/SBC/sbcmod.F90
View file @ 1398e060
......@@ -537,7 +537,11 @@ CONTAINS
!
IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
! ! ---------------------------------------- !
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* Restart: read in restart file
#if defined key_RK3
IF( ln_rstart .AND. lk_SWE ) THEN !* RK3 + SWE: Restart: read in restart file
#else
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* MLF: Restart: read in restart file
#endif
IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b ) ! i-stress
CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b ) ! j-stress
......@@ -559,9 +563,17 @@ CONTAINS
sfx_b (:,:) = sfx (:,:)
ENDIF
ENDIF
!
#if defined key_RK3
! ! ---------------------------------------- !
IF( lrst_oce .AND. lk_SWE ) THEN ! RK3: Write in the ocean restart file !
! ! ---------------------------------------- !
#else
! ! ---------------------------------------- !
IF( lrst_oce ) THEN ! Write in the ocean restart file !
IF( lrst_oce ) THEN ! MLF: Write in the ocean restart file !
! ! ---------------------------------------- !
#endif
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ', &
& 'at it= ', kt,' date= ', ndastp
......
src/OCE/SBC/sbcrnf.F90
View file @ 1398e060
......@@ -74,7 +74,7 @@ MODULE sbcrnf
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcrnf.F90 15190 2021-08-13 12:52:50Z gsamson $
!! $Id: sbcrnf.F90 14993 2021-06-14 22:35:18Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -206,13 +206,17 @@ CONTAINS
REAL(wp) :: zfact ! local scalar
!!----------------------------------------------------------------------
!
zfact = 0.5_wp
zfact = 0.5_wp * r1_rho0
!
IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN !== runoff distributed over several levels ==!
IF( ln_linssh ) THEN !* constant volume case : just apply the runoff input flow
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
DO jk = 1, nk_rnf(ji,jj)
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
#if defined key_RK3
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
#else
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj) ! MLF: rnf forcing at Kmm (n)
#endif
END DO
END_2D
ELSE !* variable volume case
......@@ -224,7 +228,11 @@ CONTAINS
END_2D
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! apply the runoff input flow
DO jk = 1, nk_rnf(ji,jj)
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
#if defined key_RK3
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
#else
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj) ! MLF: rnf forcing at Kmm (n)
#endif
END DO
END_2D
ENDIF
......@@ -233,7 +241,11 @@ CONTAINS
h_rnf (ji,jj) = e3t(ji,jj,1,Kmm) ! update h_rnf to be depth of top box
END_2D
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / e3t(ji,jj,1,Kmm)
#if defined key_RK3
phdivn(ji,jj,1) = phdivn(ji,jj,1) - rnf(ji,jj) * r1_rho0 / e3t(ji,jj,1,Kmm) ! RK3: rnf forcing at n+1/2
#else
phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / e3t(ji,jj,1,Kmm) ! MLF: rnf forcing at Kmm (n)
#endif
END_2D
ENDIF
!
......
src/OCE/SBC/sbcssm.F90
View file @ 1398e060
......@@ -64,7 +64,12 @@ CONTAINS
zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
!
! ! ---------------------------------------- !
! !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
! ! otherwise restartability and reproducibility are broken
! ! computed in traqsr only on the inner domain
CALL lbc_lnk( 'sbc_ssm', fraqsr_1lev(:,:), 'T', 1._wp )
!
! ! ---------------------------------------- !
IF( nn_fsbc == 1 ) THEN ! Instantaneous surface fields !
! ! ---------------------------------------- !
ssu_m(:,:) = uu(:,:,1,Kbb)
......@@ -74,7 +79,11 @@ CONTAINS
ENDIF
sss_m(:,:) = zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
#if defined key_RK3
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - ssh_ib(:,:) ! RK3: forcing at n+1/2
#else
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ! MLF: forcing at n (Kmm)
#endif
ELSE ; ssh_m(:,:) = ssh(:,:,Kmm)
ENDIF
!
......@@ -97,7 +106,11 @@ CONTAINS
ENDIF
sss_m(:,:) = zcoef * zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
#if defined key_RK3
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - ssh_ib(:,:) )
#else
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
#endif
ELSE ; ssh_m(:,:) = zcoef * ssh(:,:,Kmm)
ENDIF
!
......@@ -125,7 +138,11 @@ CONTAINS
ENDIF
sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
#if defined key_RK3
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - ssh_ib(:,:)
#else
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
#endif
ELSE ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm)
ENDIF
!
......
src/OCE/TRA/eosbn2.F90
View file @ 1398e060
......@@ -53,7 +53,7 @@ MODULE eosbn2
! !! * Interface
INTERFACE eos
MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
MODULE PROCEDURE eos_insitu_New, eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
END INTERFACE
!
INTERFACE eos_rab
......@@ -181,11 +181,118 @@ MODULE eosbn2
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: eosbn2.F90 15136 2021-07-23 10:07:28Z smasson $
!! $Id: eosbn2.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE eos_insitu_New( pts, Knn, prd )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) from
!! potential temperature and salinity using an equation of state
!! selected in the nameos namelist
!!
!! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
!! with prd in situ density anomaly no units
!! t TEOS10: CT or EOS80: PT Celsius
!! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu
!! z depth meters
!! rho in situ density kg/m^3
!! rho0 reference density kg/m^3
!!
!! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
!!
!! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
!!
!! ln_seos : simplified equation of state
!! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
!! linear case function of T only: rn_alpha<>0, other coefficients = 0
!! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
!! Vallis like equation: use default values of coefficients
!!
!! ** Action : compute prd , the in situ density (no units)
!!
!! References : Roquet et al, Ocean Modelling, in preparation (2014)
!! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
!! TEOS-10 Manual, 2010
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in ) :: pts ! T-S
INTEGER , INTENT(in ) :: Knn ! time-level
REAL(wp), DIMENSION(:,:,: ), INTENT( out) :: prd ! in situ density
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-insitu')
!
SELECT CASE( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Knn) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem,Knn) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal,Knn) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked)
!
END_3D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - 10._wp
zs = pts (ji,jj,jk,jp_sal,Knn) - 35._wp
zh = gdept(ji,jj,jk,Knn)
ztm = tmask(ji,jj,jk)
!
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
& - rn_nu * zt * zs
!
prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked)
END_3D
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('eos-insitu')
!
END SUBROUTINE eos_insitu_New
SUBROUTINE eos_insitu( pts, prd, pdep )
!!
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
......
src/OCE/TRA/traadv.F90
View file @ 1398e060
......@@ -9,6 +9,7 @@ MODULE traadv
!! 3.7 ! 2014-05 (G. Madec) Add 2nd/4th order cases for CEN and FCT schemes
!! - ! 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
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -43,7 +44,7 @@ MODULE traadv
IMPLICIT NONE
PRIVATE
PUBLIC tra_adv ! called by step.F90
PUBLIC tra_adv ! called by step.F90, stpmlf.F90 and stprk3_stg.F90
PUBLIC tra_adv_init ! called by nemogcm.F90
! !!* Namelist namtra_adv *
......@@ -58,43 +59,45 @@ MODULE traadv
INTEGER :: nn_ubs_v ! =2/4 : vertical choice of the order of UBS scheme
LOGICAL :: ln_traadv_qck ! QUICKEST scheme flag
INTEGER, PUBLIC :: nadv ! choice of the type of advection scheme
INTEGER :: nadv ! choice of the type of advection scheme
! ! associated indices:
INTEGER, PARAMETER, PUBLIC :: np_NO_adv = 0 ! no T-S advection
INTEGER, PARAMETER, PUBLIC :: np_CEN = 1 ! 2nd/4th order centered scheme
INTEGER, PARAMETER, PUBLIC :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme
INTEGER, PARAMETER, PUBLIC :: np_MUS = 3 ! MUSCL scheme
INTEGER, PARAMETER, PUBLIC :: np_UBS = 4 ! 3rd order Upstream Biased Scheme
INTEGER, PARAMETER, PUBLIC :: np_QCK = 5 ! QUICK scheme
INTEGER, PARAMETER :: np_NO_adv = 0 ! no T-S advection
INTEGER, PARAMETER :: np_CEN = 1 ! 2nd/4th order centered scheme
INTEGER, PARAMETER :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme
INTEGER, PARAMETER :: np_MUS = 3 ! MUSCL scheme
INTEGER, PARAMETER :: np_UBS = 4 ! 3rd order Upstream Biased Scheme
INTEGER, PARAMETER :: np_QCK = 5 ! QUICK scheme
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: traadv.F90 15073 2021-07-02 14:20:14Z clem $
!! $Id: traadv.F90 15514 2021-11-16 08:58:22Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs )
SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv ***
!!
!! ** Purpose : compute the ocean tracer advection trend.
!!
!! ** Method : - Update (uu(:,:,:,Krhs),vv(:,:,:,Krhs)) with the advection term following nadv
!! ** 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(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, 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
! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww ! 3D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww ! 3D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
LOGICAL :: lskip
LOGICAL :: lskip
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_adv')
......@@ -114,23 +117,35 @@ CONTAINS
lskip = .TRUE.
ENDIF
ENDIF
!
IF( .NOT. lskip ) THEN
! !== effective transport ==!
! !== effective advective transport ==!
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (uu,vv,ww)
zptu => pau(:,:,:)
zptv => pav(:,:,:)
zptw => paw(:,:,:)
ELSE ! MLF: advective velocity = (uu,vv,ww)
zptu => uu(:,:,:,Kmm)
zptv => vv(:,:,:,Kmm)
zptw => ww(:,:,: )
ENDIF
!
IF( ln_wave .AND. ln_sdw ) THEN
DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) )
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) )
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) )
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( zptv(ji,jj,jk) + vsd(ji,jj,jk) )
END_3D
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
zww(ji,jj,jk) = e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
zww(ji,jj,jk) = e1e2t(ji,jj) * ( zptw(ji,jj,jk) + wsd(ji,jj,jk) )
END_3D
ELSE
DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) ! eulerian transport only
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * zptu(ji,jj,jk) ! eulerian transport only
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * zptv(ji,jj,jk)
END_3D
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
zww(ji,jj,jk) = e1e2t(ji,jj) * ww(ji,jj,jk)
zww(ji,jj,jk) = e1e2t(ji,jj) * zptw(ji,jj,jk)
END_3D
ENDIF
!
......@@ -142,7 +157,7 @@ CONTAINS
ENDIF
!
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
zuu(ji,jj,jpk) = 0._wp ! no transport trough the bottom
zuu(ji,jj,jpk) = 0._wp ! no transport trough the bottom
zvv(ji,jj,jpk) = 0._wp
zww(ji,jj,jpk) = 0._wp
END_2D
......@@ -153,15 +168,17 @@ CONTAINS
IF( ln_mle ) CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm ) ! add the mle transport (if necessary)
!
! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
CALL iom_put( "uocetr_eff", zuu ) ! output effective transport
CALL iom_put( "vocetr_eff", zvv )
CALL iom_put( "wocetr_eff", zww )
IF( l_iom ) THEN
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
CALL iom_put( "uocetr_eff", zuu ) ! output effective transport
CALL iom_put( "vocetr_eff", zvv )
CALL iom_put( "wocetr_eff", zww )
ENDIF
ENDIF
!
!!gm ???
! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF
IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF
!!gm ???
!
......@@ -245,8 +262,8 @@ CONTAINS
WRITE(numout,*) ' Namelist namtra_adv : chose a advection scheme for tracers'
WRITE(numout,*) ' No advection on T & S ln_traadv_OFF = ', ln_traadv_OFF
WRITE(numout,*) ' centered scheme ln_traadv_cen = ', ln_traadv_cen
WRITE(numout,*) ' horizontal 2nd/4th order nn_cen_h = ', nn_fct_h
WRITE(numout,*) ' vertical 2nd/4th order nn_cen_v = ', nn_fct_v
WRITE(numout,*) ' horizontal 2nd/4th order nn_cen_h = ', nn_cen_h
WRITE(numout,*) ' vertical 2nd/4th order nn_cen_v = ', nn_cen_v
WRITE(numout,*) ' Flux Corrected Transport scheme ln_traadv_fct = ', ln_traadv_fct
WRITE(numout,*) ' horizontal 2nd/4th order nn_fct_h = ', nn_fct_h
WRITE(numout,*) ' vertical 2nd/4th order nn_fct_v = ', nn_fct_v
......