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 1073 additions and 603 deletions
src/OCE/DYN/dynspg_ts.F90
View file @ 4dbfc9bb
......@@ -334,14 +334,14 @@ CONTAINS
! ! ------------------ !
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(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)
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 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
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)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(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) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
END_2D
ENDIF
!
......@@ -459,8 +459,6 @@ CONTAINS
DO jn = 1, icycle ! sub-time-step loop !
! ! ==================== !
!
l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1
!
! !== Update the forcing ==! (BDY and tides)
!
IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) )
......@@ -1249,16 +1247,16 @@ CONTAINS
SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( 0, 0, 0, 0 )
zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
& + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp
zwz(ji,jj) = ( ( ht(ji,jj+1) + ht(ji+1,jj+1) ) & ! need additional () for reproducibility around NP
& + ( ht(ji,jj ) + ht(ji+1,jj ) ) ) * 0.25_wp
IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
END_2D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_2D( 0, 0, 0, 0 )
zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
& + ht(ji,jj ) + ht(ji+1,jj ) ) &
& / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
& + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) )
zwz(ji,jj) = ( ( ht(ji,jj+1) + ht(ji+1,jj+1) ) & ! need additional () for
& + ( ht(ji,jj ) + ht(ji+1,jj ) ) ) & ! reproducibility around NP
& / ( MAX( ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
& + ssmask(ji,jj ) + ssmask(ji+1,jj ), 1._wp ) )
IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
END_2D
END SELECT
......@@ -1329,24 +1327,24 @@ CONTAINS
!
CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
DO_2D( 0, 0, 0, 0 )
zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
zy1 = r1_8 * ( ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) ) & ! need additional () for
& + ( zhV(ji ,jj ) + zhV(ji+1,jj ) ) ) * r1_e1u(ji,jj) ! reproducibility around NP
zx1 = - r1_8 * ( ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) ) & ! need additional () for
& + ( zhU(ji ,jj ) + zhU(ji ,jj+1) ) ) * r1_e2v(ji,jj) ! reproducibility around NP
zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
END_2D
!
CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
DO_2D( 0, 0, 0, 0 )
zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
& + ftse(ji,jj ) * zhV(ji ,jj-1) &
& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
& + ftne(ji,jj ) * zhU(ji ,jj ) )
zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ( ftne(ji,jj ) * zhV(ji ,jj ) & ! need additional () for
& + ftnw(ji+1,jj) * zhV(ji+1,jj ) ) & ! reproducibility around NP
& + ( ftse(ji,jj ) * zhV(ji ,jj-1) &
& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) )
zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
& + ftse(ji,jj+1) * zhU(ji ,jj+1) ) &
& + ( ftnw(ji,jj ) * zhU(ji-1,jj ) &
& + ftne(ji,jj ) * zhU(ji ,jj ) ) )
END_2D
!
END SELECT
......
src/OCE/DYN/dynvor.F90
View file @ 4dbfc9bb
......@@ -22,6 +22,7 @@ MODULE dynvor
!! - ! 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.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
CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
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
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
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
ENDIF
ENDIF
CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
......@@ -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_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
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
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
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
ENDIF
ENDIF
END SELECT
!
......@@ -320,7 +339,122 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
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) ) & ! add () for
& - ( e1u(ji,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) * r1_e1e2f(ji,jj) ! NP repro
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 ) ) & ! need additional () for
& + ( zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) & ! reproducibility around NP
& * 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 ) ) & ! need additional () for
& + ( zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) ) & ! reproducibility around NP
& * 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
!!----------------------------------------------------------------------
!
......@@ -336,11 +470,11 @@ CONTAINS
SELECT CASE( kvor ) !== relative vorticity considered ==!
!
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_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) &
& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
zwz(ji,jj,jk) = ( ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! need additional () for NP repro
& - ( 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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
......@@ -364,8 +498,8 @@ CONTAINS
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
zwt(ji,jj) = r1_4 * ( ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) ) & ! need additional () for
& + ( zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) & ! reproducibility around NP
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_MET ) !* metric term
......@@ -376,8 +510,8 @@ CONTAINS
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 ,jk) + zwz(ji,jj ,jk) &
& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) &
zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) ) & ! need additional () for
& + ( zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) ) & ! reproducibility around NP
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_CME ) !* Coriolis + metric
......@@ -409,7 +543,7 @@ CONTAINS
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT
!
END SUBROUTINE vor_enT
END SUBROUTINE vor_enT_hls1
SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......@@ -440,7 +574,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
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
......@@ -462,8 +596,8 @@ CONTAINS
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 0, 1, 0 )
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)
zwz(ji,jj) = ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for
& - ( e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) * r1_e1e2f(ji,jj) ! NP repro
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 0, 1, 0 )
......@@ -477,8 +611,8 @@ CONTAINS
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 0, 1, 0 )
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)
zwz(ji,jj) = ff_f(ji,jj) + ( ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for NP repro
& - ( 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 (NOT the Coriolis term)
DO_2D( 1, 0, 1, 0 )
......@@ -502,22 +636,22 @@ CONTAINS
SELECT CASE( nn_e3f_typ ) !== potential vorticity ==!
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( 1, 0, 1, 0 )
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) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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 ; zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
ELSE ; zwz(ji,jj) = 0._wp
ENDIF
END_2D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_2D( 1, 0, 1, 0 )
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) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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 ; zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
ELSE ; zwz(ji,jj) = 0._wp
ENDIF
......@@ -573,7 +707,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
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
......@@ -594,8 +728,8 @@ CONTAINS
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 0, 1, 0 )
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)
zwz(ji,jj) = ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for NP repro
& - ( 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, 0, 1, 0 )
......@@ -609,8 +743,8 @@ CONTAINS
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 0, 1, 0 )
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)
zwz(ji,jj) = ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) &! add () for NP repro
& -( 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 (NOT the Coriolis term)
DO_2D( 1, 0, 1, 0 )
......@@ -635,22 +769,22 @@ CONTAINS
SELECT CASE( nn_e3f_typ ) !== potential vorticity ==!
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( 1, 0, 1, 0 )
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) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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 ; zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
ELSE ; zwz(ji,jj) = 0._wp
ENDIF
END_2D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_2D( 1, 0, 1, 0 )
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) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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 ; zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
ELSE ; zwz(ji,jj) = 0._wp
ENDIF
......@@ -665,10 +799,10 @@ CONTAINS
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 0, 0, 0 )
zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
zuav = r1_8 * r1_e1u(ji,jj) * ( ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) & ! need additional () for
& + ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) ) ! reproducibility around NP
zvau =-r1_8 * r1_e2v(ji,jj) * ( ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) &
& + ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
END_2D
......@@ -705,16 +839,168 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! 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
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)) :: 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 )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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 )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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) ) & ! add () for NP repro
& - ( 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 )
zwz(ji,jj) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for
& - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) & ! NP repro
& ) * 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 ) ! add () for NP repro
zua = + r1_12 * r1_e1u(ji,jj) * ( ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) ) & ! add () for
& + ( ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) ) ! NP repro
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(nn_hls)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#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),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
......@@ -737,26 +1023,22 @@ CONTAINS
SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-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)) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-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) )
ze3f = ( ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( 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
......@@ -772,8 +1054,8 @@ CONTAINS
END_2D
CASE ( np_RVO ) !* relative vorticity
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) &
& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
zwz(ji,jj,jk) = ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for NP repro
& - ( 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( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
......@@ -787,12 +1069,8 @@ CONTAINS
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj,jk) = ( 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
zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for
& - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) & ! NP repro
& ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
......@@ -837,10 +1115,10 @@ CONTAINS
ztsw_jp1 = zwz(ji ,jj ,jk) + zwz(ji-1,jj ,jk) + zwz(ji-1,jj+1,jk)
ztsw_ip1 = zwz(ji+1,jj-1,jk) + zwz(ji ,jj-1,jk) + zwz(ji ,jj ,jk)
!
zua = + r1_12 * r1_e1u(ji,jj) * ( ztne * zwy + ztnw_ip1 * zwy_ip1 &
& + ztse * zwy_jm1 + ztsw_ip1 * zwy_ip1_jm1 )
zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw_jp1 * zwx_im1_jp1 + ztse_jp1 * zwx_jp1 &
& + ztnw * zwx_im1 + ztne * zwx )
zua = + r1_12 * r1_e1u(ji,jj) * ( ( ztne * zwy + ztnw_ip1 * zwy_ip1 ) & ! need additional () for
& + ( ztse * zwy_jm1 + ztsw_ip1 * zwy_ip1_jm1 ) ) ! reproducibility around NP
zva = - r1_12 * r1_e2v(ji,jj) * ( ( ztsw_jp1 * zwx_im1_jp1 + ztse_jp1 * zwx_jp1 ) &
& + ( ztnw * zwx_im1 + ztne * zwx ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
END_3D
......@@ -862,10 +1140,10 @@ CONTAINS
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 ) )
zua = + r1_12 * r1_e1u(ji,jj) * ( ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) ) & ! add () for
& + ( ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) ) ! NP repro
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
......@@ -874,7 +1152,7 @@ CONTAINS
! ! ===============
! ! End of slab
! ! ===============
END SUBROUTINE vor_een
END SUBROUTINE vor_een_hls1
SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......@@ -904,9 +1182,125 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: 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)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
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 )
zwz(ji,jj) = ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for
& - ( e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) & ! NP reproducibility
& * 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 )
zwz(ji,jj) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! add () for
& - ( e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) & ! NP repro
& * 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 ) ) & ! add () for
& + ( ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) ) ! NP repro
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
......@@ -929,10 +1323,8 @@ CONTAINS
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj,jk) = ( (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)) ) &
zwz(ji,jj,jk) = ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) & ! need additional ()
& - ( e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) & ! for NP reproducibility
& * r1_e1e2f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
......@@ -947,11 +1339,9 @@ CONTAINS
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj,jk) = ( 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) )
zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) ) &! add () for
& - ( e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) ) ) &! NP repro
& * r1_e1e2f(ji,jj) )
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
......@@ -993,17 +1383,17 @@ CONTAINS
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 ) )
zua = + r1_12 * r1_e1u(ji,jj) * ( ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) ) & ! add () for
& + ( ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) ) ! NP repro
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
END SUBROUTINE vor_eeT_hls1
SUBROUTINE dyn_vor_init
......@@ -1115,10 +1505,10 @@ CONTAINS
SELECT CASE( nn_e3f_typ )
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_3D( 0, 0, 0, 0, 1, jpk )
e3f_0vor(ji,jj,jk) = ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) &
& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &
& + e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
& + e3t_0(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25_wp
e3f_0vor(ji,jj,jk) = ( ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
& + e3t_0(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) ) * 0.25_wp
END_3D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_3D( 0, 0, 0, 0, 1, jpk )
......@@ -1126,10 +1516,10 @@ CONTAINS
& + tmask(ji,jj ,jk) +tmask(ji+1,jj ,jk) )
!
IF( zmsk /= 0._wp ) THEN
e3f_0vor(ji,jj,jk) = ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) &
& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &
& + e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
& + e3t_0(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) / zmsk
e3f_0vor(ji,jj,jk) = ( ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) & ! need additional () for
& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) ) & ! reproducibility around NP
& + ( e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
& + e3t_0(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) ) / zmsk
ELSE ; e3f_0vor(ji,jj,jk) = 0._wp
ENDIF
END_3D
......
src/OCE/DYN/dynzdf.F90
View file @ 4dbfc9bb
......@@ -6,6 +6,7 @@ MODULE dynzdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 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.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -79,7 +80,7 @@ CONTAINS
REAL(wp) :: zWu , zWv ! - -
REAL(wp) :: zWui, zWvi ! - -
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 ! - -
!!---------------------------------------------------------------------
!
......@@ -105,315 +106,329 @@ CONTAINS
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
!
! ! time stepping except vertical diffusion
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
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)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
END_3D
ELSE ! applied on thickness weighted velocity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
& + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
& + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
END_3D
ENDIF
! ! add top/bottom friction
! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
! not lead to the effective stress seen over the whole barotropic loop.
! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
END_3D
DO_2D( 0, 0, 0, 0 ) ! Add bottom/top stress due to barotropic component only
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
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) &
& / e3u(ji,jj,iku,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)
END_2D
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) &
! ! ================= !
DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! ================= !
!
! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
!
! ! time stepping except vertical diffusion
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
DO_2Dik( 0, 0, 1, jpkm1, 1 )
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)
END_2D
ELSE ! applied on thickness weighted velocity
DO_2Dik( 0, 0, 1, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
& + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
& + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
END_2D
ENDIF
! ! add top/bottom friction
! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
! not lead to the effective stress seen over the whole barotropic loop.
! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
DO_2Dik( 0, 0, 1, jpkm1, 1 ) ! remove barotropic velocities
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
END_2D
DO_1Di( 0, 0 ) ! Add bottom/top stress due to barotropic component only
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
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) &
& / 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)
END_2D
END IF
ENDIF
!
! !== Vertical diffusion on u ==!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
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 )
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
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(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
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 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+1,jj,2) + avm(ji ,jj,2) ) &
& / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 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+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
& / ( 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) + akzu(ji,jj,jk+1) ) &
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(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+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / ( 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)
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
ENDIF
!
!
! !== Apply semi-implicit 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 ! implicit bottom friction
DO_2D( 0, 0, 0, 0 )
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa)
END_2D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
DO_2D( 0, 0, 0, 0 )
!!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
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
END_1D
IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
DO_1Di( 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)
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) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
END IF
ENDIF
!
! !== Vertical diffusion on u ==!
!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !- including terms associated with partly implicit vertical advection
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
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 )
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
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_2D
CASE DEFAULT ! iso-level lateral mixing
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(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
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_2D
END SELECT
!
zwi(:,1) = 0._wp
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
& / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
zws(ji,1) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,1) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
END_1D
ELSE !- only vertical diffusive terms
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+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
& / ( 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) + akzu(ji,jj,jk+1) ) &
& / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(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+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / ( 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)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
END SELECT
!
zwi(:,1) = 0._wp
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwd(ji,1) = 1._wp - zws(ji,1)
END_1D
ENDIF
!
!
! !== Apply semi-implicit 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 ! 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)
END_2D
END IF
ENDIF
!
! Matrix inversion starting from the first level
!-----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and a lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! The solution (the after velocity) is in puu(:,:,:,Kaa)
!-----------------------------------------------------------------------
!
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
END_1D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
DO_1Di( 0, 0 )
!!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
zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa)
END_1D
ENDIF
ENDIF
!
! Matrix inversion starting from the first level
!-----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and a lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! 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
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
! ! MLF: average of utau and utau_b
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
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
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)
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_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * puu(ji,jj,jk-1,Kaa)
END_2D
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_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)
!
DO_1Di( 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
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)
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
!
!
! !== 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_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)
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
! Matrix inversion
!-----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! The solution (after velocity) is in 2d array va
!-----------------------------------------------------------------------
!
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
! Matrix inversion
!-----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! 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) ==
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
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
! ! MLF: average of vtau and vtau_b
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
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
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)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 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)
END_3D
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * pvv(ji,jj,jk-1,Kaa)
END_2D
!
DO_1Di( 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
END_2D
! ! ================= !
END_1D ! i-k slices loop !
! ! ================= !
!
IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
......
src/OCE/DYN/sshwzv.F90
View file @ 4dbfc9bb
......@@ -175,7 +175,8 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~~~'
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
pww(:,:,:) = 0._wp ! bottom boundary condition: w=0 (set once for all)
! ! needed over the halos for the output (ww+wi) in diawri.F90
ENDIF
! !------------------------------!
! ! Now Vertical Velocity !
......@@ -244,28 +245,28 @@ CONTAINS
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
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
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
......@@ -314,7 +315,8 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
pww(:,:,:) = 0._wp ! bottom boundary condition: w=0 (set once for all)
! ! needed over the halos for the output (ww+wi) in diawri.F90
ENDIF
!
CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
......@@ -375,28 +377,28 @@ CONTAINS
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
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
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
......
src/OCE/FLO/floblk.F90
View file @ 4dbfc9bb
......@@ -105,10 +105,10 @@ CONTAINS
iloop = 0
222 DO jfl = 1, jpnfl
# if ! defined key_mpi_off
IF( iil(jfl) >= mig(Nis0) .AND. iil(jfl) <= mig(Nie0) .AND. &
ijl(jfl) >= mjg(Njs0) .AND. ijl(jfl) <= mjg(Nje0) ) THEN
iiloc(jfl) = iil(jfl) - mig(1) + 1
ijloc(jfl) = ijl(jfl) - mjg(1) + 1
IF( iil(jfl) >= mig(Nis0,nn_hls) .AND. iil(jfl) <= mig(Nie0,nn_hls) .AND. &
ijl(jfl) >= mjg(Njs0,nn_hls) .AND. ijl(jfl) <= mjg(Nje0,nn_hls) ) THEN
iiloc(jfl) = iil(jfl) - mig(1,nn_hls) + 1
ijloc(jfl) = ijl(jfl) - mjg(1,nn_hls) + 1
# else
iiloc(jfl) = iil(jfl)
ijloc(jfl) = ijl(jfl)
......
src/OCE/FLO/flodom.F90
View file @ 4dbfc9bb
......@@ -234,8 +234,8 @@ CONTAINS
zdyad = flo_dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) )
! Translation of this distances (in meter) in indexes
zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-1)
zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-1)
zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1,nn_hls)-1)
zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1,nn_hls)-1)
zgkfl(jfl) = (( gdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1,Kmm) - flzz(jfl) )* ikmfl(jfl)) &
& / ( gdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1,Kmm) &
& - gdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ,Kmm) ) &
......
src/OCE/FLO/florst.F90
View file @ 4dbfc9bb
......@@ -97,10 +97,10 @@ CONTAINS
!
IF( lk_mpp ) THEN
DO jfl = 1, jpnfl
IF( (INT(tpifl(jfl)) >= mig(Nis0)) .AND. &
&(INT(tpifl(jfl)) <= mig(Nie0)) .AND. &
&(INT(tpjfl(jfl)) >= mjg(Njs0)) .AND. &
&(INT(tpjfl(jfl)) <= mjg(Nje0)) ) THEN
IF( (INT(tpifl(jfl)) >= mig(Nis0,nn_hls)) .AND. &
&(INT(tpifl(jfl)) <= mig(Nie0,nn_hls)) .AND. &
&(INT(tpjfl(jfl)) >= mjg(Njs0,nn_hls)) .AND. &
&(INT(tpjfl(jfl)) <= mjg(Nje0,nn_hls)) ) THEN
iperproc(narea) = iperproc(narea)+1
ENDIF
END DO
......
src/OCE/FLO/flowri.F90
View file @ 4dbfc9bb
......@@ -103,8 +103,8 @@ CONTAINS
IF( lk_mpp ) THEN
iafloc = mi1( iafl )
ibfloc = mj1( ibfl )
iafloc = mi1( iafl, nn_hls )
ibfloc = mj1( ibfl, nn_hls )
IF( Nis0 <= iafloc .AND. iafloc <= Nie0 .AND. &
& Njs0 <= ibfloc .AND. ibfloc <= Nje0 ) THEN
......
src/OCE/ICB/icbclv.F90
View file @ 4dbfc9bb
......@@ -132,8 +132,8 @@ CONTAINS
!
newpt%lon = glamt(ji,jj) ! at t-point (centre of the cell)
newpt%lat = gphit(ji,jj)
newpt%xi = REAL( mig(ji), wp ) - ( nn_hls - 1 )
newpt%yj = REAL( mjg(jj), wp ) - ( nn_hls - 1 )
newpt%xi = REAL( mig(ji,nn_hls), wp ) - ( nn_hls - 1 )
newpt%yj = REAL( mjg(jj,nn_hls), wp ) - ( nn_hls - 1 )
!
newpt%uvel = 0._wp ! initially at rest
newpt%vvel = 0._wp
......
src/OCE/ICB/icbdyn.F90
View file @ 4dbfc9bb
......@@ -197,10 +197,10 @@ CONTAINS
IF( ii == ii0 .AND. ij == ij0 ) RETURN ! berg remains in the same cell
!
! map into current processor
ii0 = mi1( ii0 )
ij0 = mj1( ij0 )
ii = mi1( ii )
ij = mj1( ij )
ii0 = mi1( ii0, nn_hls )
ij0 = mj1( ij0, nn_hls )
ii = mi1( ii , nn_hls )
ij = mj1( ij , nn_hls )
!
! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
IF ( ln_M2016 .AND. ln_icb_grd ) THEN
......
src/OCE/ICB/icbini.F90
View file @ 4dbfc9bb
......@@ -140,7 +140,7 @@ CONTAINS
DO_2D( 1, 1, 1, 1 )
src_calving_hflx(ji,jj) = narea
src_calving (ji,jj) = nicbpack * mjg(jj) + mig(ji)
src_calving (ji,jj) = nicbpack * mjg(jj,nn_hls) + mig(ji,nn_hls)
END_2D
CALL lbc_lnk( 'icbini', src_calving_hflx, 'T', 1._wp )
CALL lbc_lnk( 'icbini', src_calving , 'T', 1._wp )
......@@ -156,7 +156,7 @@ CONTAINS
i2 = INT( i3/nicbpack )
i1 = i3 - i2*nicbpack
i3 = INT( src_calving_hflx(ji,jj) )
IF( i1 == mig(ji) .AND. i3 == narea ) THEN
IF( i1 == mig(ji,nn_hls) .AND. i3 == narea ) THEN
IF( nicbdi < 0 ) THEN ; nicbdi = ji
ELSE ; nicbei = ji
ENDIF
......@@ -172,7 +172,7 @@ CONTAINS
i2 = INT( i3/nicbpack )
i1 = i3 - i2*nicbpack
i3 = INT( src_calving_hflx(ji,jj) )
IF( i2 == mjg(jj) .AND. i3 == narea ) THEN
IF( i2 == mjg(jj,nn_hls) .AND. i3 == narea ) THEN
IF( nicbdj < 0 ) THEN ; nicbdj = jj
ELSE ; nicbej = jj
ENDIF
......@@ -361,8 +361,8 @@ CONTAINS
rn_test_box(1) < glamt(ji,jj) .AND. glamt(ji,jj) < rn_test_box(2) .AND. &
rn_test_box(3) < gphit(ji,jj) .AND. gphit(ji,jj) < rn_test_box(4) ) THEN
localberg%mass_scaling = rn_mass_scaling(iberg)
localpt%xi = REAL( mig(ji) - (nn_hls-1), wp )
localpt%yj = REAL( mjg(jj) - (nn_hls-1), wp )
localpt%xi = REAL( mig(ji,nn_hls) - (nn_hls-1), wp )
localpt%yj = REAL( mjg(jj,nn_hls) - (nn_hls-1), wp )
CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon )
localpt%mass = rn_initial_mass (iberg)
localpt%thickness = rn_initial_thickness(iberg)
......
src/OCE/ICB/icblbc.F90
View file @ 4dbfc9bb
......@@ -90,9 +90,9 @@ CONTAINS
this => first_berg
DO WHILE( ASSOCIATED(this) )
pt => this%current_point
IF( pt%xi > REAL(mig(nicbei),wp) + 0.5_wp ) THEN
IF( pt%xi > REAL(mig(nicbei,nn_hls),wp) + 0.5_wp ) THEN
pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp
ELSE IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp ) THEN
ELSE IF( pt%xi < REAL(mig(nicbdi,nn_hls),wp) - 0.5_wp ) THEN
pt%xi = ricb_left + MOD(pt%xi, 1._wp )
ENDIF
this => this%next
......@@ -125,10 +125,10 @@ CONTAINS
DO WHILE( ASSOCIATED(this) )
pt => this%current_point
ijne = INT( pt%yj + 0.5 )
IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
IF( pt%yj > REAL(mjg(nicbej,nn_hls),wp) + 0.5_wp ) THEN
!
iine = INT( pt%xi + 0.5 )
ipts = nicbfldpts (mi1(iine))
ipts = nicbfldpts (mi1(iine,nn_hls))
!
! moving across the cut line means both position and
! velocity must change
......@@ -228,7 +228,7 @@ CONTAINS
this => first_berg
DO WHILE (ASSOCIATED(this))
pt => this%current_point
IF( ipe_E >= 0 .AND. pt%xi > REAL(mig(nicbei),wp) + 0.5_wp - (nn_hls-1) ) THEN
IF( ipe_E >= 0 .AND. pt%xi > REAL(mig(nicbei,nn_hls),wp) + 0.5_wp - (nn_hls-1) ) THEN
tmpberg => this
this => this%next
ibergs_to_send_e = ibergs_to_send_e + 1
......@@ -241,7 +241,7 @@ CONTAINS
! now pack it into buffer and delete from list
CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e)
CALL icb_utl_delete(first_berg, tmpberg)
ELSE IF( ipe_W >= 0 .AND. pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp - (nn_hls-1) ) THEN
ELSE IF( ipe_W >= 0 .AND. pt%xi < REAL(mig(nicbdi,nn_hls),wp) - 0.5_wp - (nn_hls-1) ) THEN
tmpberg => this
this => this%next
ibergs_to_send_w = ibergs_to_send_w + 1
......@@ -320,7 +320,7 @@ CONTAINS
this => first_berg
DO WHILE (ASSOCIATED(this))
pt => this%current_point
IF( ipe_N >= 0 .AND. pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN
IF( ipe_N >= 0 .AND. pt%yj > REAL(mjg(nicbej,nn_hls),wp) + 0.5_wp - (nn_hls-1) ) THEN
tmpberg => this
this => this%next
ibergs_to_send_n = ibergs_to_send_n + 1
......@@ -330,7 +330,7 @@ CONTAINS
ENDIF
CALL icb_pack_into_buffer( tmpberg, obuffer_n, ibergs_to_send_n)
CALL icb_utl_delete(first_berg, tmpberg)
ELSE IF( ipe_S >= 0 .AND. pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp - (nn_hls-1) ) THEN
ELSE IF( ipe_S >= 0 .AND. pt%yj < REAL(mjg(nicbdj,nn_hls),wp) - 0.5_wp - (nn_hls-1) ) THEN
tmpberg => this
this => this%next
ibergs_to_send_s = ibergs_to_send_s + 1
......@@ -441,10 +441,10 @@ CONTAINS
this => first_berg
DO WHILE (ASSOCIATED(this))
pt => this%current_point
IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp - (nn_hls-1) .OR. &
pt%xi > REAL(mig(nicbei),wp) + 0.5_wp - (nn_hls-1) .OR. &
pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp - (nn_hls-1) .OR. &
pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN
IF( pt%xi < REAL(mig(nicbdi,nn_hls),wp) - 0.5_wp - (nn_hls-1) .OR. &
pt%xi > REAL(mig(nicbei,nn_hls),wp) + 0.5_wp - (nn_hls-1) .OR. &
pt%yj < REAL(mjg(nicbdj,nn_hls),wp) - 0.5_wp - (nn_hls-1) .OR. &
pt%yj > REAL(mjg(nicbej,nn_hls),wp) + 0.5_wp - (nn_hls-1) ) THEN
i = i + 1
WRITE(numicb,*) 'berg lost in halo: ', this%number(:)
WRITE(numicb,*) ' ', nimpp, njmpp
......@@ -514,8 +514,8 @@ CONTAINS
DO WHILE (ASSOCIATED(this))
pt => this%current_point
iine = INT( pt%xi + 0.5 ) + (nn_hls-1)
iproc = nicbflddest(mi1(iine))
IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN
iproc = nicbflddest(mi1(iine,nn_hls))
IF( pt%yj > REAL(mjg(nicbej,nn_hls),wp) + 0.5_wp - (nn_hls-1) ) THEN
IF( iproc == ifldproc ) THEN
!
IF( iproc /= narea ) THEN
......@@ -593,9 +593,9 @@ CONTAINS
pt => this%current_point
iine = INT( pt%xi + 0.5 ) + (nn_hls-1)
ijne = INT( pt%yj + 0.5 ) + (nn_hls-1)
ipts = nicbfldpts (mi1(iine))
iproc = nicbflddest(mi1(iine))
IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN
ipts = nicbfldpts (mi1(iine,nn_hls))
iproc = nicbflddest(mi1(iine,nn_hls))
IF( pt%yj > REAL(mjg(nicbej,nn_hls),wp) + 0.5_wp - (nn_hls-1) ) THEN
IF( iproc == ifldproc ) THEN
!
! moving across the cut line means both position and
......
src/OCE/ICB/icbrst.F90
View file @ 4dbfc9bb
......@@ -90,8 +90,8 @@ CONTAINS
ii = INT( localpt%xi + 0.5 ) + ( nn_hls-1 )
ij = INT( localpt%yj + 0.5 ) + ( nn_hls-1 )
! Only proceed if this iceberg is on the local processor (excluding halos).
IF ( ii >= mig(Nis0) .AND. ii <= mig(Nie0) .AND. &
& ij >= mjg(Njs0) .AND. ij <= mjg(Nje0) ) THEN
IF ( ii >= mig(Nis0,nn_hls) .AND. ii <= mig(Nie0,nn_hls) .AND. &
& ij >= mjg(Njs0,nn_hls) .AND. ij <= mjg(Nje0,nn_hls) ) THEN
CALL iom_get( ncid, jpdom_unknown, 'number', zdata(:) , ktime=jn, kstart=(/1/), kcount=(/nkounts/) )
localberg%number(:) = INT(zdata(:))
......@@ -244,16 +244,16 @@ CONTAINS
! global attributes
IF( lk_mpp ) THEN
! Set domain parameters (assume jpdom_local_full)
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number_total' , jpnij )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number' , narea-1 )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/ 1 , 2 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_global' , (/ Ni0glo , Nj0glo /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_local' , (/ Ni_0 , Nj_0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_first' , (/ mig0(Nis0), mjg0(Njs0) /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_last' , (/ mig0(Nie0), mjg0(Nje0) /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_start', (/ 0 , 0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_end' , (/ 0 , 0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_type' , 'BOX' )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number_total' , jpnij )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_number' , narea-1 )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/ 1 , 2 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_global' , (/ Ni0glo , Nj0glo /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_size_local' , (/ Ni_0 , Nj_0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_first' , (/ mig(Nis0,0), mjg(Njs0,0) /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_position_last' , (/ mig(Nie0,0), mjg(Nje0,0) /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_start', (/ 0 , 0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_halo_size_end' , (/ 0 , 0 /) )
nret = NF90_PUT_ATT( ncid, NF90_GLOBAL, 'DOMAIN_type' , 'BOX' )
ENDIF
IF (associated(first_berg)) then
......
src/OCE/ICB/icbthm.F90
View file @ 4dbfc9bb
......@@ -31,6 +31,8 @@ MODULE icbthm
PUBLIC icb_thm ! routine called in icbstp.F90 module
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: icbthm.F90 15088 2021-07-06 13:03:34Z acc $
......@@ -112,9 +114,9 @@ CONTAINS
zxi = pt%xi ! position in (i,j) referential
zyj = pt%yj
ii = INT( zxi + 0.5 ) ! T-cell of the berg
ii = mi1( ii + (nn_hls-1) )
ii = mi1( ii + (nn_hls-1), nn_hls )
ij = INT( zyj + 0.5 )
ij = mj1( ij + (nn_hls-1) )
ij = mj1( ij + (nn_hls-1), nn_hls )
zVol = zT * zW * zL
! Environment
......@@ -287,8 +289,8 @@ CONTAINS
! now use melt and associated heat flux in ocean (or not)
!
IF(.NOT. ln_passive_mode ) THEN
emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:)
emp (A2D(0)) = emp (A2D(0)) - berg_grid%floating_melt(A2D(0))
qns (:,:) = qns (:,:) + berg_grid%calving_hflx (A2D(0))
ENDIF
!
END SUBROUTINE icb_thm
......
src/OCE/ICB/icbutl.F90
View file @ 4dbfc9bb
......@@ -57,6 +57,7 @@ MODULE icbutl
PUBLIC icb_utl_heat ! routine called in icbdia module
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -101,7 +102,7 @@ CONTAINS
CALL lbc_lnk_icb( 'icbutl', ua_e , 'U', -1._wp, 1, 1 )
CALL lbc_lnk_icb( 'icbutl', va_e , 'V', -1._wp, 1, 1 )
#if defined key_si3
hi_e(1:jpi, 1:jpj) = hm_i (:,:)
hi_e(A2D(0)) = hm_i (:,:) ! clem: something is wrong here (hm_i defined in the interior only) but I do not what to do
ui_e(1:jpi, 1:jpj) = u_ice(:,:)
vi_e(1:jpi, 1:jpj) = v_ice(:,:)
!
......@@ -312,18 +313,18 @@ CONTAINS
!
IF (TRIM(cd_type) == 'T' ) THEN
ierr = 0
IF ( kii < mig( 1 ) ) THEN ; ierr = ierr + 1
ELSEIF( kii >= mig(jpi) ) THEN ; ierr = ierr + 1
IF ( kii < mig( 1 ,nn_hls) ) THEN ; ierr = ierr + 1
ELSEIF( kii >= mig(jpi,nn_hls) ) THEN ; ierr = ierr + 1
ENDIF
!
IF ( kij < mjg( 1 ) ) THEN ; ierr = ierr + 1
ELSEIF( kij >= mjg(jpj) ) THEN ; ierr = ierr + 1
IF ( kij < mjg( 1 ,nn_hls) ) THEN ; ierr = ierr + 1
ELSEIF( kij >= mjg(jpj,nn_hls) ) THEN ; ierr = ierr + 1
ENDIF
!
IF ( ierr > 0 ) THEN
WRITE(numicb,*) 'bottom left corner T point out of bound'
WRITE(numicb,*) pi, kii, mig( 1 ), mig(jpi)
WRITE(numicb,*) pj, kij, mjg( 1 ), mjg(jpj)
WRITE(numicb,*) pi, kii, mig( 1,nn_hls ), mig(jpi,nn_hls)
WRITE(numicb,*) pj, kij, mjg( 1,nn_hls ), mjg(jpj,nn_hls)
WRITE(numicb,*) pmsk
CALL FLUSH(numicb)
CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error).' , &
......@@ -335,13 +336,13 @@ CONTAINS
! find position in this processor. Prevent near edge problems (see #1389)
! (PM) will be useless if extra halo is used in NEMO
!
IF ( kii <= mig(1)-1 ) THEN ; kii = 0
ELSEIF( kii > mig(jpi) ) THEN ; kii = jpi
ELSE ; kii = mi1(kii)
IF ( kii <= mig(1,nn_hls)-1 ) THEN ; kii = 0
ELSEIF( kii > mig(jpi,nn_hls) ) THEN ; kii = jpi
ELSE ; kii = mi1(kii,nn_hls)
ENDIF
IF ( kij <= mjg(1)-1 ) THEN ; kij = 0
ELSEIF( kij > mjg(jpj) ) THEN ; kij = jpj
ELSE ; kij = mj1(kij)
IF ( kij <= mjg(1,nn_hls)-1 ) THEN ; kij = 0
ELSEIF( kij > mjg(jpj,nn_hls) ) THEN ; kij = jpj
ELSE ; kij = mj1(kij,nn_hls)
ENDIF
!
! define mask array
......@@ -462,8 +463,8 @@ CONTAINS
zj = pj - REAL(ij,wp)
! conversion to local domain (no need to do a sanity check already done in icbpos)
ii = mi1(ii) + (nn_hls-1)
ij = mj1(ij) + (nn_hls-1)
ii = mi1(ii,nn_hls) + (nn_hls-1)
ij = mj1(ij,nn_hls) + (nn_hls-1)
!
IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN
IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NE quadrant
......
src/OCE/IOM/iom.F90
View file @ 4dbfc9bb
......@@ -1202,6 +1202,7 @@ CONTAINS
CHARACTER(LEN=1) :: cl_type ! local value of cd_type
LOGICAL :: ll_only3rd ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
INTEGER :: inlev ! number of levels for 3D data
INTEGER :: ihls ! local value of the halo size
REAL(dp) :: gma, gmi
!---------------------------------------------------------------------
CHARACTER(LEN=lc) :: context
......@@ -1298,7 +1299,7 @@ CONTAINS
ENDIF
ELSE ! not a 1D array as pv_r1d requires jpdom_unknown
! we do not read the overlap and the extra-halos -> from Nis0 to Nie0 and from Njs0 to Nje0
IF( idom == jpdom_global ) istart(1:2) = (/ mig0(Nis0), mjg0(Njs0) /)
IF( idom == jpdom_global ) istart(1:2) = (/ mig(Nis0,0), mjg(Njs0,0) /)
icnt(1:2) = (/ Ni_0, Nj_0 /)
IF( PRESENT(pv_r3d) ) THEN
IF( idom == jpdom_auto_xy ) THEN
......@@ -1327,14 +1328,26 @@ CONTAINS
IF( irankpv == 1 ) ishape(1:1) = SHAPE(pv_r1d)
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2) ! index of the array to be read
ctmp1 = 'd'
ELSE
IF( irankpv == 2 ) THEN
ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0 )) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
ENDIF
IF( irankpv == 3 ) THEN
ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:)) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN ! array with 0 halo
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0 ! index of the array to be read
ctmp1 = 'd(:,:'
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN ! array with nn_hls halos
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0 ! index of the array to be read
ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0'
ELSEIF( ishape(1) == Ni_0+2 .AND. ishape(2) == Nj_0+2 ) THEN ! nn_hls = 2 and array with 1 halo
ix1 = 2 ; ix2 = Ni_0+1 ; iy1 = 2 ; iy2 = Nj_0+1 ! index of the array to be read
ctmp1 = 'd(2:Ni_0+1,2:Ni_0+1'
ELSE
CALL ctl_stop( 'iom_get_123d: should have been an impossible case...' )
ENDIF
ishape(1:2) = (/ Ni_0, Nj_0 /) ! update and force ishape to match the inner domain
IF( irankpv == 3 ) ctmp1 = TRIM(ctmp1)//',:'
ctmp1 = TRIM(ctmp1)//')'
ENDIF
DO jl = 1, irankpv
WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
......@@ -1347,11 +1360,6 @@ CONTAINS
!-
IF( idvar > 0 .AND. istop == nstop ) THEN ! no additional errors until this point...
!
! find the right index of the array to be read
IF( idom /= jpdom_unknown ) THEN ; ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSE ; ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2)
ENDIF
CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
IF( istop == nstop ) THEN ! no additional errors until this point...
......@@ -1359,13 +1367,15 @@ CONTAINS
cl_type = 'T'
IF( PRESENT(cd_type) ) cl_type = cd_type
zsgn = 1._wp
IF( PRESENT(psgn ) ) zsgn = psgn
!--- overlap areas and extra hallows (mpp)
IF( PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
!--- halos and NP folding (NP folding to be done even if no halos)
IF( idom /= jpdom_unknown .AND. cl_type /= 'Z' .AND. ( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) ) THEN
zsgn = 1._wp
IF( PRESENT(psgn ) ) zsgn = psgn
IF( PRESENT(pv_r2d) ) THEN
CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
ELSEIF( PRESENT(pv_r3d) ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
ENDIF
ENDIF
!
ELSE
......@@ -1391,13 +1401,13 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'XIOS RST READ (3D): ',TRIM(cdvar)
CALL xios_recv_field( trim(cdvar), pv_r3d(:, :, :))
IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill)
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
ENDIF
ELSEIF( PRESENT(pv_r2d) ) THEN
IF(lwp) WRITE(numout,*) 'XIOS RST READ (2D): ', TRIM(cdvar)
CALL xios_recv_field( trim(cdvar), pv_r2d(:, :))
IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk('iom', pv_r2d, cl_type, zsgn, kfillmode = kfill)
CALL lbc_lnk('iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
ENDIF
ELSEIF( PRESENT(pv_r1d) ) THEN
IF(lwp) WRITE(numout,*) 'XIOS RST READ (1D): ', TRIM(cdvar)
......@@ -2322,11 +2332,11 @@ CONTAINS
LOGICAL, INTENT(IN) :: ldxios, ldrxios
!!----------------------------------------------------------------------
!
CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0)
CALL iom_set_domain_attr("grid_"//cdgrd,ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig(Nis0,0)-1,jbegin=mjg(Njs0,0)-1,ni=Ni_0,nj=Nj_0)
CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni=jpi, data_jbegin = -nn_hls, data_nj=jpj)
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ni_glo = Ni0glo, nj_glo = Nj0glo, &
& ibegin = mig0(Nis0) - 1, jbegin = mjg0(Njs0) - 1, ni = Ni_0, nj = Nj_0)
& ibegin = mig(Nis0,0) - 1, jbegin = mjg(Njs0,0) - 1, ni = Ni_0, nj = Nj_0)
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", data_dim=2, data_ibegin = 0, data_ni=Ni_0, data_jbegin = 0, data_nj=Nj_0)
IF( ln_tile ) THEN
......@@ -2336,14 +2346,14 @@ CONTAINS
idb(jn) = -nn_hls ! Tile data offset (halo size)
END DO
! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
idb(:) = 0
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
......@@ -2406,7 +2416,7 @@ CONTAINS
END SELECT
!
z_fld(:,:) = 1._wp
CALL lbc_lnk( 'iom', z_fld, cdgrd, -1.0_wp ) ! Working array for location of northfold
CALL lbc_lnk( 'iom', z_fld, cdgrd, -1.0_wp, ldfull = .TRUE. ) ! Working array for location of northfold
!
! Cell vertices that can be defined
DO_2D( 0, 0, 0, 0 )
......@@ -2452,8 +2462,8 @@ CONTAINS
!
! CALL dom_ngb( -168.53_wp, 65.03_wp, ix, iy, 'T' ) ! i-line that passes through Bering Strait: Reference latitude (used in plots)
CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) ! i-line that passes near the North Pole : Reference latitude (used in plots)
CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig(Nis0,0)-1, jbegin=mjg(Njs0,0)-1, ni=Ni_0, nj=Nj_0)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin=0, data_ni=Ni_0, data_jbegin=0, data_nj=Nj_0)
CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp), &
& latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))
CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
......
src/OCE/IOM/iom_nf90.F90
View file @ 4dbfc9bb
......@@ -40,6 +40,8 @@ MODULE iom_nf90
MODULE PROCEDURE iom_nf90_rp0123d_dp
END INTERFACE
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: iom_nf90.F90 14433 2021-02-11 08:06:49Z smasson $
......@@ -144,16 +146,16 @@ CONTAINS
END SELECT
CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo)
! global attributes
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_number_total' , jpnij ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_number' , narea-1 ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/ 1 , 2 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_size_global' , (/ Ni0glo , Nj0glo /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_size_local' , (/ Ni_0 , Nj_0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_position_first' , (/ mig0(Nis0), mjg0(Njs0) /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_position_last' , (/ mig0(Nie0), mjg0(Nje0) /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_halo_size_start', (/ 0 , 0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_halo_size_end' , (/ 0 , 0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_type' , 'BOX' ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_number_total' , jpnij ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_number' , narea-1 ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_dimensions_ids' , (/ 1 , 2 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_size_global' , (/ Ni0glo , Nj0glo /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_size_local' , (/ Ni_0 , Nj_0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_position_first' , (/ mig(Nis0,0), mjg(Njs0,0) /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_position_last' , (/ mig(Nie0,0), mjg(Nje0,0) /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_halo_size_start', (/ 0 , 0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_halo_size_end' , (/ 0 , 0 /) ), clinfo)
CALL iom_nf90_check(NF90_PUT_ATT(if90id,NF90_GLOBAL, 'DOMAIN_type' , 'BOX' ), clinfo)
ELSE !* the file should be open for read mode so it must exist...
CALL ctl_stop( TRIM(clinfo), ' should be impossible case...' )
ENDIF
......@@ -544,7 +546,7 @@ CONTAINS
INTEGER :: idvar ! variable id
INTEGER :: jd ! dimension loop counter
INTEGER :: ix1, ix2, iy1, iy2 ! subdomain indexes
INTEGER, DIMENSION(4) :: idimsz ! dimensions size
INTEGER, DIMENSION(3) :: ishape ! dimensions size
INTEGER, DIMENSION(4) :: idimid ! dimensions id
CHARACTER(LEN=256) :: clinfo ! info character
INTEGER :: if90id ! nf90 file identifier
......@@ -627,11 +629,9 @@ CONTAINS
itype = NF90_DOUBLE
ENDIF
IF( PRESENT(pv_r0d) ) THEN
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, iom_file(kiomid)%nvid(idvar) ), clinfo )
ELSE
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
ENDIF
lchunk = .false.
IF( snc4set%luse .AND. idims == 4 ) lchunk = .true.
......@@ -673,23 +673,13 @@ CONTAINS
ENDIF
! on what kind of domain must the data be written?
IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
IF( idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
! write dimension variables if it is not already done
! =============
! trick: is defined to 0 => dimension variable are defined but not yet written
IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN ! time_counter = 0
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(A2D(0)) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(A2D(0)) ), clinfo )
SELECT CASE (iom_file(kiomid)%comp)
CASE ('OCE')
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, gdept_1d ), clinfo )
......@@ -704,6 +694,19 @@ CONTAINS
iom_file(kiomid)%dimsz(1, 4) = 1 ! so we don't enter this IF case any more...
IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
ENDIF
IF( PRESENT(pv_r2d) ) ishape(1:2) = SHAPE(pv_r2d)
IF( PRESENT(pv_r3d) ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN ! array with 0 halo
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN ! array with nn_hls halos
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSEIF( ishape(1) == Ni_0+2 .AND. ishape(2) == Nj_0+2 ) THEN ! nn_hls = 2 and array with 1 halo
ix1 = 2 ; ix2 = Ni_0+1 ; iy1 = 2 ; iy2 = Nj_0+1
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
ENDIF
! write the data
......@@ -712,7 +715,7 @@ CONTAINS
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d ), clinfo )
ELSEIF( PRESENT(pv_r1d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:) ), clinfo )
ELSEIF( PRESENT(pv_r2d) ) THEN
ELSEIF( PRESENT(pv_r2d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2) ), clinfo )
ELSEIF( PRESENT(pv_r3d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
......
src/OCE/IOM/prtctl.F90
View file @ 4dbfc9bb
......@@ -119,11 +119,11 @@ CONTAINS
!! clinfo3 : additional information
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: ktab2d_1, ktab3d_1, ktab4d_1, ktab2d_2, ktab3d_2
REAL(2*wp), DIMENSION(A2D_T(ktab2d_1)) , INTENT(in), OPTIONAL :: tab2d_1
REAL(2*wp), DIMENSION(A2D_T(ktab3d_1),:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(2*wp), DIMENSION(A2D_T(ktab4d_1),:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(2*wp), DIMENSION(A2D_T(ktab2d_2)) , INTENT(in), OPTIONAL :: tab2d_2
REAL(2*wp), DIMENSION(A2D_T(ktab3d_2),:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(2*wp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_1
REAL(2*wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(2*wp), DIMENSION(:,:,:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(2*wp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_2
REAL(2*wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask1
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask2
CHARACTER(len=*), DIMENSION(:) , INTENT(in), OPTIONAL :: clinfo ! information about the tab3d array
......@@ -135,11 +135,20 @@ CONTAINS
CHARACTER(len=30) :: cl1, cl2
CHARACTER(len=6) :: clfmt
INTEGER :: jn, jl, kdir
INTEGER :: iis, iie, jjs, jje
INTEGER :: ipi1, ipi2, ipim1, ipim2
INTEGER :: isht1, isht2, ishtm1, ishtm2
INTEGER :: ii1s, ii1e, jj1s, jj1e
INTEGER :: ii2s, ii2e, jj2s, jj2e
INTEGER :: iim1s, iim1e, jjm1s, jjm1e
INTEGER :: iim2s, iim2e, jjm2s, jjm2e
INTEGER :: itra, inum
REAL(2*wp) :: zsum1, zsum2, zvctl1, zvctl2
!!----------------------------------------------------------------------
!
IF( ( ktab2d_1 * ktab3d_1 * ktab4d_1 * ktab2d_2 * ktab3d_2 ) /= 0 ) THEN
CALL ctl_stop( 'prt_ctl is a debugging toll that should not be used with tiles' )
ENDIF
! Arrays, scalars initialization
cl1 = ''
cl2 = ''
......@@ -154,22 +163,55 @@ CONTAINS
IF( wp == sp ) clfmt = 'D23.16' ! 16 significant numbers
IF( wp == dp ) clfmt = 'D41.34' ! 34 significant numbers
IF( PRESENT(tab2d_1) ) ipi1 = SIZE(tab2d_1,1)
IF( PRESENT(tab3d_1) ) ipi1 = SIZE(tab3d_1,1)
IF( PRESENT(tab4d_1) ) ipi1 = SIZE(tab4d_1,1)
isht1 = ( jpi - ipi1 ) / 2
ipi2 = -1 ! default definition
IF( PRESENT(tab2d_2) ) ipi2 = SIZE(tab2d_2,1)
IF( PRESENT(tab3d_2) ) ipi2 = SIZE(tab3d_2,1)
isht2 = ( jpi - ipi2 ) / 2
ipim1 = -1 ! default definition
IF( PRESENT(mask1) ) ipim1 = SIZE(mask1,1)
ishtm1 = ( jpi - ipim1 ) / 2
ipim2 = -1 ! default definition
IF( PRESENT(mask2) ) ipim2 = SIZE(mask2,1)
ishtm2 = ( jpi - ipim2 ) / 2
! Loop over each sub-domain, i.e. the total number of processors ijsplt
DO jl = 1, SIZE(nall_ictls)
ii1s = MAX( nall_ictls(jl), Nis0 ) - isht1
ii1e = MIN( nall_ictle(jl), Nie0 ) - isht1
jj1s = MAX( nall_jctls(jl), Njs0 ) - isht1
jj1e = MIN( nall_jctle(jl), Nje0 ) - isht1
! define shoter names...
iis = MAX( nall_ictls(jl), ntsi )
iie = MIN( nall_ictle(jl), ntei )
jjs = MAX( nall_jctls(jl), ntsj )
jje = MIN( nall_jctle(jl), ntej )
ii2s = MAX( nall_ictls(jl), Nis0 ) - isht2
ii2e = MIN( nall_ictle(jl), Nie0 ) - isht2
jj2s = MAX( nall_jctls(jl), Njs0 ) - isht2
jj2e = MIN( nall_jctle(jl), Nje0 ) - isht2
iim1s = MAX( nall_ictls(jl), Nis0 ) - ishtm1
iim1e = MIN( nall_ictle(jl), Nie0 ) - ishtm1
jjm1s = MAX( nall_jctls(jl), Njs0 ) - ishtm1
jjm1e = MIN( nall_jctle(jl), Nje0 ) - ishtm1
iim2s = MAX( nall_ictls(jl), Nis0 ) - ishtm2
iim2e = MIN( nall_ictle(jl), Nie0 ) - ishtm2
jjm2s = MAX( nall_jctls(jl), Njs0 ) - ishtm2
jjm2e = MIN( nall_jctle(jl), Nje0 ) - ishtm2
IF( PRESENT(clinfo) ) THEN ; inum = numprt_top(jl)
ELSE ; inum = numprt_oce(jl)
ENDIF
! Compute the sum control only where the tile domain and control print area overlap
IF( iie >= iis .AND. jje >= jjs ) THEN
IF( ii1e >= ii1s .AND. jj1e >= jj1s ) THEN
DO jn = 1, itra
IF( PRESENT(clinfo3) ) THEN
......@@ -188,32 +230,42 @@ CONTAINS
! 2D arrays
IF( PRESENT(tab2d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(iis:iie,jjs:jje,1) )
ELSE ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) )
IF( PRESENT(mask1) ) THEN
zsum1 = SUM( tab2d_1(ii1s:ii1e,jj1s:jj1e) * mask1(iim1s:iim1e,jjm1s:jjm1e,1) )
ELSE
zsum1 = SUM( tab2d_1(ii1s:ii1e,jj1s:jj1e) )
ENDIF
ENDIF
IF( PRESENT(tab2d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(iis:iie,jjs:jje,1) )
ELSE ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) )
IF( PRESENT(mask2) ) THEN
zsum2 = SUM( tab2d_2(ii2s:ii2e,jj2s:jj2e) * mask2(iim2s:iim2e,jjm2s:jjm2e,1) )
ELSE
zsum2 = SUM( tab2d_2(ii2s:ii2e,jj2s:jj2e) )
ENDIF
ENDIF
! 3D arrays
IF( PRESENT(tab3d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask1) ) THEN
zsum1 = SUM( tab3d_1(ii1s:ii1e,jj1s:jj1e,1:kdir) * mask1(iim1s:iim1e,jjm1s:jjm1e,1:kdir) )
ELSE
zsum1 = SUM( tab3d_1(ii1s:ii1e,jj1s:jj1e,1:kdir) )
ENDIF
ENDIF
IF( PRESENT(tab3d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask2) ) THEN
zsum2 = SUM( tab3d_2(ii2s:ii2e,jj2s:jj2e,1:kdir) * mask2(iim2s:iim2e,jjm2s:jjm2e,1:kdir) )
ELSE
zsum2 = SUM( tab3d_2(ii2s:ii2e,jj2s:jj2e,1:kdir) )
ENDIF
ENDIF
! 4D arrays
IF( PRESENT(tab4d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) )
IF( PRESENT(mask1) ) THEN
zsum1 = SUM( tab4d_1(ii1s:ii1e,jj1s:jj1e,1:kdir,jn) * mask1(iim1s:iim1e,jjm1s:jjm1e,1:kdir) )
ELSE
zsum1 = SUM( tab4d_1(ii1s:ii1e,jj1s:jj1e,1:kdir,jn) )
ENDIF
ENDIF
......@@ -460,22 +512,22 @@ CONTAINS
!
idg = MAXVAL( (/ nall_ictls(jl), nall_ictle(jl), nall_jctls(jl), nall_jctle(jl) /) ) ! temporary use of idg
idg = INT(LOG10(REAL(idg,wp))) + 1 ! how many digits do we use?
idg2 = MAXVAL( (/ mig0(nall_ictls(jl)), mig0(nall_ictle(jl)), mjg0(nall_jctls(jl)), mjg0(nall_jctle(jl)) /) )
idg2 = MAXVAL( (/ mig(nall_ictls(jl),0), mig(nall_ictle(jl),0), mjg(nall_jctls(jl),0), mjg(nall_jctle(jl),0) /) )
idg2 = INT(LOG10(REAL(idg2,wp))) + 1 ! how many digits do we use?
WRITE(clfmt2, "('(18x, 13a1, a9, i', i1, ', a2, i',i1,', a2, 13a1)')") idg, idg2
WRITE(clfmt3, "('(18x, a1, ', i2,'x, a1)')") 13+9+idg+2+idg2+2+13 - 2
WRITE(clfmt4, "('(', i2,'x, a9, i', i1,', a2, i', i1,', a2, ', i2,'x, a9, i', i1,', a2, i', i1,', a2)')") &
& 18-7, idg, idg2, 13+9+idg+2+idg2+2+13 - (2+idg+2+idg2+2+8), idg, idg2
WRITE(inum,clfmt2) ('-', ji=1,13), ' jctle = ', nall_jctle(jl), ' (', mjg0(nall_jctle(jl)), ') ', ('-', ji=1,13)
WRITE(inum,clfmt2) ('-', ji=1,13), ' jctle = ', nall_jctle(jl), ' (', mjg(nall_jctle(jl),0), ') ', ('-', ji=1,13)
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt4) ' ictls = ', nall_ictls(jl), ' (', mig0(nall_ictls(jl)), ') ', &
& ' ictle = ', nall_ictle(jl), ' (', mig0(nall_ictle(jl)), ') '
WRITE(inum,clfmt4) ' ictls = ', nall_ictls(jl), ' (', mig(nall_ictls(jl),0), ') ', &
& ' ictle = ', nall_ictle(jl), ' (', mig(nall_ictle(jl),0), ') '
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt3) '|', '|'
WRITE(inum,clfmt2) ('-', ji=1,13), ' jctls = ', nall_jctls(jl), ' (', mjg0(nall_jctls(jl)), ') ', ('-', ji=1,13)
WRITE(inum,clfmt2) ('-', ji=1,13), ' jctls = ', nall_jctls(jl), ' (', mjg(nall_jctls(jl),0), ') ', ('-', ji=1,13)
WRITE(inum,*)
WRITE(inum,*)
!
......
src/OCE/ISF/isfcpl.F90
View file @ 4dbfc9bb
......@@ -736,7 +736,8 @@ CONTAINS
END IF
!
! update isfpts structure
sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
sisfpts(kpts) = isfcons(mig(ki,nn_hls), mjg(kj,nn_hls), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, &
& glamt(ki,kj), gphit(ki,kj), ifind )
!
END SUBROUTINE update_isfpts
!
......@@ -761,8 +762,8 @@ CONTAINS
IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
!
! fill the correction array
DO jj = mj0(ijg),mj1(ijg)
DO ji = mi0(iig),mi1(iig)
DO jj = mj0(ijg,nn_hls),mj1(ijg,nn_hls)
DO ji = mi0(iig,nn_hls),mi1(iig,nn_hls)
! correct the vol_flx and corresponding heat/salt flx in the closest cell
risfcpl_cons_vol(ji,jj,kk) = risfcpl_cons_vol(ji,jj,kk ) + pvolinc
risfcpl_cons_tsc(ji,jj,kk,jp_sal) = risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
......
src/OCE/LBC/lbc_lnk_call_generic.h90
View file @ 4dbfc9bb
......@@ -27,7 +27,7 @@
& , pt21, cdna21, psgn21, pt22, cdna22, psgn22, pt23, cdna23, psgn23, pt24, cdna24, psgn24 &
& , pt25, cdna25, psgn25, pt26, cdna26, psgn26, pt27, cdna27, psgn27, pt28, cdna28, psgn28 &
& , pt29, cdna29, psgn29, pt30, cdna30, psgn30 &
& , kfillmode, pfillval, khls, lsend, lrecv, ld4only )
& , kfillmode, pfillval, lsend, lrecv, ld4only, ldfull )
!!---------------------------------------------------------------------
CHARACTER(len=*) , INTENT(in ) :: cdname ! name of the calling subroutine
REAL(PRECISION), DIMENSION(DIMS) , TARGET, CONTIGUOUS, INTENT(inout) :: pt1 ! arrays on which the lbc is applied
......@@ -50,9 +50,9 @@
& psgn30
INTEGER , OPTIONAL , INTENT(in ) :: kfillmode ! filling method for halo over land (default = constant)
REAL(PRECISION) , OPTIONAL , INTENT(in ) :: pfillval ! background value (used at closed boundaries)
INTEGER , OPTIONAL , INTENT(in ) :: khls ! halo size, default = nn_hls
LOGICAL, DIMENSION(8), OPTIONAL , INTENT(in ) :: lsend, lrecv ! indicate how communications are to be carried out
LOGICAL , OPTIONAL , INTENT(in ) :: ld4only ! if .T., do only 4-neighbour comm (ignore corners)
LOGICAL , OPTIONAL , INTENT(in ) :: ldfull ! .true. if we also update the last line of the inner domain
!!
INTEGER :: kfld ! number of elements that will be attributed
TYPE(PTR_4d_/**/PRECISION), DIMENSION(30) :: ptab_ptr ! pointer array
......@@ -96,15 +96,11 @@
IF( PRESENT(psgn29) ) CALL load_ptr_/**/XD/**/_/**/PRECISION( pt29, cdna29, psgn29, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
IF( PRESENT(psgn30) ) CALL load_ptr_/**/XD/**/_/**/PRECISION( pt30, cdna30, psgn30, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
!
#if ! defined key_mpi2
IF( nn_comm == 1 ) THEN
CALL lbc_lnk_pt2pt( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, khls, lsend, lrecv, ld4only )
CALL lbc_lnk_pt2pt( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ld4only, ldfull )
ELSE
CALL lbc_lnk_neicoll( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, khls, lsend, lrecv, ld4only )
CALL lbc_lnk_neicoll( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ld4only, ldfull )
ENDIF
#else
CALL lbc_lnk_pt2pt( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, khls, lsend, lrecv, ld4only )
#endif
!
END SUBROUTINE lbc_lnk_call_/**/XD/**/_/**/PRECISION
......