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 branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN – NEMO

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2005

Last change on this file since 2005 was 2005, checked in by mlelod, 14 years ago

ticket: #663 MLF: second part (local compatibility essentially)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.7 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_update
35   USE agrif_opa_interp
36   USE in_out_manager  ! I/O manager
37   USE lbclnk          ! lateral boundary condition (or mpp link)
38   USE prtctl          ! Print control
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC    dyn_nxt   ! routine called by step.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47   !!-------------------------------------------------------------------------
48   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
49   !! $Id$
50   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!-------------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE dyn_nxt ( kt )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dyn_nxt  ***
58      !!                   
59      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
60      !!             condition on the after velocity, achieved the time stepping
61      !!             by applying the Asselin filter on now fields and swapping
62      !!             the fields.
63      !!
64      !! ** Method  : * After velocity is compute using a leap-frog scheme:
65      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
66      !!             Note that with flux form advection and variable volume layer
67      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
68      !!             velocity.
69      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
70      !!             the time stepping has already been done in dynspg module
71      !!
72      !!              * Apply lateral boundary conditions on after velocity
73      !!             at the local domain boundaries through lbc_lnk call,
74      !!             at the radiative open boundaries (lk_obc=T),
75      !!             at the relaxed   open boundaries (lk_bdy=T), and
76      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
77      !!
78      !!              * Apply the time filter applied and swap of the dynamics
79      !!             arrays to start the next time step:
80      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
81      !!                (un,vn) = (ua,va).
82      !!             Note that with flux form advection and variable volume layer
83      !!             (lk_vvl=T), the time filter is applied on thickness weighted
84      !!             velocity.
85      !!
86      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
87      !!               un,vn   now horizontal velocity of next time-step
88      !!----------------------------------------------------------------------
89#if defined key_vvl
90      USE oce, ONLY :   ze3u_f => ta   ! use ta as 3D workspace
91      USE oce, ONLY :   ze3v_f => sa   ! use sa as 3D workspace
92#endif
93      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
94      !!
95      INTEGER  ::   ji, jj, jk   ! dummy loop indices
96#if ! defined key_dynspg_flt
97      REAL(wp) ::   z2dt         ! temporary scalar
98#endif
99      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
100      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
101      REAL(wp) ::   zuf   , zvf              !    -         -
102      REAL(wp) ::   zec                      !    -         -
103      REAL(wp) ::   zv_t_ij  , zv_t_ip1j     !     -        -
104      REAL(wp) ::   zv_t_ijp1                !     -        -
105      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      ! temporary 2D workspace
106      !!----------------------------------------------------------------------
107
108      IF( kt == nit000 ) THEN
109         IF(lwp) WRITE(numout,*)
110         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
111         IF(lwp) WRITE(numout,*) '~~~~~~~'
112      ENDIF
113
114#if defined key_dynspg_flt
115      !
116      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
117      ! -------------
118
119      ! Update after velocity on domain lateral boundaries      (only local domain required)
120      ! --------------------------------------------------
121      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
122      CALL lbc_lnk( va, 'V', -1. ) 
123      !
124#else
125      ! Next velocity :   Leap-frog time stepping
126      ! -------------
127      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
128      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
129      !
130      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
131         DO jk = 1, jpkm1
132            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
133            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
134         END DO
135      ELSE                                            ! applied on thickness weighted velocity
136         DO jk = 1, jpkm1
137            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
138               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
139               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
140            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
141               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
142               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
143         END DO
144      ENDIF
145
146
147      ! Update after velocity on domain lateral boundaries
148      ! --------------------------------------------------     
149      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
150      CALL lbc_lnk( va, 'V', -1. ) 
151      !
152# if defined key_obc
153      !                                !* OBC open boundaries
154      CALL obc_dyn( kt )
155      !
156      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
157         ! Flather boundary condition : - Update sea surface height on each open boundary
158         !                                       sshn   (= after ssh   ) for explicit case
159         !                                       sshn_b (= after ssha_b) for time-splitting case
160         !                              - Correct the barotropic velocities
161         CALL obc_dyn_bt( kt )
162         !
163!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
164         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
165         !
166         IF( ln_vol_cst )   CALL obc_vol( kt )
167         !
168         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
169      ENDIF
170      !
171# elif defined key_bdy 
172      !                                !* BDY open boundaries
173      !RB all this part should be in a specific routine
174      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
175         !
176         CALL bdy_dyn_frs( kt )
177         !
178         IF( ln_bdy_dyn_fla ) THEN
179            ua_e(:,:) = 0.e0
180            va_e(:,:) = 0.e0
181            ! Set these variables for use in bdy_dyn_fla
182            hur_e(:,:) = hur(:,:)
183            hvr_e(:,:) = hvr(:,:)
184            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
185               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
186               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
187            END DO
188            ua_e(:,:) = ua_e(:,:) * hur(:,:)
189            va_e(:,:) = va_e(:,:) * hvr(:,:)
190            DO jk = 1 , jpkm1
191               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
192               va(:,:,jk) = va(:,:,jk) - va_e(:,:)
193            END DO
194            CALL bdy_dta_bt( kt+1, 0)
195            CALL bdy_dyn_fla( sshn_b )
196            CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
197            CALL lbc_lnk( va_e, 'V', -1. )   !
198            DO jk = 1 , jpkm1
199               ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
200               va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
201            END DO
202         ENDIF
203         !
204      ENDIF
205# endif
206      !
207# if defined key_agrif
208      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
209# endif
210#endif
211
212      ! Time filter and swap of dynamics arrays
213      ! ------------------------------------------
214      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
215         DO jk = 1, jpkm1
216            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
217            vn(:,:,jk) = va(:,:,jk)
218         END DO
219      ELSE                                             !* Leap-Frog : Asselin filter and swap
220         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
221            DO jk = 1, jpkm1                             
222               DO jj = 1, jpj
223                  DO ji = 1, jpi   
224                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
225                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
226                     !
227                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
228                     vb(ji,jj,jk) = zvf
229                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
230                     vn(ji,jj,jk) = va(ji,jj,jk)
231                  END DO
232               END DO
233            END DO
234         ELSE                                                ! applied on thickness weighted velocity
235            ! Before scale factors at (t/u/v)-points (actually "now filtered" and futur "before")
236            ! ======================================
237            ! Scale factor at t-points
238            ! ------------------------
239            fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * fse3t_m(:,:,:)
240            ! Add volume filter correction only at the first level of t-point scale factors
241            zec = atfp * rdt / rau0
242            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
243            ! Scale factor at (u/v)-points
244            ! ------------------------
245            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
246            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
247            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
248            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
249            ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
250            DO jk = 1, jpkm1
251               DO jj = 1, jpjm1
252                  DO ji = 1, jpim1
253                     zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
254                     zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
255                     zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
256                     ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
257                     ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
258                  END DO
259               END DO
260            END DO
261            CALL lbc_lnk( ze3u_f, 'U', 1. )               ! lateral boundary conditions
262            CALL lbc_lnk( ze3v_f, 'U', 1. )
263            ! Add initial scale factor to scale factor anomaly
264            ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
265            ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
266           
267            DO jk = 1, jpkm1
268               DO jj = 1, jpj
269                  DO ji = 1, jpim
270                     zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
271                     zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
272                     zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
273                     zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
274                     zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
275                     zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
276                     !
277                     zuf  = ( zue3n + atfp * ( zue3b  - 2.e0 * zue3n  + zue3a  ) ) / ze3u_f(ji,jj,jk)
278                     zvf  = ( zve3n + atfp * ( zve3b  - 2.e0 * zve3n  + zve3a  ) ) / ze3v_f(ji,jj,jk)
279                     !
280                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
281                     vb(ji,jj,jk) = zvf
282                     fse3u_b(ji,jj,jk) = ze3u_f(ji,jj,jk)    ! e3u_b <-- filtered scale factor
283                     fse3v_b(ji,jj,jk) = ze3v_f(ji,jj,jk)
284                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
285                     vn(ji,jj,jk) = va(ji,jj,jk)
286                  END DO
287               END DO
288            END DO
289            CALL lbc_lnk( ub, 'U', -1. )         ! local domain boundaries
290            CALL lbc_lnk( vb, 'V', -1. ) 
291         ENDIF
292      ENDIF
293
294#if defined key_agrif
295      ! Update velocity at AGRIF zoom boundaries
296      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
297#endif     
298
299      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
300         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
301      !
302   END SUBROUTINE dyn_nxt
303
304   !!=========================================================================
305END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.