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
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1390 additions and 614 deletions
src/OCE/DIU/diu_bulk.F90
View file @ 005ca9c1
...@@ -64,7 +64,7 @@ CONTAINS ...@@ -64,7 +64,7 @@ CONTAINS
IF( ln_diurnal ) THEN IF( ln_diurnal ) THEN
! !
ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) ) ALLOCATE( x_dsst(A2D(0)), x_solfrac(A2D(0)) )
! !
x_solfrac = 0._wp ! Initialise the solar fraction x_solfrac = 0._wp ! Initialise the solar fraction
x_dsst = 0._wp x_dsst = 0._wp
...@@ -92,25 +92,25 @@ CONTAINS ...@@ -92,25 +92,25 @@ CONTAINS
!! ** Reference : Refinements to a prognostic scheme of skin sea surface !! ** Reference : Refinements to a prognostic scheme of skin sea surface
!! temperature, Takaya et al, JGR, 2010 !! temperature, Takaya et al, JGR, 2010
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: psolflux ! solar flux (Watts) REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts) REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2) REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3) REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pLa ! Langmuir number REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m) REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m) REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pmu ! mu parameter REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
! !
INTEGER :: ji,jj INTEGER :: ji,jj
LOGICAL :: ll_calcfrac LOGICAL :: ll_calcfrac
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel ! friction velocity REAL(wp), DIMENSION(A2D(0)) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla REAL(wp), DIMENSION(A2D(0)) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(jpi,jpj) :: z_abflux ! absorbed flux REAL(wp), DIMENSION(A2D(0)) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(jpi,jpj) :: z_fla ! Langmuir function value REAL(wp), DIMENSION(A2D(0)) :: z_fla ! Langmuir function value
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! Set optional arguments to their defaults ! Set optional arguments to their defaults
...@@ -129,14 +129,14 @@ CONTAINS ...@@ -129,14 +129,14 @@ CONTAINS
! If not done already, calculate the solar fraction ! If not done already, calculate the solar fraction
IF ( kt==nit000 ) THEN IF ( kt==nit000 ) THEN
DO_2D( 1, 1, 1, 1 ) DO_2D( 0, 0, 0, 0 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) & IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( smask0(ji,jj) == 1._wp ) ) &
& x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) ) & x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
END_2D END_2D
ENDIF ENDIF
! convert solar flux and heat flux to absorbed flux ! convert solar flux and heat flux to absorbed flux
WHERE ( tmask(:,:,1) == 1._wp) WHERE ( smask0(:,:) == 1._wp)
z_abflux(:,:) = ( x_solfrac(:,:) * psolflux (:,:)) + pqflux(:,:) z_abflux(:,:) = ( x_solfrac(:,:) * psolflux (:,:)) + pqflux(:,:)
ELSEWHERE ELSEWHERE
z_abflux(:,:) = 0._wp z_abflux(:,:) = 0._wp
...@@ -147,8 +147,8 @@ CONTAINS ...@@ -147,8 +147,8 @@ CONTAINS
ENDWHERE ENDWHERE
! Calculate the friction velocity ! Calculate the friction velocity
WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) ) WHERE ( (ptauflux(:,:) /= 0) .AND. ( smask0(:,:) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) ) z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(A2D(0)) )
ELSEWHERE ELSEWHERE
z_fvel(:,:) = 0._wp z_fvel(:,:) = 0._wp
ENDWHERE ENDWHERE
...@@ -157,7 +157,7 @@ CONTAINS ...@@ -157,7 +157,7 @@ CONTAINS
! Calculate the Langmuir function value ! Calculate the Langmuir function value
WHERE ( tmask(:,:,1) == 1.) WHERE ( smask0(:,:) == 1.)
z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) ) z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )
ELSEWHERE ELSEWHERE
z_fla(:,:) = 0._wp z_fla(:,:) = 0._wp
...@@ -176,16 +176,16 @@ CONTAINS ...@@ -176,16 +176,16 @@ CONTAINS
IMPLICIT NONE IMPLICIT NONE
! Function definition ! Function definition
REAL(wp), DIMENSION(jpi,jpj) :: t_imp REAL(wp), DIMENSION(A2D(0)) :: t_imp
! Dummy variables ! Dummy variables
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst ! Delta SST REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux ! Heat forcing REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel ! Friction velocity REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla ! Langmuir number REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu ! Structure parameter REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick ! Layer thickness REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pthick ! Layer thickness
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho ! Water density REAL(wp), DIMENSION(jpi,jpj),INTENT(IN) :: prho ! Water density
! Local variables ! Local variables
REAL(wp) :: z_olength ! Obukhov length REAL(wp) :: z_olength ! Obukhov length
...@@ -198,10 +198,10 @@ CONTAINS ...@@ -198,10 +198,10 @@ CONTAINS
INTEGER :: ji,jj INTEGER :: ji,jj
DO_2D( 1, 1, 1, 1 ) DO_2D( 0, 0, 0, 0 )
! Only calculate outside tmask ! Only calculate outside tmask
IF ( tmask(ji,jj,1) /= 1._wp ) THEN IF ( smask0(ji,jj) /= 1._wp ) THEN
t_imp(ji,jj) = 0._wp t_imp(ji,jj) = 0._wp
CYCLE CYCLE
END IF END IF
......
src/OCE/DIU/diu_coolskin.F90
View file @ 005ca9c1
...@@ -59,7 +59,7 @@ MODULE diu_coolskin ...@@ -59,7 +59,7 @@ MODULE diu_coolskin
!! ** Reference : !! ** Reference :
!! !!
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) ) ALLOCATE( x_csdsst(A2D(0)), x_csthick(A2D(0)) )
x_csdsst = 0. x_csdsst = 0.
x_csthick = 0. x_csthick = 0.
! !
...@@ -77,16 +77,16 @@ MODULE diu_coolskin ...@@ -77,16 +77,16 @@ MODULE diu_coolskin
!! ** Reference : !! ** Reference :
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! Dummy variables ! Dummy variables
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux ! Heat (non-solar)(Watts) REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: psqflux ! Heat (non-solar)(Watts)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux ! Wind stress (kg/ m s^2) REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: pstauflux ! Wind stress (kg/ m s^2)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho ! Water density (kg/m^3) REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho ! Water density (kg/m^3)
REAL(wp), INTENT(IN) :: pDt ! Time-step REAL(wp), INTENT(IN) :: pDt ! Time-step
! Local variables ! Local variables
REAL(wp), DIMENSION(jpi,jpj) :: z_fv ! Friction velocity REAL(wp), DIMENSION(A2D(0)) :: z_fv ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: z_gamma ! Dimensionless function of wind speed REAL(wp), DIMENSION(A2D(0)) :: z_gamma ! Dimensionless function of wind speed
REAL(wp), DIMENSION(jpi,jpj) :: z_lamda ! Sauders (dimensionless) proportionality constant REAL(wp), DIMENSION(A2D(0)) :: z_lamda ! Sauders (dimensionless) proportionality constant
REAL(wp), DIMENSION(jpi,jpj) :: z_wspd ! Wind speed (m/s) REAL(wp), DIMENSION(A2D(0)) :: z_wspd ! Wind speed (m/s)
REAL(wp) :: z_ztx ! Temporary u wind stress REAL(wp) :: z_ztx ! Temporary u wind stress
REAL(wp) :: z_zty ! Temporary v wind stress REAL(wp) :: z_zty ! Temporary v wind stress
REAL(wp) :: z_zmod ! Temporary total wind stress REAL(wp) :: z_zmod ! Temporary total wind stress
...@@ -96,10 +96,10 @@ MODULE diu_coolskin ...@@ -96,10 +96,10 @@ MODULE diu_coolskin
! !
IF( .NOT. (ln_blk .OR. ln_abl) ) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing") IF( .NOT. (ln_blk .OR. ln_abl) ) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
! !
DO_2D( 1, 1, 1, 1 ) DO_2D( 0, 0, 0, 0 )
! !
! Calcualte wind speed from wind stress and friction velocity ! Calcualte wind speed from wind stress and friction velocity
IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN IF( smask0(ji,jj) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) ) z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) ) z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
ELSE ELSE
...@@ -108,28 +108,28 @@ MODULE diu_coolskin ...@@ -108,28 +108,28 @@ MODULE diu_coolskin
ENDIF ENDIF
! !
! Calculate gamma function which is dependent upon wind speed ! Calculate gamma function which is dependent upon wind speed
IF( tmask(ji,jj,1) == 1. ) THEN IF( smask0(ji,jj) == 1. ) THEN
IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5 IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10. IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6. IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
ENDIF ENDIF
! !
! Calculate lamda function ! Calculate lamda function
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v ) z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
ELSE ELSE
z_lamda(ji,jj) = 0. z_lamda(ji,jj) = 0.
ENDIF ENDIF
! !
! Calculate the cool skin thickness - only when heat flux is out of the ocean ! Calculate the cool skin thickness - only when heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj) x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
ELSE ELSE
x_csthick(ji,jj) = 0. x_csthick(ji,jj) = 0.
ENDIF ENDIF
! !
! Calculate the cool skin correction - only when the heat flux is out of the ocean ! Calculate the cool skin correction - only when the heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN IF( smask0(ji,jj) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
ELSE ELSE
x_csdsst(ji,jj) = 0. x_csdsst(ji,jj) = 0.
......
src/OCE/DIU/step_diu.F90
View file @ 005ca9c1
...@@ -46,8 +46,7 @@ MODULE step_diu ...@@ -46,8 +46,7 @@ MODULE step_diu
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER :: jk ! dummy loop indices INTEGER :: jk ! dummy loop indices
INTEGER :: indic ! error indicator if < 0 INTEGER :: indic ! error indicator if < 0
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel_bkginc, z_hflux_bkginc INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
!! --------------------------------------------------------------------- !! ---------------------------------------------------------------------
IF(ln_diurnal_only) THEN IF(ln_diurnal_only) THEN
......
src/OCE/DOM/dom_oce.F90
View file @ 005ca9c1
...@@ -27,6 +27,8 @@ MODULE dom_oce ...@@ -27,6 +27,8 @@ MODULE dom_oce
PUBLIC dom_oce_alloc ! Called from nemogcm.F90 PUBLIC dom_oce_alloc ! Called from nemogcm.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! time & space domain namelist !! time & space domain namelist
!! ---------------------------- !! ----------------------------
...@@ -74,10 +76,10 @@ MODULE dom_oce ...@@ -74,10 +76,10 @@ MODULE dom_oce
INTEGER :: nn_ltile_i, nn_ltile_j INTEGER :: nn_ltile_i, nn_ltile_j
! Domain tiling ! Domain tiling
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsj_a ! INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsj_a !
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntej_a ! INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntej_a !
LOGICAL, PUBLIC :: l_istiled ! whether tiling is currently active or not LOGICAL, PUBLIC :: l_istiled ! whether tiling is currently active or not
! !: domain MPP decomposition parameters ! !: domain MPP decomposition parameters
...@@ -85,32 +87,30 @@ MODULE dom_oce ...@@ -85,32 +87,30 @@ MODULE dom_oce
INTEGER , PUBLIC :: narea !: number for local area (starting at 1) = MPI rank + 1 INTEGER , PUBLIC :: narea !: number for local area (starting at 1) = MPI rank + 1
INTEGER, PUBLIC :: nidom !: IOIPSL things... INTEGER, PUBLIC :: nidom !: IOIPSL things...
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig !: local ==> global domain, including halos (jpiglo), i-index INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mig !: local ==> global domain, i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg !: local ==> global domain, including halos (jpjglo), j-index INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mjg !: local ==> global domain, j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig0 !: local ==> global domain, excluding halos (Ni0glo), i-index INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mi0, mi1 !: global ==> local domain, i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg0 !: local ==> global domain, excluding halos (Nj0glo), j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mi0, mi1 !: global, including halos (jpiglo) ==> local domain i-index
! !: (mi0=1 and mi1=0 if global index not in local domain) ! !: (mi0=1 and mi1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mj0, mj1 !: global, including halos (jpjglo) ==> local domain j-index INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mj0, mj1 !: global ==> local domain, j-index
! !: (mj0=1 and mj1=0 if global index not in local domain) ! !: (mj0=1 and mj1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nfimpp, nfproc, nfjpi INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: nfimpp, nfproc, nfjpi, nfni_0
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! horizontal curvilinear coordinate and scale factors !! horizontal curvilinear coordinate and scale factors
!! --------------------------------------------------------------------- !! ---------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m] REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m] REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m] REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m] REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
! !
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
! !
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s]
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! vertical coordinate and scale factors !! vertical coordinate and scale factors
...@@ -130,75 +130,77 @@ MODULE dom_oce ...@@ -130,75 +130,77 @@ MODULE dom_oce
LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate
LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
! ! reference scale factors ! ! reference scale factors
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m]
! ! time-dependent scale factors (domvvl) ! ! time-dependent scale factors (domvvl)
#if defined key_qco || defined key_linssh #if defined key_qco || defined key_linssh
#else #else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
#endif #endif
! ! time-dependent ratio ssh / h_0 (domqco) ! ! time-dependent ratio ssh / h_0 (domqco)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
! ! reference depths of cells ! ! reference depths of cells
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
! ! time-dependent depths of cells (domvvl) ! ! time-dependent depths of cells (domvvl)
#if defined key_qco || defined key_linssh #if defined key_qco || defined key_linssh
#else #else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: gdept, gdepw
#endif #endif
! ! reference heights of ocean water column and its inverse ! ! reference heights of ocean water column and its inverse
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
! ! time-dependent heights of ocean water column (domvvl) ! ! time-dependent heights of ocean water column (domvvl)
#if defined key_qco || defined key_linssh #if defined key_qco || defined key_linssh
#else #else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: t-points [m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht !: t-points [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m] REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m]
#endif #endif
INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1)
INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1) INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1)
!! 1D reference vertical coordinate !! 1D reference vertical coordinate
!! =-----------------====------ !! =-----------------====------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m) REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep, bathy REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: risfdep, bathy
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! masks, top and bottom ocean point position !! masks, top and bottom ocean point position
!! --------------------------------------------------------------------- !! ---------------------------------------------------------------------
!!gm Proposition of new name for top/bottom vertical indices !!gm Proposition of new name for top/bottom vertical indices
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF) ! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF)
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level ! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level
!!gm !!gm
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF) INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0 !: surface mask at T-pts on inner domain
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0_i !: equivalent of tmask_i for inner domain
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! calendar variables !! calendar variables
...@@ -326,7 +328,7 @@ CONTAINS ...@@ -326,7 +328,7 @@ CONTAINS
ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii) ) ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii) )
! !
ii = ii+1 ii = ii+1
ALLOCATE( tmask_i(jpi,jpj) , & ALLOCATE( tmask_i(jpi,jpj) , smask0(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(0)) , smask0_i(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(0)) , &
& ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
& mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , mbkf(jpi,jpj) , STAT=ierr(ii) ) & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , mbkf(jpi,jpj) , STAT=ierr(ii) )
! !
......
src/OCE/DOM/dommsk.F90
View file @ 005ca9c1
...@@ -144,7 +144,8 @@ CONTAINS ...@@ -144,7 +144,8 @@ CONTAINS
tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
END_3D END_3D
ENDIF ENDIF
smask0(:,:) = tmask(A2D(0),1)
! Ocean/land mask at u-, v-, and f-points (computed from tmask) ! Ocean/land mask at u-, v-, and f-points (computed from tmask)
! ---------------------------------------- ! ----------------------------------------
! NB: at this point, fmask is designed for free slip lateral boundary condition ! NB: at this point, fmask is designed for free slip lateral boundary condition
...@@ -194,7 +195,8 @@ CONTAINS ...@@ -194,7 +195,8 @@ CONTAINS
! -------------------- ! --------------------
! !
CALL dom_uniq( tmask_i, 'T' ) CALL dom_uniq( tmask_i, 'T' )
tmask_i(:,:) = ssmask(:,:) * tmask_i(:,:) tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)
smask0_i(:,:) = tmask_i(A2D(0))
! Lateral boundary conditions on velocity (modify fmask) ! Lateral boundary conditions on velocity (modify fmask)
! --------------------------------------- ! ---------------------------------------
......
src/OCE/DOM/domtile.F90
View file @ 005ca9c1
...@@ -52,16 +52,6 @@ CONTAINS ...@@ -52,16 +52,6 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
IF( ln_tile .AND. nn_hls /= 2 ) CALL ctl_stop('dom_tile_init: Tiling is only supported for nn_hls = 2') IF( ln_tile .AND. nn_hls /= 2 ) CALL ctl_stop('dom_tile_init: Tiling is only supported for nn_hls = 2')
ntile = 0 ! Initialise to full domain
nijtile = 1
ntsi = Nis0
ntsj = Njs0
ntei = Nie0
ntej = Nje0
nthl = 0
nthr = 0
nthb = 0
ntht = 0
l_istiled = .FALSE. l_istiled = .FALSE.
IF( ln_tile ) THEN ! Calculate tile domain indices IF( ln_tile ) THEN ! Calculate tile domain indices
......
src/OCE/DOM/domvvl.F90
View file @ 005ca9c1
...@@ -282,8 +282,8 @@ CONTAINS ...@@ -282,8 +282,8 @@ CONTAINS
IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1 ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1
ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls
frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp frq_rst_e3t( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 0.0_wp
frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt frq_rst_hdv( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 1.e0_wp / rn_Dt
ENDIF ENDIF
ENDIF ENDIF
ENDIF ENDIF
......
src/OCE/DOM/domzgr.F90
View file @ 005ca9c1
...@@ -130,14 +130,14 @@ CONTAINS ...@@ -130,14 +130,14 @@ CONTAINS
! !
zmsk(:,:) = 1._wp ! default: no closed boundaries zmsk(:,:) = 1._wp ! default: no closed boundaries
IF( .NOT. l_Iperio ) THEN ! E-W closed: IF( .NOT. l_Iperio ) THEN ! E-W closed:
zmsk( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = 0._wp ! first column of inner global domain at 0 zmsk( mi0( 1+nn_hls,nn_hls):mi1( 1+nn_hls,nn_hls),:) = 0._wp ! first column of inner global domain at 0
zmsk( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0._wp ! last column of inner global domain at 0 zmsk( mi0(jpiglo-nn_hls,nn_hls):mi1(jpiglo-nn_hls,nn_hls),:) = 0._wp ! last column of inner global domain at 0
ENDIF ENDIF
IF( .NOT. l_Jperio ) THEN ! S closed: IF( .NOT. l_Jperio ) THEN ! S closed:
zmsk(:,mj0( 1+nn_hls):mj1( 1+nn_hls) ) = 0._wp ! first line of inner global domain at 0 zmsk(:,mj0( 1+nn_hls,nn_hls):mj1( 1+nn_hls,nn_hls) ) = 0._wp ! first line of inner global domain at 0
ENDIF ENDIF
IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed: IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed:
zmsk(:,mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls) ) = 0._wp ! last line of inner global domain at 0 zmsk(:,mj0(jpjglo-nn_hls,nn_hls):mj1(jpjglo-nn_hls,nn_hls) ) = 0._wp ! last line of inner global domain at 0
ENDIF ENDIF
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) ) k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
......
src/OCE/DOM/dtatsd.F90
View file @ 005ca9c1
...@@ -161,8 +161,8 @@ CONTAINS ...@@ -161,8 +161,8 @@ CONTAINS
ij0 = 101 + nn_hls ; ij1 = 109 + nn_hls ! Reduced T & S in the Alboran Sea ij0 = 101 + nn_hls ; ij1 = 109 + nn_hls ! Reduced T & S in the Alboran Sea
ii0 = 141 + nn_hls - 1 ; ii1 = 155 + nn_hls - 1 ii0 = 141 + nn_hls - 1 ; ii1 = 155 + nn_hls - 1
IF( sf_tsd(jp_tem)%ln_tint .OR. irec_n(jp_tem) /= irec_b(jp_tem) ) THEN IF( sf_tsd(jp_tem)%ln_tint .OR. irec_n(jp_tem) /= irec_b(jp_tem) ) THEN
DO jj = mj0(ij0), mj1(ij1) DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0), mi1(ii1) DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
...@@ -172,8 +172,8 @@ CONTAINS ...@@ -172,8 +172,8 @@ CONTAINS
ENDIF ENDIF
! !
IF( sf_tsd(jp_sal)%ln_tint .OR. irec_n(jp_sal) /= irec_b(jp_sal) ) THEN IF( sf_tsd(jp_sal)%ln_tint .OR. irec_n(jp_sal) /= irec_b(jp_sal) ) THEN
DO jj = mj0(ij0), mj1(ij1) DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0), mi1(ii1) DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
...@@ -185,9 +185,9 @@ CONTAINS ...@@ -185,9 +185,9 @@ CONTAINS
! !
ij0 = 87 + nn_hls ; ij1 = 96 + nn_hls ! Reduced temperature in Red Sea ij0 = 87 + nn_hls ; ij1 = 96 + nn_hls ! Reduced temperature in Red Sea
ii0 = 148 + nn_hls - 1 ; ii1 = 160 + nn_hls - 1 ii0 = 148 + nn_hls - 1 ; ii1 = 160 + nn_hls - 1
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 14:20 ) = 6.0_wp
ENDIF ENDIF
ENDIF ENDIF
!!gm end !!gm end
......
src/OCE/DYN/dynadv.F90
View file @ 005ca9c1
...@@ -7,6 +7,7 @@ MODULE dynadv ...@@ -7,6 +7,7 @@ MODULE dynadv
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase !! 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 !! 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.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 ...@@ -50,7 +51,7 @@ MODULE dynadv
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
CONTAINS 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 *** !! *** ROUTINE dyn_adv ***
!! !!
...@@ -64,7 +65,6 @@ CONTAINS ...@@ -64,7 +65,6 @@ CONTAINS
!! (see dynvor module). !! (see dynvor module).
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time step and level indices 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(:,:,:), 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. REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum Eq.
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -73,12 +73,23 @@ CONTAINS ...@@ -73,12 +73,23 @@ CONTAINS
! !
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==! SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 ) != vector form =! CASE( np_VEC_c2 ) != vector form =!
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! horizontal gradient of kinetic energy ! !* horizontal gradient of kinetic energy
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) ! vertical advection 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 =! CASE( np_FLX_c2 ) != flux form =!
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 2nd order centered scheme CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw ) !* 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) 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 END SELECT
! !
IF( ln_timing ) CALL timing_stop( 'dyn_adv' ) IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
View file @ 005ca9c1
...@@ -6,6 +6,7 @@ MODULE dynadv_cen2 ...@@ -6,6 +6,7 @@ MODULE dynadv_cen2
!!====================================================================== !!======================================================================
!! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code !! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option !! 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 ...@@ -35,7 +36,7 @@ MODULE dynadv_cen2
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
CONTAINS 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 *** !! *** ROUTINE dyn_adv_cen2 ***
!! !!
...@@ -51,15 +52,17 @@ CONTAINS ...@@ -51,15 +52,17 @@ CONTAINS
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend !! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kmm, Krhs ! ocean time-step and level indices 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(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 REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
! !
INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zzv ! local scalars REAL(wp) :: zzu, zzfu_kp1 ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu REAL(wp) :: zzv, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw 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(:,:,:) , 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 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
...@@ -71,8 +74,9 @@ CONTAINS ...@@ -71,8 +74,9 @@ CONTAINS
ENDIF ENDIF
! !
IF( l_trddyn ) THEN ! trends: store the input trends IF( l_trddyn ) THEN ! trends: store the input trends
zfu_uw(:,:,:) = puu(:,:,:,Krhs) ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zfv_vw(:,:,:) = pvv(:,:,:,Krhs) zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF ENDIF
! !
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
...@@ -89,84 +93,86 @@ CONTAINS ...@@ -89,84 +93,86 @@ CONTAINS
! !
DO jk = 1, jpkm1 ! horizontal transport DO jk = 1, jpkm1 ! horizontal transport
DO_2D( 1, 1, 1, 1 ) 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) zfu(ji,jj) = 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) zfv(ji,jj) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D END_2D
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point) 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) ) 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 ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,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 ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,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,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,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 END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes 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) & puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) & & + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm) & / 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) & pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) & & + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm) & / e3v(ji,jj,jk,Kmm)
END_2D END_2D
END DO END DO
! !
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:) zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:) zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm ) CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zfu_t(:,:,:) = puu(:,:,:,Krhs) zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF ENDIF
! !
IF( PRESENT( no_zad ) ) THEN !== No vertical advection ==! (except if linear free surface) ! !== Vertical advection ==!
! == !
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0 ! ! surface vertical fluxes
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) IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
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) DO_2D( 0, 0, 0, 0 )
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) & 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)
& / e3u(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)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) & END_2D
& / e3v(ji,jj,1,Kmm) ELSE ! non linear free: surface advective fluxes set to zero
END_2D DO_2D( 0, 0, 0, 0 )
ENDIF zfu_uw(ji,jj) = 0._wp
! zfv_vw(ji,jj) = 0._wp
ELSE !== Vertical advection ==! 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 DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp
END_2D END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0 DO_2D( 0, 0, 0, 0 )
DO_2D( 0, 0, 0, 0 ) ! ! vertical flux at level k+1
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) zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,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) zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
END_2D ! ! divergence of vertical momentum flux
ENDIF puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
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) & / 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) & 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) & / e3v(ji,jj,jk,Kmm)
END_3D ! ! store vertical flux for next level calculation
! zfu_uw(ji,jj) = zzfu_kp1
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic zfv_vw(ji,jj) = zzfv_kp1
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:) END_2D
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:) END DO
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm ) !
ENDIF jk = jpkm1
! ! Control print DO_2D( 0, 0, 0, 0 )
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, & puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) & / 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 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 END SUBROUTINE dyn_adv_cen2
......
src/OCE/DYN/dynadv_ubs.F90
View file @ 005ca9c1
...@@ -6,6 +6,7 @@ MODULE dynadv_ubs ...@@ -6,6 +6,7 @@ MODULE dynadv_ubs
!!====================================================================== !!======================================================================
!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code !! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option !! 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
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -29,7 +30,8 @@ MODULE dynadv_ubs ...@@ -29,7 +30,8 @@ MODULE dynadv_ubs
REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred
PUBLIC dyn_adv_ubs ! routine called by step.F90 PUBLIC dyn_adv_ubs ! routine called by dynadv.F90
PUBLIC dyn_adv_ubs_hls1 ! routine called by dynadv.F90
!! * Substitutions !! * Substitutions
# include "do_loop_substitute.h90" # include "do_loop_substitute.h90"
...@@ -41,7 +43,243 @@ MODULE dynadv_ubs ...@@ -41,7 +43,243 @@ MODULE dynadv_ubs
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
CONTAINS CONTAINS
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : The scheme is the one implemeted in ROMS. It depends
!! on two parameter gamma1 and gamma2. The former control the
!! upstream baised part of the scheme and the later the centred
!! part: gamma1 = 0 pure centered (no diffusive part)
!! = 1/4 Quick scheme
!! = 1/3 3rd order Upstream biased scheme
!! gamma2 = 0 2nd order finite differencing
!! = 1/32 4th order finite differencing
!! For stability reasons, the first term of the fluxes which cor-
!! responds to a second order centered scheme is evaluated using
!! the now velocity (centered in time) while the second term which
!! is the diffusive part of the scheme, is evaluated using the
!! before velocity (forward in time).
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! 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 , Kbb, Kmm, Krhs ! ocean time-step and level indices
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, zui, zfuj, zl_u, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zvj, zfvi, zl_v, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(2)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(2)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(2),2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(2),2) :: zlv_vv, zlv_vu
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
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
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)
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( 2, 2, 2, 2 )
zfu(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
!
DO_2D( 1, 1, 1, 1 ) ! laplacian
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility (north fold)
!! zlu_uu(ji,jj,1) = ( ( puu(ji+1,jj ,jk,Kbb) + puu(ji-1,jj ,jk,Kbb) ) - 2._wp * puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
!! zlv_vv(ji,jj,1) = ( ( pvv(ji ,jj+1,jk,Kbb) + pvv(ji ,jj-1,jk,Kbb) ) - 2._wp * pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
zlu_uu(ji,jj,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) ) * umask(ji ,jj ,jk)
zlv_vv(ji,jj,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,1) = ( puu(ji ,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( puu(ji ,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,1) = ( pvv(ji+1,jj ,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( pvv(ji ,jj ,jk,Kbb) - pvv(ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
!
!! zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) + zfu(ji-1,jj ) ) - 2._wp * zfu(ji,jj) ) * umask(ji ,jj ,jk)
!! zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) + zfv(ji ,jj-1) ) - 2._wp * zfv(ji,jj) ) * vmask(ji ,jj ,jk)
zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) - zfu(ji ,jj ) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( zfu(ji-1,jj ) - zfu(ji ,jj ) &
& ) ) * umask(ji ,jj ,jk)
zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) - zfv(ji ,jj ) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( zfv(ji ,jj-1) - zfv(ji ,jj ) &
& ) ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,2) = ( zfu(ji ,jj+1) - zfu(ji ,jj ) ) * fmask(ji ,jj ,jk) &
& - ( zfu(ji ,jj ) - zfu(ji ,jj-1) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,2) = ( zfv(ji+1,jj ) - zfv(ji ,jj ) ) * fmask(ji ,jj ,jk) &
& - ( zfv(ji ,jj ) - zfv(ji-1,jj ) ) * fmask(ji-1,jj ,jk)
END_2D
!
! ! ====================== !
! ! Horizontal advection !
! ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj) = 0.25_wp * zfu(ji,jj)
zfv(ji,jj) = 0.25_wp * zfv(ji,jj)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
!
IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,1)
ELSE ; zl_u = zlu_uu(ji+1,jj,1)
ENDIF
IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,1)
ELSE ; zl_v = zlv_vv(ji,jj+1,1)
ENDIF
!
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj ) - gamma2 * ( zlu_uu(ji,jj,2) + zlu_uu(ji+1,jj ,2) ) ) &
& * ( zui - gamma1 * zl_u )
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji ,jj+1) - gamma2 * ( zlv_vv(ji,jj,2) + zlv_vv(ji ,jj+1,2) ) ) &
& * ( zvj - gamma1 * zl_v )
!
zfuj = ( zfu(ji,jj) + zfu(ji ,jj+1) )
zfvi = ( zfv(ji,jj) + zfv(ji+1,jj ) )
IF( zfuj > 0 ) THEN ; zl_v = zlv_vu(ji ,jj,1)
ELSE ; zl_v = zlv_vu(ji+1,jj,1)
ENDIF
IF( zfvi > 0 ) THEN ; zl_u = zlu_uv(ji,jj ,1)
ELSE ; zl_u = zlu_uv(ji,jj+1,1)
ENDIF
!
zfv_f(ji ,jj ) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,2) + zlv_vu(ji+1,jj ,2) ) ) &
& * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
zfu_f(ji ,jj ) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,2) + zlu_uv(ji ,jj+1,2) ) ) &
& * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
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) - 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 ) - 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
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
! ! ==================== !
! ! Vertical advection !
! ! ==================== !
!
#define zfu_uw zfu_t
#define zfv_vw zfv_t
#define zfw zfu
!
! ! 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, 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
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) - 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 ! compute last level (zzfu_kp1 = 0)
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
!
#undef zfu_uw
#undef zfv_vw
#undef zfw
! ! 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
SUBROUTINE dyn_adv_ubs_hls1( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs *** !! *** ROUTINE dyn_adv_ubs ***
!! !!
...@@ -73,7 +311,6 @@ CONTAINS ...@@ -73,7 +311,6 @@ CONTAINS
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices 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(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 REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
! !
...@@ -89,7 +326,7 @@ CONTAINS ...@@ -89,7 +326,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection' IF(lwp) WRITE(numout,*) 'dyn_adv_ubs_hls1 : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF ENDIF
ENDIF ENDIF
...@@ -230,63 +467,44 @@ CONTAINS ...@@ -230,63 +467,44 @@ CONTAINS
! ! Vertical advection ! ! ! Vertical advection !
! ! ==================== ! ! ! ==================== !
! !
! ! ======================== ! DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
IF( PRESENT( no_zad ) ) THEN ! No vertical advection ! (except if linear free surface) zfu_uw(ji,jj,jpk) = 0._wp
! ! ======================== ! ------ zfv_vw(ji,jj,jpk) = 0._wp
! zfu_uw(ji,jj, 1 ) = 0._wp
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0 zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 0, 0, 0 ) END_2D
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) IF( ln_linssh ) THEN ! constant volume : advection through the surface
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) DO_2D( 0, 0, 0, 0 )
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) & 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)
& / e3u(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)
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_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 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
! !
END SUBROUTINE dyn_adv_ubs 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' )
!
END SUBROUTINE dyn_adv_ubs_hls1
!!============================================================================== !!==============================================================================
END MODULE dynadv_ubs END MODULE dynadv_ubs
src/OCE/DYN/dynatf.F90
View file @ 005ca9c1
...@@ -328,29 +328,29 @@ CONTAINS ...@@ -328,29 +328,29 @@ CONTAINS
! !
IF ( iom_use("utau") ) THEN IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj)) ALLOCATE(zutau(A2D(0)))
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj) jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa) zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D END_2D
CALL iom_put( "utau", zutau(:,:) ) CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau) DEALLOCATE(zutau)
ELSE ELSE
CALL iom_put( "utau", utau(:,:) ) CALL iom_put( "utau", utau(A2D(0)) )
ENDIF ENDIF
ENDIF ENDIF
! !
IF ( iom_use("vtau") ) THEN IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj)) ALLOCATE(zvtau(A2D(0)))
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj) jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa) zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D END_2D
CALL iom_put( "vtau", zvtau(:,:) ) CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau) DEALLOCATE(zvtau)
ELSE ELSE
CALL iom_put( "vtau", vtau(:,:) ) CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF ENDIF
ENDIF ENDIF
! !
......
src/OCE/DYN/dynatf_qco.F90
View file @ 005ca9c1
...@@ -245,29 +245,29 @@ CONTAINS ...@@ -245,29 +245,29 @@ CONTAINS
! !
IF ( iom_use("utau") ) THEN IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj)) ALLOCATE(zutau(A2D(0)))
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj) jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa) zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D END_2D
CALL iom_put( "utau", zutau(:,:) ) CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau) DEALLOCATE(zutau)
ELSE ELSE
CALL iom_put( "utau", utau(:,:) ) CALL iom_put( "utau", utau(A2D(0)) )
ENDIF ENDIF
ENDIF ENDIF
! !
IF ( iom_use("vtau") ) THEN IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj)) ALLOCATE(zvtau(A2D(0)))
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj) jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa) zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D END_2D
CALL iom_put( "vtau", zvtau(:,:) ) CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau) DEALLOCATE(zvtau)
ELSE ELSE
CALL iom_put( "vtau", vtau(:,:) ) CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF ENDIF
ENDIF ENDIF
! !
......
src/OCE/DYN/dynkeg.F90
View file @ 005ca9c1
...@@ -6,7 +6,8 @@ MODULE dynkeg ...@@ -6,7 +6,8 @@ MODULE dynkeg
!! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code !! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
!! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
!! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module !! 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 ...@@ -27,7 +28,8 @@ MODULE dynkeg
IMPLICIT NONE IMPLICIT NONE
PRIVATE 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_C2 = 0 !: 2nd order centered scheme (standard scheme)
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983 INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
...@@ -44,6 +46,123 @@ MODULE dynkeg ...@@ -44,6 +46,123 @@ MODULE dynkeg
CONTAINS CONTAINS
SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs ) 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 *** !! *** ROUTINE dyn_keg ***
!! !!
...@@ -86,7 +205,7 @@ CONTAINS ...@@ -86,7 +205,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*) 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,*) '~~~~~~~' IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF ENDIF
ENDIF ENDIF
...@@ -147,7 +266,7 @@ CONTAINS ...@@ -147,7 +266,7 @@ CONTAINS
! !
IF( ln_timing ) CALL timing_stop('dyn_keg') IF( ln_timing ) CALL timing_stop('dyn_keg')
! !
END SUBROUTINE dyn_keg END SUBROUTINE dyn_keg_hls1
!!====================================================================== !!======================================================================
END MODULE dynkeg END MODULE dynkeg
src/OCE/DYN/dynspg_ts.F90
View file @ 005ca9c1
...@@ -334,14 +334,14 @@ CONTAINS ...@@ -334,14 +334,14 @@ CONTAINS
! ! ------------------ ! ! ! ------------------ !
IF( ln_bt_fw ) THEN IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D END_2D
ELSE ELSE
zztmp = r1_rho0 * r1_2 zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 ) DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
END_2D END_2D
ENDIF ENDIF
! !
......
src/OCE/DYN/dynvor.F90
View file @ 005ca9c1
...@@ -22,6 +22,7 @@ MODULE dynvor ...@@ -22,6 +22,7 @@ MODULE dynvor
!! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation
!! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity !! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity
!! 4.2 ! 2020-12 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T) !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -172,11 +173,20 @@ CONTAINS ...@@ -172,11 +173,20 @@ CONTAINS
CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF ENDIF
CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
IF( nn_hls == 1 ) THEN
CALL vor_eeT_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_eeT_hls1( 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_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ELSE
CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 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 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 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 CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ENDIF ENDIF
CASE( np_ENE ) !* energy conserving scheme CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
...@@ -199,11 +209,20 @@ CONTAINS ...@@ -199,11 +209,20 @@ CONTAINS
IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 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 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 CASE( np_EEN ) !* energy and enstrophy conserving scheme
IF( nn_hls == 1 ) THEN
CALL vor_een_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_een_hls1( 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_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ELSE
CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 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 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 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 CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ENDIF ENDIF
END SELECT END SELECT
! !
...@@ -320,7 +339,122 @@ CONTAINS ...@@ -320,7 +339,122 @@ CONTAINS
! !
INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwt ! 2D workspace REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
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:vor_enT : vorticity term: t-point energy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
SELECT CASE( kvor ) !== relative vorticity considered ==!
!
CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
ALLOCATE( zwz(A2D(1)) )
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
!
END SELECT
!
SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
& + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) &
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_MET ) !* metric term
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
& * e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
& + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) &
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
& + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
& * e3t(ji,jj,jk,Kmm)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
END SELECT
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 0, 0, 0 )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
& + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
!
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
& + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
END_2D
! ! ===============
END DO ! End of slab
! ! ===============
!
SELECT CASE( kvor ) ! deallocate zwz if necessary
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT
!
END SUBROUTINE vor_enT
SUBROUTINE vor_enT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_enT ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and t-point evaluation of vorticity (planetary and relative).
!! conserves the horizontal kinetic energy.
!! The general trend of momentum is increased due to the vorticity
!! term which is given by:
!! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ]
!! vorv = 1/bv mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f mj[un] ]
!! where rvor is the relative vorticity at f-point
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
...@@ -336,7 +470,7 @@ CONTAINS ...@@ -336,7 +470,7 @@ CONTAINS
SELECT CASE( kvor ) !== relative vorticity considered ==! SELECT CASE( kvor ) !== relative vorticity considered ==!
! !
CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
ALLOCATE( zwz(A2D(nn_hls),jpk) ) ALLOCATE( zwz(A2D(1),jpk) )
DO jk = 1, jpkm1 ! Horizontal slab DO jk = 1, jpkm1 ! Horizontal slab
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
...@@ -409,7 +543,7 @@ CONTAINS ...@@ -409,7 +543,7 @@ CONTAINS
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz ) CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT END SELECT
! !
END SUBROUTINE vor_enT END SUBROUTINE vor_enT_hls1
SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs ) SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
...@@ -440,7 +574,7 @@ CONTAINS ...@@ -440,7 +574,7 @@ CONTAINS
! !
INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
...@@ -573,7 +707,7 @@ CONTAINS ...@@ -573,7 +707,7 @@ CONTAINS
! !
INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
...@@ -705,16 +839,176 @@ CONTAINS ...@@ -705,16 +839,176 @@ CONTAINS
INTEGER :: ierr ! local integer INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, ze3f ! local scalars REAL(wp) :: zmsk, ze3f ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: z1_e3f REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
#if defined key_loop_fusion #if defined key_loop_fusion
REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1 REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1 REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1 REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
#else #else
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#endif #endif
REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
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:vor_een : vorticity term: energy and enstrophy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
#if defined key_qco || defined key_linssh
DO_2D( 1, 1, 1, 1 ) ! == reciprocal of e3 at F-point (key_qco)
z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
END_2D
#else
SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
& + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
ELSE ; z1_e3f(ji,jj) = 0._wp
ENDIF
END_2D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
& + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
ELSE ; z1_e3f(ji,jj) = 0._wp
ENDIF
END_2D
END SELECT
#endif
!
SELECT CASE( kvor ) !== vorticity considered ==!
!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
CASE ( np_MET ) !* metric term
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
END_2D
ENDIF
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
END SELECT
!
! !== horizontal fluxes ==!
DO_2D( 1, 1, 1, 1 )
zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
END_2D
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 1, 0, 1 )
ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
END_2D
!
DO_2D( 0, 0, 0, 0 )
zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
END_2D
END DO
! ! ===============
! ! End of slab
! ! ===============
END SUBROUTINE vor_een
SUBROUTINE vor_een_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_een ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and the Arakawa and Lamb (1980) flux form formulation : conserves
!! both the horizontal kinetic energy and the potential enstrophy
!! when horizontal divergence is zero (see the NEMO documentation)
!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!
!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, ze3f ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
#if defined key_loop_fusion
REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
#else
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#endif
REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
...@@ -874,7 +1168,7 @@ CONTAINS ...@@ -874,7 +1168,7 @@ CONTAINS
! ! =============== ! ! ===============
! ! End of slab ! ! End of slab
! ! =============== ! ! ===============
END SUBROUTINE vor_een END SUBROUTINE vor_een_hls1
SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs ) SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
...@@ -904,9 +1198,129 @@ CONTAINS ...@@ -904,9 +1198,129 @@ CONTAINS
INTEGER :: ierr ! local integer INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars REAL(wp) :: zmsk, z1_e3t ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
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:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
!
SELECT CASE( kvor ) !== vorticity considered ==!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
& - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
& * r1_e1e2f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
CASE ( np_MET ) !* metric term
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( ff_f(ji,jj) + ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
& - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
& * r1_e1e2f(ji,jj) )
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
END_2D
ENDIF
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
END SELECT
!
!
! !== horizontal fluxes ==!
DO_2D( 1, 1, 1, 1 )
zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
END_2D
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 1, 0, 1 )
z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t
ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t
ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t
END_2D
!
DO_2D( 0, 0, 0, 0 )
zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
END_2D
! ! ===============
END DO ! End of slab
! ! ===============
END SUBROUTINE vor_eeT
SUBROUTINE vor_eeT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_eeT ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and the Arakawa and Lamb (1980) vector form formulation using
!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
!! The change consists in
!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!
!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
...@@ -1003,7 +1417,7 @@ CONTAINS ...@@ -1003,7 +1417,7 @@ CONTAINS
! ! =============== ! ! ===============
END DO ! End of slab END DO ! End of slab
! ! =============== ! ! ===============
END SUBROUTINE vor_eeT END SUBROUTINE vor_eeT_hls1
SUBROUTINE dyn_vor_init SUBROUTINE dyn_vor_init
......
src/OCE/DYN/dynzdf.F90
View file @ 005ca9c1
...@@ -6,6 +6,7 @@ MODULE dynzdf ...@@ -6,6 +6,7 @@ MODULE dynzdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code !! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -79,7 +80,7 @@ CONTAINS ...@@ -79,7 +80,7 @@ CONTAINS
REAL(wp) :: zWu , zWv ! - - REAL(wp) :: zWu , zWv ! - -
REAL(wp) :: zWui, zWvi ! - - REAL(wp) :: zWui, zWvi ! - -
REAL(wp) :: zWus, zWvs ! - - REAL(wp) :: zWus, zWvs ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwd, zws ! 3D workspace REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwd, zws ! 2D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
!!--------------------------------------------------------------------- !!---------------------------------------------------------------------
! !
...@@ -105,315 +106,329 @@ CONTAINS ...@@ -105,315 +106,329 @@ CONTAINS
ztrdv(:,:,:) = pvv(:,:,:,Krhs) ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF ENDIF
! !
! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) ! ! ================= !
! DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! time stepping except vertical diffusion ! ! ================= !
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity !
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk) !
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk) ! ! time stepping except vertical diffusion
END_3D IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
ELSE ! applied on thickness weighted velocity DO_2Dik( 0, 0, 1, jpkm1, 1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) & pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
& + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) & END_2D
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) ELSE ! applied on thickness weighted velocity
pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) & DO_2Dik( 0, 0, 1, jpkm1, 1 )
& + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) & puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
END_3D & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
ENDIF pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
! ! add top/bottom friction & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom. & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does END_2D
! not lead to the effective stress seen over the whole barotropic loop. ENDIF
! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) ! ! add top/bottom friction
IF( ln_drgimp .AND. ln_dynspg_ts ) THEN ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk) ! not lead to the effective stress seen over the whole barotropic loop.
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
END_3D IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
DO_2D( 0, 0, 0, 0 ) ! Add bottom/top stress due to barotropic component only DO_2Dik( 0, 0, 1, jpkm1, 1 ) ! remove barotropic velocities
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) & END_2D
& / e3u(ji,jj,iku,Kaa) DO_1Di( 0, 0 ) ! Add bottom/top stress due to barotropic component only
pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) & iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
& / e3v(ji,jj,ikv,Kaa) ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
END_2D puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) &
IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
DO_2D( 0, 0, 0, 0 )
iku = miku(ji,jj) ! top ocean level at u- and v-points
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
& / e3u(ji,jj,iku,Kaa) & / e3u(ji,jj,iku,Kaa)
pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) & pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
& / e3v(ji,jj,ikv,Kaa) & / e3v(ji,jj,ikv,Kaa)
END_2D END_1D
END IF IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
ENDIF DO_1Di( 0, 0 )
! iku = miku(ji,jj) ! top ocean level at u- and v-points
! !== Vertical diffusion on u ==! ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
! puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
! !* Matrix construction & / e3u(ji,jj,iku,Kaa)
IF( ln_zad_Aimp ) THEN !! pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) &
SELECT CASE( nldf_dyn ) & / e3v(ji,jj,ikv,Kaa)
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) END_1D
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) END IF
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point ENDIF
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & !
& / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk ) ! !== Vertical diffusion on u ==!
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & !
& / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1) !
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua ! !* Matrix construction
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua IF( ln_zad_Aimp ) THEN !- including terms associated with partly implicit vertical advection
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) SELECT CASE( nldf_dyn )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) DO_2Dik( 0, 0, 1, jpkm1, 1 )
END_3D z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
CASE DEFAULT ! iso-level lateral mixing zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
& / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk ) zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
& / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1) zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) END_2D
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) CASE DEFAULT ! iso-level lateral mixing
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) DO_2Dik( 0, 0, 1, jpkm1, 1 )
END_3D z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
END SELECT zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
zwi(ji,jj,1) = 0._wp zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
& / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa) zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp ) zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) ) zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
END_2D zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
ELSE END_2D
SELECT CASE( nldf_dyn ) END SELECT
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) !
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) zwi(:,1) = 0._wp
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & DO_1Di( 0, 0 ) !* Surface boundary conditions
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) zwi(ji,1) = 0._wp
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) & / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
zwi(ji,jj,jk) = zzwi zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
zws(ji,jj,jk) = zzws zws(ji,1) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws zwd(ji,1) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
END_3D END_1D
CASE DEFAULT ! iso-level lateral mixing ELSE !- only vertical diffusive terms
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) SELECT CASE( nldf_dyn )
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zwi(ji,jj,jk) = zzwi zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
zws(ji,jj,jk) = zzws & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
zwd(ji,jj,jk) = 1._wp - zzwi - zzws zwi(ji,jk) = zzwi
END_3D zws(ji,jk) = zzws
END SELECT zwd(ji,jk) = 1._wp - zzwi - zzws
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions END_2D
zwi(ji,jj,1) = 0._wp CASE DEFAULT ! iso-level lateral mixing
zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) DO_2Dik( 0, 0, 1, jpkm1, 1 )
END_2D zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
ENDIF & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
! zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
! & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
! !== Apply semi-implicit bottom friction ==! zwi(ji,jk) = zzwi
! zws(ji,jk) = zzws
! Only needed for semi-implicit bottom friction setup. The explicit zwd(ji,jk) = 1._wp - zzwi - zzws
! bottom friction has been included in "u(v)a" which act as the R.H.S END_2D
! column vector of the tri-diagonal matrix equation END SELECT
! !
IF ( ln_drgimp ) THEN ! implicit bottom friction zwi(:,1) = 0._wp
DO_2D( 0, 0, 0, 0 ) DO_1Di( 0, 0 ) !* Surface boundary conditions
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points zwd(ji,1) = 1._wp - zws(ji,1)
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) & END_1D
& / e3u(ji,jj,iku,Kaa) ENDIF
END_2D !
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit) !
DO_2D( 0, 0, 0, 0 ) ! !== Apply semi-implicit bottom friction ==!
!!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed !
iku = miku(ji,jj) ! ocean top level at u- and v-points ! Only needed for semi-implicit bottom friction setup. The explicit
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) & ! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
IF ( ln_drgimp ) THEN ! implicit bottom friction
DO_1Di( 0, 0 )
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa) & / e3u(ji,jj,iku,Kaa)
END_2D END_1D
END IF IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
ENDIF DO_1Di( 0, 0 )
! !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
! Matrix inversion starting from the first level iku = miku(ji,jj) ! ocean top level at u- and v-points
!----------------------------------------------------------------------- zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) & / e3u(ji,jj,iku,Kaa)
! END_1D
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) ENDIF
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) ENDIF
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) !
! ( ... )( ... ) ( ... ) ! Matrix inversion starting from the first level
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) !-----------------------------------------------------------------------
! ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
! m is decomposed in the product of an upper and a lower triangular matrix !
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! The solution (the after velocity) is in puu(:,:,:,Kaa) ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
!----------------------------------------------------------------------- ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ! ( ... )( ... ) ( ... )
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) !
END_3D ! m is decomposed in the product of an upper and a lower triangular matrix
! ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! ! The solution (the after velocity) is in puu(:,:,:,Kaa)
!-----------------------------------------------------------------------
!
DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
END_2D
!
DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3 #if defined key_RK3
! ! RK3: use only utau (not utau_b) ! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj) & puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1) & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else #else
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) & ! ! MLF: average of utau and utau_b
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1) puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif #endif
END_2D END_1D
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) DO_2Dik( 0, 0, 2, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * puu(ji,jj,jk-1,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
END_3D
!
! !== Vertical diffusion on v ==!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,1) = 0._wp
zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
& / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
END_2D
ELSE
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,1) = 0._wp
zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
END_2D END_2D
ENDIF !
! DO_1Di( 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
! !== Apply semi-implicit top/bottom friction ==! puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
! END_1D
! Only needed for semi-implicit bottom friction setup. The explicit DO_2Dik( 0, 0, jpk-2, 1, -1 )
! bottom friction has been included in "u(v)a" which act as the R.H.S puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
! column vector of the tri-diagonal matrix equation
!
IF( ln_drgimp ) THEN
DO_2D( 0, 0, 0, 0 )
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
END_2D END_2D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN !
DO_2D( 0, 0, 0, 0 ) !
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) ! !== Vertical diffusion on v ==!
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) & !
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_2D
CASE DEFAULT ! iso-level lateral mixing
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_2D
END SELECT
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
& / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
zws(ji,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
END_1D
ELSE
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
CASE DEFAULT ! iso-level lateral mixing
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
END SELECT
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zwd(ji,1) = 1._wp - zws(ji,1)
END_1D
ENDIF
!
! !== Apply semi-implicit top/bottom friction ==!
!
! Only needed for semi-implicit bottom friction setup. The explicit
! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
IF( ln_drgimp ) THEN
DO_1Di( 0, 0 )
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa) & / e3v(ji,jj,ikv,Kaa)
END_2D END_1D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
DO_1Di( 0, 0 )
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
ENDIF
ENDIF ENDIF
ENDIF
! Matrix inversion
! Matrix inversion !-----------------------------------------------------------------------
!----------------------------------------------------------------------- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) !
! ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) ! ( ... )( ... ) ( ... )
! ( ... )( ... ) ( ... ) ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) !
! ! m is decomposed in the product of an upper and lower triangular matrix
! m is decomposed in the product of an upper and lower triangular matrix ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi ! The solution (after velocity) is in 2d array va
! The solution (after velocity) is in 2d array va !-----------------------------------------------------------------------
!----------------------------------------------------------------------- !
! DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) END_2D
END_3D !
! DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3 #if defined key_RK3
! ! RK3: use only vtau (not vtau_b) ! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj) & pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1) & / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else #else
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) & ! ! MLF: average of vtau and vtau_b
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1) pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtauV(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif #endif
END_2D END_1D
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) DO_2Dik( 0, 0, 2, jpkm1, 1 )
pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * pvv(ji,jj,jk-1,Kaa)
END_3D END_2D
! !
DO_2D( 0, 0, 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! DO_1Di( 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
END_2D END_1D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) DO_2Dik( 0, 0, jpk-2, 1, -1 )
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
END_3D END_2D
! ! ================= !
END_1D ! i-k slices loop !
! ! ================= !
! !
IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:) ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
......
src/OCE/DYN/sshwzv.F90
View file @ 005ca9c1
...@@ -244,28 +244,28 @@ CONTAINS ...@@ -244,28 +244,28 @@ CONTAINS
! inside computational domain (cosmetic) ! inside computational domain (cosmetic)
DO jk = 1, jpkm1 DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- ! IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls) DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls)
DO jj = 1, jpj DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_east ) THEN ! --- East --- ! IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_south ) THEN ! --- South --- ! IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls) DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_north ) THEN ! --- North --- ! IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
...@@ -375,28 +375,28 @@ CONTAINS ...@@ -375,28 +375,28 @@ CONTAINS
! inside computational domain (cosmetic) ! inside computational domain (cosmetic)
DO jk = 1, jpkm1 DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- ! IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls) DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls)
DO jj = 1, jpj DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_east ) THEN ! --- East --- ! IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_south ) THEN ! --- South --- ! IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls) DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
END DO END DO
ENDIF ENDIF
IF( lk_north ) THEN ! --- North --- ! IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp pww(ji,jj,jk) = 0._wp
END DO END DO
......
src/OCE/FLO/floblk.F90
View file @ 005ca9c1
...@@ -105,10 +105,10 @@ CONTAINS ...@@ -105,10 +105,10 @@ CONTAINS
iloop = 0 iloop = 0
222 DO jfl = 1, jpnfl 222 DO jfl = 1, jpnfl
# if ! defined key_mpi_off # if ! defined key_mpi_off
IF( iil(jfl) >= mig(Nis0) .AND. iil(jfl) <= mig(Nie0) .AND. & IF( iil(jfl) >= mig(Nis0,nn_hls) .AND. iil(jfl) <= mig(Nie0,nn_hls) .AND. &
ijl(jfl) >= mjg(Njs0) .AND. ijl(jfl) <= mjg(Nje0) ) THEN ijl(jfl) >= mjg(Njs0,nn_hls) .AND. ijl(jfl) <= mjg(Nje0,nn_hls) ) THEN
iiloc(jfl) = iil(jfl) - mig(1) + 1 iiloc(jfl) = iil(jfl) - mig(1,nn_hls) + 1
ijloc(jfl) = ijl(jfl) - mjg(1) + 1 ijloc(jfl) = ijl(jfl) - mjg(1,nn_hls) + 1
# else # else
iiloc(jfl) = iil(jfl) iiloc(jfl) = iil(jfl)
ijloc(jfl) = ijl(jfl) ijloc(jfl) = ijl(jfl)
......