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 !!------------------------------------------------------------------------- !!------------------------------------------------------------------------- !! dyn_nxt : obtain the next (after) horizontal velocity !!------------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE dynspg_oce ! type of surface pressure gradient USE dynadv ! dynamics: vector invariant versus flux form USE domvvl ! variable volume USE obc_oce ! ocean open boundary conditions USE obcdyn ! open boundary condition for momentum (obc_dyn routine) USE obcdyn_bt ! 2D open boundary condition for momentum (obc_dyn_bt routine) USE obcvol ! ocean open boundary condition (obc_vol routines) USE bdy_oce ! unstructured open boundary conditions USE bdydta ! unstructured open boundary conditions USE bdydyn ! unstructured open boundary conditions USE agrif_opa_update USE agrif_opa_interp USE in_out_manager ! I/O manager USE lbclnk ! lateral boundary condition (or mpp link) USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC dyn_nxt ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" !!------------------------------------------------------------------------- !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) !! $Id$ !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!------------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_nxt ( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_nxt *** !! !! ** Purpose : Compute the after horizontal velocity. Apply the boundary !! condition on the after velocity, achieved the time stepping !! by applying the Asselin filter on now fields and swapping !! the fields. !! !! ** Method : * After velocity is compute using a leap-frog scheme: !! (ua,va) = (ub,vb) + 2 rdt (ua,va) !! Note that with flux form advection and variable volume layer !! (lk_vvl=T), the leap-frog is applied on thickness weighted !! velocity. !! Note also that in filtered free surface (lk_dynspg_flt=T), !! the time stepping has already been done in dynspg module !! !! * Apply lateral boundary conditions on after velocity !! at the local domain boundaries through lbc_lnk call, !! at the radiative open boundaries (lk_obc=T), !! at the relaxed open boundaries (lk_bdy=T), and !! 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) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] !! (un,vn) = (ua,va). !! Note that with flux form advection and variable volume layer !! (lk_vvl=T), the time filter is applied on thickness weighted !! velocity. !! !! ** 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 #if ! defined key_dynspg_flt REAL(wp) :: z2dt ! temporary scalar #endif REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar REAL(wp) :: zve3a , zve3n , zve3b ! - - REAL(wp) :: ze3u_b, ze3u_n, ze3u_a ! - - REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - - REAL(wp) :: zuf , zvf ! - - !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~' ENDIF #if defined key_dynspg_flt ! ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine ! ------------- ! Update after velocity on domain lateral boundaries (only local domain required) ! -------------------------------------------------- CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries CALL lbc_lnk( va, 'V', -1. ) ! #else ! Next velocity : Leap-frog time stepping ! ------------- z2dt = 2. * rdt ! Euler or leap-frog time step IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity DO jk = 1, jpkm1 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) END DO ELSE ! applied on thickness weighted velocity DO jk = 1, jpkm1 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & & / fse3u_a(:,:,jk) * umask(:,:,jk) va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & & / fse3v_a(:,:,jk) * vmask(:,:,jk) END DO ENDIF ! Update after velocity on domain lateral boundaries ! -------------------------------------------------- CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries CALL lbc_lnk( va, 'V', -1. ) ! # if defined key_obc ! !* OBC open boundaries IF( lk_obc ) CALL obc_dyn( kt ) ! IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN ! Flather boundary condition : - Update sea surface height on each open boundary ! sshn (= after ssh ) for explicit case ! sshn_b (= after ssha_b) for time-splitting case ! - Correct the barotropic velocities CALL obc_dyn_bt( kt ) ! !!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called CALL lbc_lnk( sshn, 'T', 1. ) ! Boundary conditions on sshn ! IF( ln_vol_cst ) CALL obc_vol( kt ) ! IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) ENDIF ! # elif defined key_bdy ! !* BDY open boundaries !RB all this part should be in a specific routine IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN ! except for filtered option ! CALL bdy_dyn_frs( kt ) ! IF( ln_bdy_dyn_fla ) THEN ua_e(:,:) = 0.e0 va_e(:,:) = 0.e0 ! Set these variables for use in bdy_dyn_fla hur_e(:,:) = hur(:,:) hvr_e(:,:) = hvr(:,:) DO jk = 1, jpkm1 !! Vertically integrated momentum trends ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) END DO ua_e(:,:) = ua_e(:,:) * hur(:,:) va_e(:,:) = va_e(:,:) * hvr(:,:) DO jk = 1 , jpkm1 ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) va(:,:,jk) = va(:,:,jk) - va_e(:,:) END DO CALL bdy_dta_bt( kt+1, 0) CALL bdy_dyn_fla( sshn_b ) CALL lbc_lnk( ua_e, 'U', -1. ) ! Boundary points should be updated CALL lbc_lnk( va_e, 'V', -1. ) ! DO jk = 1 , jpkm1 ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk) va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk) END DO ENDIF ! ENDIF # endif ! # if defined key_agrif CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries # endif #endif ! Time filter and swap of dynamics arrays ! ------------------------------------------ IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap DO jk = 1, jpkm1 un(:,:,jk) = ua(:,:,jk) ! un <-- ua vn(:,:,jk) = va(:,:,jk) END DO ELSE !* Leap-Frog : Asselin filter and swap IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(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 ! applied on thickness weighted velocity DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ze3u_a = fse3u_a(ji,jj,jk) ze3v_a = fse3v_a(ji,jj,jk) ze3u_n = fse3u_n(ji,jj,jk) ze3v_n = fse3v_n(ji,jj,jk) ze3u_b = fse3u_b(ji,jj,jk) ze3v_b = fse3v_b(ji,jj,jk) ! zue3a = ua(ji,jj,jk) * ze3u_a zve3a = va(ji,jj,jk) * ze3v_a zue3n = un(ji,jj,jk) * ze3u_n zve3n = vn(ji,jj,jk) * ze3v_n zue3b = ub(ji,jj,jk) * ze3u_b zve3b = vb(ji,jj,jk) * ze3v_b ! zuf = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) zvf = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(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 ENDIF ENDIF #if defined key_agrif ! Update velocity at AGRIF zoom boundaries IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn( kt ) #endif IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) ! END SUBROUTINE dyn_nxt !!========================================================================= END MODULE dynnxt