New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynnxt.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2524

Last change on this file since 2524 was 2486, checked in by rblod, 14 years ago

Correct Agrif inconstency for ssh, nemo_v3_2 version, see ticket 669
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.5 KB
Line 
1MODULE dynnxt
2   !!=========================================================================
3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
5   !!=========================================================================
6   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
7   !!                 !  1990-10  (C. Levy, G. Madec)
8   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
9   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
10   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
11   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
12   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
14   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
17   !!-------------------------------------------------------------------------
18 
19   !!-------------------------------------------------------------------------
20   !!   dyn_nxt      : obtain the next (after) horizontal velocity
21   !!-------------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE dynspg_oce      ! type of surface pressure gradient
25   USE dynadv          ! dynamics: vector invariant versus flux form
26   USE domvvl          ! variable volume
27   USE obc_oce         ! ocean open boundary conditions
28   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
29   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
30   USE obcvol          ! ocean open boundary condition (obc_vol routines)
31   USE bdy_oce         ! unstructured open boundary conditions
32   USE bdydta          ! unstructured open boundary conditions
33   USE bdydyn          ! unstructured open boundary conditions
34   USE agrif_opa_interp
35   USE in_out_manager  ! I/O manager
36   USE lbclnk          ! lateral boundary condition (or mpp link)
37   USE prtctl          ! Print control
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC    dyn_nxt   ! routine called by step.F90
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46   !!-------------------------------------------------------------------------
47   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
48   !! $Id$
49   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!-------------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE dyn_nxt ( kt )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE dyn_nxt  ***
57      !!                   
58      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
59      !!             condition on the after velocity, achieved the time stepping
60      !!             by applying the Asselin filter on now fields and swapping
61      !!             the fields.
62      !!
63      !! ** Method  : * After velocity is compute using a leap-frog scheme:
64      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
65      !!             Note that with flux form advection and variable volume layer
66      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
67      !!             velocity.
68      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
69      !!             the time stepping has already been done in dynspg module
70      !!
71      !!              * Apply lateral boundary conditions on after velocity
72      !!             at the local domain boundaries through lbc_lnk call,
73      !!             at the radiative open boundaries (lk_obc=T),
74      !!             at the relaxed   open boundaries (lk_bdy=T), and
75      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
76      !!
77      !!              * Apply the time filter applied and swap of the dynamics
78      !!             arrays to start the next time step:
79      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
80      !!                (un,vn) = (ua,va).
81      !!             Note that with flux form advection and variable volume layer
82      !!             (lk_vvl=T), the time filter is applied on thickness weighted
83      !!             velocity.
84      !!
85      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
86      !!               un,vn   now horizontal velocity of next time-step
87      !!----------------------------------------------------------------------
88      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
89      !!
90      INTEGER  ::   ji, jj, jk   ! dummy loop indices
91#if ! defined key_dynspg_flt
92      REAL(wp) ::   z2dt         ! temporary scalar
93#endif
94      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
95      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
96      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
97      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -
98      REAL(wp) ::   zuf   , zvf              !    -         -
99      !!----------------------------------------------------------------------
100
101      IF( kt == nit000 ) THEN
102         IF(lwp) WRITE(numout,*)
103         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
104         IF(lwp) WRITE(numout,*) '~~~~~~~'
105      ENDIF
106
107#if defined key_dynspg_flt
108      !
109      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
110      ! -------------
111
112      ! Update after velocity on domain lateral boundaries      (only local domain required)
113      ! --------------------------------------------------
114      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
115      CALL lbc_lnk( va, 'V', -1. ) 
116      !
117#else
118      ! Next velocity :   Leap-frog time stepping
119      ! -------------
120      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
121      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
122      !
123      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
124         DO jk = 1, jpkm1
125            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
126            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
127         END DO
128      ELSE                                            ! applied on thickness weighted velocity
129         DO jk = 1, jpkm1
130            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
131               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
132               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
133            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
134               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
135               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
136         END DO
137      ENDIF
138
139
140      ! Update after velocity on domain lateral boundaries
141      ! --------------------------------------------------     
142      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
143      CALL lbc_lnk( va, 'V', -1. ) 
144      !
145# if defined key_obc
146      !                                !* OBC open boundaries
147      IF( lk_obc )   CALL obc_dyn( kt )
148      !
149      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
150         ! Flather boundary condition : - Update sea surface height on each open boundary
151         !                                       sshn   (= after ssh   ) for explicit case
152         !                                       sshn_b (= after ssha_b) for time-splitting case
153         !                              - Correct the barotropic velocities
154         CALL obc_dyn_bt( kt )
155         !
156!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
157         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
158         !
159         IF( ln_vol_cst )   CALL obc_vol( kt )
160         !
161         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
162      ENDIF
163      !
164# elif defined key_bdy 
165      !                                !* BDY open boundaries
166      !RB all this part should be in a specific routine
167      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
168         !
169         CALL bdy_dyn_frs( kt )
170         !
171         IF( ln_bdy_dyn_fla ) THEN
172            ua_e(:,:) = 0.e0
173            va_e(:,:) = 0.e0
174            ! Set these variables for use in bdy_dyn_fla
175            hur_e(:,:) = hur(:,:)
176            hvr_e(:,:) = hvr(:,:)
177            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
178               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
179               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
180            END DO
181            ua_e(:,:) = ua_e(:,:) * hur(:,:)
182            va_e(:,:) = va_e(:,:) * hvr(:,:)
183            DO jk = 1 , jpkm1
184               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
185               va(:,:,jk) = va(:,:,jk) - va_e(:,:)
186            END DO
187            CALL bdy_dta_bt( kt+1, 0)
188            CALL bdy_dyn_fla( sshn_b )
189            CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
190            CALL lbc_lnk( va_e, 'V', -1. )   !
191            DO jk = 1 , jpkm1
192               ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
193               va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
194            END DO
195         ENDIF
196         !
197      ENDIF
198# endif
199      !
200# if defined key_agrif
201      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
202# endif
203#endif
204
205      ! Time filter and swap of dynamics arrays
206      ! ------------------------------------------
207      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
208         DO jk = 1, jpkm1
209            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
210            vn(:,:,jk) = va(:,:,jk)
211         END DO
212      ELSE                                             !* Leap-Frog : Asselin filter and swap
213         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
214            DO jk = 1, jpkm1                             
215               DO jj = 1, jpj
216                  DO ji = 1, jpi   
217                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
218                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
219                     !
220                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
221                     vb(ji,jj,jk) = zvf
222                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
223                     vn(ji,jj,jk) = va(ji,jj,jk)
224                  END DO
225               END DO
226            END DO
227         ELSE                                                ! applied on thickness weighted velocity
228            DO jk = 1, jpkm1
229               DO jj = 1, jpj
230                  DO ji = 1, jpi
231                     ze3u_a = fse3u_a(ji,jj,jk)
232                     ze3v_a = fse3v_a(ji,jj,jk)
233                     ze3u_n = fse3u_n(ji,jj,jk)
234                     ze3v_n = fse3v_n(ji,jj,jk)
235                     ze3u_b = fse3u_b(ji,jj,jk)
236                     ze3v_b = fse3v_b(ji,jj,jk)
237                     !
238                     zue3a = ua(ji,jj,jk) * ze3u_a
239                     zve3a = va(ji,jj,jk) * ze3v_a
240                     zue3n = un(ji,jj,jk) * ze3u_n
241                     zve3n = vn(ji,jj,jk) * ze3v_n
242                     zue3b = ub(ji,jj,jk) * ze3u_b
243                     zve3b = vb(ji,jj,jk) * ze3v_b
244                     !
245                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
246                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
247                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
248                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
249                     !
250                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
251                     vb(ji,jj,jk) = zvf
252                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
253                     vn(ji,jj,jk) = va(ji,jj,jk)
254                  END DO
255               END DO
256            END DO
257         ENDIF
258      ENDIF
259
260      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
261         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
262      !
263   END SUBROUTINE dyn_nxt
264
265   !!=========================================================================
266END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.