MODULE dynnxt !!========================================================================= !! *** MODULE dynnxt *** !! Ocean dynamics: time stepping !!========================================================================= !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code !! ! 1990-10 (C. Levy, G. Madec) !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions !! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0 !! 8.2 ! 1997-04 (A. Weaver) Euler forward step !! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes !! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification !! 4.0 ! 2018-07 (G. Madec) 1- z-star (s-star) only !! 2- remove dom_vvl_interpol !!------------------------------------------------------------------------- !!------------------------------------------------------------------------- !! dyn_nxt : obtain the next (after) horizontal velocity !!------------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE sbc_oce ! Surface boundary condition: ocean fields USE sbcrnf ! river runoffs USE sbcisf ! ice shelf USE phycst ! physical constants USE dynadv ! dynamics: vector invariant versus flux form USE dynspg_ts ! surface pressure gradient: split-explicit scheme USE domvvl ! variable volume USE bdy_oce , ONLY : ln_bdy USE bdydta ! ocean open boundary conditions USE bdydyn ! ocean open boundary conditions USE bdyvol ! ocean open boundary condition (bdy_vol routines) USE trd_oce ! trends: ocean variables USE trddyn ! trend manager: dynamics USE trdken ! trend manager: kinetic energy ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE lbclnk ! lateral boundary condition (or mpp link) USE lib_mpp ! MPP library USE prtctl ! Print control USE timing ! Timing #if defined key_agrif USE agrif_oce_interp #endif IMPLICIT NONE PRIVATE PUBLIC dyn_nxt ! routine called by step.F90 !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL licence (./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_nxt( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_nxt *** !! !! ** Purpose : Finalize after horizontal velocity. Apply the boundary !! condition on the after velocity, achieve the time stepping !! by applying the Asselin filter on now fields and swapping !! the fields. !! !! ** Method : * Ensure after velocities transport matches time splitting !! estimate (ln_dynspg_ts=T) !! !! * Apply lateral boundary conditions on after velocity !! at the local domain boundaries through lbc_lnk call, !! at the one-way open boundaries (ln_bdy=T), !! at the AGRIF zoom boundaries (lk_agrif=T) !! !! * Apply the time filter applied and swap of the dynamics !! arrays to start the next time step: !! (ub,vb) = (un,vn) + rn_atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] !! (un,vn) = (ua,va). !! Note that with flux form advection and non linear free surface, !! the time filter is applied on thickness weighted velocity. !! As a result, dyn_nxt MUST be called after tra_nxt. !! !! ** Action : ub,vb filtered before horizontal velocity of next time-step !! un,vn now horizontal velocity of next time-step !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ikt ! local integers REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva REAL(wp), DIMENSION(jpi,jpj) :: z_ssh_h0, zsshu, zsshv !!---------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('dyn_nxt') IF( ln_dynspg_ts ) ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) ) IF( l_trddyn ) ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) ) ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~' ENDIF IF ( ln_dynspg_ts ) THEN ! Ensure below that barotropic velocities match time splitting estimate ! Compute actual transport and replace it with ts estimate at "after" time step zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) DO jk = 2, jpkm1 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) END DO DO jk = 1, jpkm1 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) END DO ! IF( .NOT.ln_bt_fw ) THEN ! Remove advective velocity from "now velocities" ! prior to asselin filtering ! In the forward case, this is done below after asselin filtering ! so that asselin contribution is removed at the same time DO jk = 1, jpkm1 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) ) * umask(:,:,jk) vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) ) * vmask(:,:,jk) END DO ENDIF ENDIF ! Update after velocity on domain lateral boundaries ! -------------------------------------------------- # if defined key_agrif CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries # endif ! CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. ) !* local domain boundaries ! ! !* BDY open boundaries IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) !!$ Do we need a call to bdy_vol here?? ! IF( l_trddyn ) THEN !* prepare the atf trend computation + some diagnostics IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) ! Kinetic energy and Conversion IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends IF( ln_dynadv_vec ) THEN zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt ELSE zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt ENDIF CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin filter CALL iom_put( "vtrd_tot", zva ) ENDIF zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter zva(:,:,:) = vn(:,:,:) ! (caution: the Asselin filter trends computation will be shifted by 1 timestep) ENDIF ! Time filter and swap of dynamics arrays ! --------------------------------------- IF( l_1st_euler ) THEN !== Euler at 1st time-step ==! (swap only) DO jk = 1, jpkm1 un(:,:,jk) = ua(:,:,jk) ! un <-- ua vn(:,:,jk) = va(:,:,jk) END DO IF( .NOT.ln_linssh ) THEN ! e3._n <-- e3._a DO jk = 1, jpkm1 e3t_n(:,:,jk) = e3t_a(:,:,jk) e3u_n(:,:,jk) = e3u_a(:,:,jk) e3v_n(:,:,jk) = e3v_a(:,:,jk) END DO ENDIF ! ELSE !== Leap-Frog ==! (Asselin filter and swap) ! ! ! =============! IF( ln_linssh ) THEN ! Fixed volume ! ! ! =============! DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) ! ub(ji,jj,jk) = zuf ! ub <-- filtered velocity vb(ji,jj,jk) = zvf un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua vn(ji,jj,jk) = va(ji,jj,jk) END DO END DO END DO ! ! ================! ELSE ! Variable volume ! ! ! ================! ! Before scale factor at t-points (used as a now filtered scale factor until the swap) DO jk = 1, jpkm1 e3t_b(:,:,jk) = e3t_n(:,:,jk) + rn_atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) END DO ! Add volume filter correction: compatibility with tracer advection scheme ! => time filter + conservation correction (only at the first level) zcoef = rn_atfp * rn_Dt * r1_rho0 ! e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) ! IF ( ln_rnf ) THEN IF( ln_rnf_depth ) THEN DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too DO jj = 1, jpj DO ji = 1, jpi IF( jk <= nk_rnf(ji,jj) ) THEN e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) & & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) ENDIF END DO END DO END DO ELSE e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:) )*tmask(:,:,1) ENDIF ENDIF ! IF ( ln_isf ) THEN ! if ice shelf melting DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too DO jj = 1, jpj DO ji = 1, jpi IF( misfkt(ji,jj) <= jk .AND. jk <= misfkb(ji,jj) ) THEN e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) ENDIF END DO END DO END DO ENDIF ! ! IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity ! !* ssh at u- and v-points) DO jj = 2, jpjm1 ; DO ji = 2, jpim1 zsshu(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) zsshv(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) ! ! ! !* e3u and e3v z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:) ! u-point DO jk = 1, jpkm1 e3u_b (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) END DO z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:) ! v-point DO jk = 1, jpkm1 e3v_b (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) END DO ! DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) ! ub(ji,jj,jk) = zuf ! ub <-- filtered velocity vb(ji,jj,jk) = zvf un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua vn(ji,jj,jk) = va(ji,jj,jk) END DO END DO END DO ! ELSE ! Asselin filter applied on thickness weighted velocity ! ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) ! ! !* ssh at u- and v-points) DO jj = 2, jpjm1 ; DO ji = 2, jpim1 zsshu(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) zsshv(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) ! ! ! !* e3u and e3v z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:) ! u-point DO jk = 1, jpkm1 ze3u_f(:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) END DO z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:) ! v-point DO jk = 1, jpkm1 ze3u_f(:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) END DO ! DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) ! zuf = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) zvf = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) ! ub(ji,jj,jk) = zuf ! ub <-- filtered velocity vb(ji,jj,jk) = zvf un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua vn(ji,jj,jk) = va(ji,jj,jk) END DO END DO END DO e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) ! DEALLOCATE( ze3u_f , ze3v_f ) ENDIF ! ENDIF ! IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN ! Revert "before" velocities to time split estimate ! Doing it here also means that asselin filter contribution is removed zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) DO jk = 2, jpkm1 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) END DO DO jk = 1, jpkm1 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) END DO ENDIF ! ENDIF ! end Leap-Frog time stepping ! ! Set "now" and "before" barotropic velocities for next time step: ! JC: Would be more clever to swap variables than to make a full vertical integration ! ! IF(.NOT.ln_linssh ) THEN DO jj = 2, jpjm1 ; DO ji = 2, jpim1 zsshu(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) zsshv(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) ! hu_b (:,:) = hu_0(:,:) + zsshu(:,:) hv_b (:,:) = hv_0(:,:) + zsshv(:,:) r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) ENDIF ! ! ! un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) DO jk = 2, jpkm1 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) END DO un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) ! IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents CALL iom_put( "ubar", un_b(:,:) ) CALL iom_put( "vbar", vn_b(:,:) ) ENDIF IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_Dt zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_Dt CALL trd_dyn( zua, zva, jpdyn_atf, kt ) ENDIF ! IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) ! IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) IF( l_trddyn ) DEALLOCATE( zua, zva ) IF( ln_timing ) CALL timing_stop('dyn_nxt') ! END SUBROUTINE dyn_nxt !!========================================================================= END MODULE dynnxt