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/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

Implementation of the merge of TRA/TRP : first guess, see ticket #842
  • Property svn:keywords set to Id
File size: 17.6 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   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
18   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
19   !!-------------------------------------------------------------------------
20 
21   !!-------------------------------------------------------------------------
22   !!   dyn_nxt      : obtain the next (after) horizontal velocity
23   !!-------------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE phycst          ! physical constants
28   USE dynspg_oce      ! type of surface pressure gradient
29   USE dynadv          ! dynamics: vector invariant versus flux form
30   USE domvvl          ! variable volume
31   USE obc_oce         ! ocean open boundary conditions
32   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
33   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
34   USE obcvol          ! ocean open boundary condition (obc_vol routines)
35   USE bdy_oce         ! unstructured open boundary conditions
36   USE bdydta          ! unstructured open boundary conditions
37   USE bdydyn          ! unstructured open boundary conditions
38   USE in_out_manager  ! I/O manager
39   USE lbclnk          ! lateral boundary condition (or mpp link)
40   USE lib_mpp         ! MPP library
41   USE prtctl          ! Print control
42#if defined key_agrif
43   USE agrif_opa_interp
44#endif
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC    dyn_nxt   ! routine called by step.F90
50
51   !! * Substitutions
52#  include "domzgr_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
55   !! $Id$
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE dyn_nxt ( kt )
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE dyn_nxt  ***
63      !!                   
64      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
65      !!             condition on the after velocity, achieved the time stepping
66      !!             by applying the Asselin filter on now fields and swapping
67      !!             the fields.
68      !!
69      !! ** Method  : * After velocity is compute using a leap-frog scheme:
70      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
71      !!             Note that with flux form advection and variable volume layer
72      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
73      !!             velocity.
74      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
75      !!             the time stepping has already been done in dynspg module
76      !!
77      !!              * Apply lateral boundary conditions on after velocity
78      !!             at the local domain boundaries through lbc_lnk call,
79      !!             at the radiative open boundaries (lk_obc=T),
80      !!             at the relaxed   open boundaries (lk_bdy=T), and
81      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and variable volume layer
88      !!             (lk_vvl=T), the time filter is applied on thickness weighted
89      !!             velocity.
90      !!
91      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
95      USE oce     , ONLY:   tsa             ! tsa used as 2 3D workspace
96      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
97      !
98      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
99      !
100      INTEGER  ::   ji, jj, jk   ! dummy loop indices
101#if ! defined key_dynspg_flt
102      REAL(wp) ::   z2dt         ! temporary scalar
103#endif
104      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
105      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
106      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
107      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
108      !!----------------------------------------------------------------------
109
110      IF( wrk_in_use(2, 1,2,3) ) THEN
111         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
112      ENDIF
113      !
114      ze3u_f => tsa(:,:,:,1) 
115      ze3v_f => tsa(:,:,:,2) 
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
120         IF(lwp) WRITE(numout,*) '~~~~~~~'
121      ENDIF
122
123#if defined key_dynspg_flt
124      !
125      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
126      ! -------------
127
128      ! Update after velocity on domain lateral boundaries      (only local domain required)
129      ! --------------------------------------------------
130      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
131      CALL lbc_lnk( va, 'V', -1. ) 
132      !
133#else
134      ! Next velocity :   Leap-frog time stepping
135      ! -------------
136      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
137      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
138      !
139      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
142            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
143         END DO
144      ELSE                                            ! applied on thickness weighted velocity
145         DO jk = 1, jpkm1
146            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
147               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
148               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
149            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
150               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
151               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
152         END DO
153      ENDIF
154
155
156      ! Update after velocity on domain lateral boundaries
157      ! --------------------------------------------------     
158      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
159      CALL lbc_lnk( va, 'V', -1. ) 
160      !
161# if defined key_obc
162      !                                !* OBC open boundaries
163      CALL obc_dyn( kt )
164      !
165      IF( .NOT. lk_dynspg_flt ) THEN
166         ! Flather boundary condition : - Update sea surface height on each open boundary
167         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
168         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
169         !                              - Correct the barotropic velocities
170         CALL obc_dyn_bt( kt )
171         !
172!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
173         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
174         !
175         IF( ln_vol_cst )   CALL obc_vol( kt )
176         !
177         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
178      ENDIF
179      !
180# elif defined key_bdy 
181      !                                !* BDY open boundaries
182      IF( .NOT. lk_dynspg_flt ) THEN
183         CALL bdy_dyn_frs( kt )
184#  if ! defined key_vvl
185         ua_e(:,:) = 0.e0
186         va_e(:,:) = 0.e0
187         ! Set these variables for use in bdy_dyn_fla
188         hur_e(:,:) = hur(:,:)
189         hvr_e(:,:) = hvr(:,:)
190         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
191            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
192            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
193         END DO
194         ua_e(:,:) = ua_e(:,:) * hur(:,:)
195         va_e(:,:) = va_e(:,:) * hvr(:,:)
196         DO jk = 1 , jpkm1
197            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
198            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
199         END DO
200         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
201         CALL bdy_dyn_fla( sshn_b )
202         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
203         CALL lbc_lnk( va_e, 'V', -1. )   !
204         DO jk = 1 , jpkm1
205            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
206            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
207         END DO
208#  endif
209      ENDIF
210# endif
211      !
212# if defined key_agrif
213      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
214# endif
215#endif
216
217      ! Time filter and swap of dynamics arrays
218      ! ------------------------------------------
219      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
220         DO jk = 1, jpkm1
221            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
222            vn(:,:,jk) = va(:,:,jk)
223         END DO
224      ELSE                                             !* Leap-Frog : Asselin filter and swap
225         !                                ! =============!
226         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
227            !                             ! =============!
228            DO jk = 1, jpkm1                             
229               DO jj = 1, jpj
230                  DO ji = 1, jpi   
231                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
232                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
233                     !
234                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
235                     vb(ji,jj,jk) = zvf
236                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
237                     vn(ji,jj,jk) = va(ji,jj,jk)
238                  END DO
239               END DO
240            END DO
241            !                             ! ================!
242         ELSE                             ! Variable volume !
243            !                             ! ================!
244            ! Before scale factor at t-points
245            ! -------------------------------
246            DO jk = 1, jpkm1
247               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
248                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
249                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
250            ENDDO
251            ! Add volume filter correction only at the first level of t-point scale factors
252            zec = atfp * rdt / rau0
253            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
254            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
255            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
256            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
257            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
258            !
259            IF( ln_dynadv_vec ) THEN
260               ! Before scale factor at (u/v)-points
261               ! -----------------------------------
262               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
263               DO jk = 1, jpkm1
264                  DO jj = 1, jpjm1
265                     DO ji = 1, jpim1
266                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
267                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
268                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
269                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
270                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
271                     END DO
272                  END DO
273               END DO
274               ! lateral boundary conditions
275               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
276               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
277               ! Add initial scale factor to scale factor anomaly
278               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
279               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
280               ! Leap-Frog - Asselin filter and swap: applied on velocity
281               ! -----------------------------------
282               DO jk = 1, jpkm1
283                  DO jj = 1, jpj
284                     DO ji = 1, jpi
285                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
286                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
287                        !
288                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
289                        vb(ji,jj,jk) = zvf
290                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
291                        vn(ji,jj,jk) = va(ji,jj,jk)
292                     END DO
293                  END DO
294               END DO
295               !
296            ELSE
297               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
298               !-----------------------------------------------
299               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
300               DO jk = 1, jpkm1
301                  DO jj = 1, jpjm1
302                     DO ji = 1, jpim1
303                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
304                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
305                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
306                        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) )
307                        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) )
308                     END DO
309                  END DO
310               END DO
311               ! lateral boundary conditions
312               CALL lbc_lnk( ze3u_f, 'U', 1. )
313               CALL lbc_lnk( ze3v_f, 'V', 1. )
314               ! Add initial scale factor to scale factor anomaly
315               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
316               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
317               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
318               ! -----------------------------------             ===========================
319               DO jk = 1, jpkm1
320                  DO jj = 1, jpj
321                     DO ji = 1, jpim1
322                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
323                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
324                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
325                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
326                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
327                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
328                        !
329                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
330                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
331                        !
332                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
333                        vb(ji,jj,jk) = zvf
334                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
335                        vn(ji,jj,jk) = va(ji,jj,jk)
336                     END DO
337                  END DO
338               END DO
339               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
340               fse3v_b(:,:,:) = ze3v_f(:,:,:)
341               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
342               CALL lbc_lnk( vb, 'V', -1. )
343            ENDIF
344            !
345         ENDIF
346         !
347      ENDIF
348
349      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
350         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
351      !
352      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
353      !
354   END SUBROUTINE dyn_nxt
355
356   !!=========================================================================
357END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.