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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
  • Property svn:keywords set to Id
File size: 16.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   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
18   !!-------------------------------------------------------------------------
19 
20   !!-------------------------------------------------------------------------
21   !!   dyn_nxt      : obtain the next (after) horizontal velocity
22   !!-------------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE sbc_oce         ! Surface boundary condition: ocean fields
26   USE phycst          ! physical constants
27   USE dynspg_oce      ! type of surface pressure gradient
28   USE dynadv          ! dynamics: vector invariant versus flux form
29   USE domvvl          ! variable volume
30   USE obc_oce         ! ocean open boundary conditions
31   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
32   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
33   USE obcvol          ! ocean open boundary condition (obc_vol routines)
34   USE bdy_oce         ! unstructured open boundary conditions
35   USE bdydta          ! unstructured open boundary conditions
36   USE bdydyn          ! unstructured open boundary conditions
37   USE in_out_manager  ! I/O manager
38   USE lbclnk          ! lateral boundary condition (or mpp link)
39   USE lib_mpp         ! MPP library
40   USE prtctl          ! Print control
41#if defined key_agrif
42   USE agrif_opa_interp
43#endif
44
45   IMPLICIT NONE
46   PRIVATE
47
48   PUBLIC    dyn_nxt   ! routine called by step.F90
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52   !!-------------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
56   !!-------------------------------------------------------------------------
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 oce, ONLY :   ze3u_f => ta   ! use ta as 3D workspace
95      USE oce, ONLY :   ze3v_f => sa   ! use sa as 3D workspace
96      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
97      USE wrk_nemo, ONLY:   zs_t => wrk_2d_1, zs_u_1 => wrk_2d_2, &
98                          zs_v_1 => wrk_2d_3
99      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
100      !!
101      INTEGER  ::   ji, jj, jk   ! dummy loop indices
102#if ! defined key_dynspg_flt
103      REAL(wp) ::   z2dt         ! temporary scalar
104#endif
105      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
106      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
107      REAL(wp) ::   zuf   , zvf              !    -         -
108      REAL(wp) ::   zec                      !    -         -
109      REAL(wp) ::   zv_t_ij  , zv_t_ip1j     !     -        -
110      REAL(wp) ::   zv_t_ijp1                !     -        -
111      !!----------------------------------------------------------------------
112
113      IF( wrk_in_use(2, 1,2,3) ) THEN
114         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
115      ENDIF
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 ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
166         ! Flather boundary condition : - Update sea surface height on each open boundary
167         !                                       sshn   (= after ssh   ) for explicit case
168         !                                       sshn_b (= after ssha_b) for time-splitting case
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( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
183         CALL bdy_dyn_frs( kt )
184      ENDIF
185# endif
186      !
187# if defined key_agrif
188      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
189# endif
190#endif
191
192      ! Time filter and swap of dynamics arrays
193      ! ------------------------------------------
194      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
195         DO jk = 1, jpkm1
196            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
197            vn(:,:,jk) = va(:,:,jk)
198         END DO
199      ELSE                                             !* Leap-Frog : Asselin filter and swap
200         !                                ! =============!
201         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
202            !                             ! =============!
203            DO jk = 1, jpkm1                             
204               DO jj = 1, jpj
205                  DO ji = 1, jpi   
206                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
207                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
208                     !
209                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
210                     vb(ji,jj,jk) = zvf
211                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
212                     vn(ji,jj,jk) = va(ji,jj,jk)
213                  END DO
214               END DO
215            END DO
216            !                             ! ================!
217         ELSE                             ! Variable volume !
218            !                             ! ================!
219            ! Before scale factor at t-points
220            ! -------------------------------
221            DO jk = 1, jpkm1
222               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
223                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
224                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
225            ENDDO
226            ! Add volume filter correction only at the first level of t-point scale factors
227            zec = atfp * rdt / rau0
228            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
229            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
230            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
231            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
232            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
233            !
234            IF( ln_dynadv_vec ) THEN
235               ! Before scale factor at (u/v)-points
236               ! -----------------------------------
237               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
238               DO jk = 1, jpkm1
239                  DO jj = 1, jpjm1
240                     DO ji = 1, jpim1
241                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
242                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
243                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
244                        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) )
245                        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) )
246                     END DO
247                  END DO
248               END DO
249               ! lateral boundary conditions
250               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
251               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
252               ! Add initial scale factor to scale factor anomaly
253               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
254               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
255               ! Leap-Frog - Asselin filter and swap: applied on velocity
256               ! -----------------------------------
257               DO jk = 1, jpkm1
258                  DO jj = 1, jpj
259                     DO ji = 1, jpi
260                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
261                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
262                        !
263                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
264                        vb(ji,jj,jk) = zvf
265                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
266                        vn(ji,jj,jk) = va(ji,jj,jk)
267                     END DO
268                  END DO
269               END DO
270               !
271            ELSE
272               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
273               !-----------------------------------------------
274               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
275               DO jk = 1, jpkm1
276                  DO jj = 1, jpjm1
277                     DO ji = 1, jpim1
278                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
279                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
280                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
281                        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) )
282                        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) )
283                     END DO
284                  END DO
285               END DO
286               ! lateral boundary conditions
287               CALL lbc_lnk( ze3u_f, 'U', 1. )
288               CALL lbc_lnk( ze3v_f, 'V', 1. )
289               ! Add initial scale factor to scale factor anomaly
290               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
291               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
292               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
293               ! -----------------------------------             ===========================
294               DO jk = 1, jpkm1
295                  DO jj = 1, jpj
296                     DO ji = 1, jpim1
297                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
298                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
299                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
300                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
301                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
302                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
303                        !
304                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
305                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
306                        !
307                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
308                        vb(ji,jj,jk) = zvf
309                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
310                        vn(ji,jj,jk) = va(ji,jj,jk)
311                     END DO
312                  END DO
313               END DO
314               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
315               fse3v_b(:,:,:) = ze3v_f(:,:,:)
316               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
317               CALL lbc_lnk( vb, 'V', -1. )
318            ENDIF
319            !
320         ENDIF
321         !
322      ENDIF
323
324      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
325         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
326      !
327      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays.')
328      !
329   END SUBROUTINE dyn_nxt
330
331   !!=========================================================================
332END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.