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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 4438

Last change on this file since 4438 was 4438, checked in by trackstand2, 10 years ago

Protect debug SUM output in dynnxt and sshwzv with ARPDBGSUM
  • Property svn:keywords set to Id
File size: 19.4 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   !! * Control permutation of array indices
52#  include "oce_ftrans.h90"
53#  include "dom_oce_ftrans.h90"
54#  include "sbc_oce_ftrans.h90"
55#  include "domvvl_ftrans.h90"
56#  include "obc_oce_ftrans.h90"
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE dyn_nxt ( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE dyn_nxt  ***
70      !!                   
71      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
72      !!             condition on the after velocity, achieved the time stepping
73      !!             by applying the Asselin filter on now fields and swapping
74      !!             the fields.
75      !!
76      !! ** Method  : * After velocity is compute using a leap-frog scheme:
77      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
78      !!             Note that with flux form advection and variable volume layer
79      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
80      !!             velocity.
81      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
82      !!             the time stepping has already been done in dynspg module
83      !!
84      !!              * Apply lateral boundary conditions on after velocity
85      !!             at the local domain boundaries through lbc_lnk call,
86      !!             at the radiative open boundaries (lk_obc=T),
87      !!             at the relaxed   open boundaries (lk_bdy=T), and
88      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
89      !!
90      !!              * Apply the time filter applied and swap of the dynamics
91      !!             arrays to start the next time step:
92      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
93      !!                (un,vn) = (ua,va).
94      !!             Note that with flux form advection and variable volume layer
95      !!             (lk_vvl=T), the time filter is applied on thickness weighted
96      !!             velocity.
97      !!
98      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
99      !!               un,vn   now horizontal velocity of next time-step
100      !!----------------------------------------------------------------------
101      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
102      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
103      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
104      USE arpdebugging, ONLY: dump_array
105      !! DCSE_NEMO: need additional directives for renamed module variables
106!FTRANS ze3u_f :I :I :z
107!FTRANS ze3v_f :I :I :z
108      !
109      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
110      !
111      INTEGER  ::   ji, jj, jk   ! dummy loop indices
112#if ! defined key_dynspg_flt
113      REAL(wp) ::   z2dt         ! temporary scalar
114#endif
115      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
116      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
117      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
118      !!----------------------------------------------------------------------
119
120      IF( wrk_in_use(2, 1,2,3) ) THEN
121         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
122      ENDIF
123
124      IF( kt == nit000 ) THEN
125         IF(lwp) WRITE(numout,*)
126         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
127         IF(lwp) WRITE(numout,*) '~~~~~~~'
128      ENDIF
129
130!      CALL dump_array(kt, 'ua_nxt_start',ua(:,:,1),withHalos=.TRUE.)
131
132#if defined key_dynspg_flt
133      !
134      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
135      ! -------------
136
137      ! Update after velocity on domain lateral boundaries      (only local domain required)
138      ! --------------------------------------------------
139      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
140      CALL lbc_lnk( va, 'V', -1. ) 
141      !
142#else
143      ! Next velocity :   Leap-frog time stepping
144      ! -------------
145      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
146      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
147      !
148      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
149         DO jk = 1, jpkm1
150            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
151            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
152         END DO
153      ELSE                                            ! applied on thickness weighted velocity
154         DO jk = 1, jpkm1
155            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
156               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
157               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
158            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
159               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
160               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
161         END DO
162      ENDIF
163
164
165      ! Update after velocity on domain lateral boundaries
166      ! --------------------------------------------------     
167      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
168      CALL lbc_lnk( va, 'V', -1. ) 
169      !
170#if defined ARPDBGSUM
171      WRITE (*,*) 'ARPDBG: dyn_nxt: after lbc update sum(ua)=',SUM(ua),' at step=',kt
172#endif
173
174# if defined key_obc
175      !                                !* OBC open boundaries
176      CALL obc_dyn( kt )
177      !
178      IF( .NOT. lk_dynspg_flt ) THEN
179         ! Flather boundary condition : - Update sea surface height on each open boundary
180         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
181         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
182         !                              - Correct the barotropic velocities
183         CALL obc_dyn_bt( kt )
184         !
185!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
186         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
187         !
188         IF( ln_vol_cst )   CALL obc_vol( kt )
189         !
190         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
191      ENDIF
192      !
193# elif defined key_bdy 
194      !                                !* BDY open boundaries
195      IF( .NOT. lk_dynspg_flt ) THEN
196         CALL bdy_dyn_frs( kt )
197#  if ! defined key_vvl
198         ua_e(:,:) = 0.e0
199         va_e(:,:) = 0.e0
200         ! Set these variables for use in bdy_dyn_fla
201         hur_e(:,:) = hur(:,:)
202         hvr_e(:,:) = hvr(:,:)
203         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
204            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
205            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
206         END DO
207         ua_e(:,:) = ua_e(:,:) * hur(:,:)
208         va_e(:,:) = va_e(:,:) * hvr(:,:)
209         DO jk = 1 , jpkm1
210            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
211            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
212         END DO
213         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
214         CALL bdy_dyn_fla( sshn_b )
215         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
216         CALL lbc_lnk( va_e, 'V', -1. )   !
217         DO jk = 1 , jpkm1
218            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
219            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
220         END DO
221#  endif
222      ENDIF
223# endif
224      !
225# if defined key_agrif
226      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
227# endif
228#endif
229
230#if defined ARPDBGSUM
231      WRITE(*,*) 'ARPDBG: dyn_nxt b4 time filter SUM(ua)=',SUM(ua),'at step=',kt
232#endif
233
234      ! Time filter and swap of dynamics arrays
235      ! ------------------------------------------
236      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
237#if defined key_z_first
238         DO jj = 1, jpj
239            DO ji = 1, jpi
240               DO jk = 1, jpkm1
241                  un(ji,jj,jk) = ua(ji,jj,jk)                ! un <-- ua
242                  vn(ji,jj,jk) = va(ji,jj,jk)
243               END DO
244            END DO
245         END DO
246#else
247         DO jk = 1, jpkm1
248            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
249            vn(:,:,jk) = va(:,:,jk)
250         END DO
251#endif
252      ELSE                                             !* Leap-Frog : Asselin filter and swap
253         !                                ! =============!
254         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
255            !                             ! =============!
256#if defined key_z_first
257            DO jj = 1, jpj
258               DO ji = 1, jpi   
259                  DO jk = 1, jpkm1                             
260#else
261            DO jk = 1, jpkm1                             
262               DO jj = 1, jpj
263                  DO ji = 1, jpi   
264#endif
265                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
266                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
267                     !
268                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
269                     vb(ji,jj,jk) = zvf
270                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
271                     vn(ji,jj,jk) = va(ji,jj,jk)
272                  END DO
273               END DO
274            END DO
275            !                             ! ================!
276         ELSE                             ! Variable volume !
277            !                             ! ================!
278            ! Before scale factor at t-points
279            ! -------------------------------
280            DO jk = 1, jpkm1
281               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
282                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
283                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
284            ENDDO
285            ! Add volume filter correction only at the first level of t-point scale factors
286            zec = atfp * rdt / rau0
287#if defined key_z_first
288            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask_1(:,:)
289#else
290            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
291#endif
292            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
293            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
294            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
295            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
296            !
297            IF( ln_dynadv_vec ) THEN
298               ! Before scale factor at (u/v)-points
299               ! -----------------------------------
300               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
301#if defined key_z_first
302               DO jj = 1, jpjm1
303                  DO ji = 1, jpim1
304                     DO jk = 1, jpkm1
305#else
306               DO jk = 1, jpkm1
307                  DO jj = 1, jpjm1
308                     DO ji = 1, jpim1
309#endif
310                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
311                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
312                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
313                        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) )
314                        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) )
315                     END DO
316                  END DO
317               END DO
318               ! lateral boundary conditions
319               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
320               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
321               ! Add initial scale factor to scale factor anomaly
322               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
323               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
324               ! Leap-Frog - Asselin filter and swap: applied on velocity
325               ! -----------------------------------
326#if defined key_z_first
327               DO jj = 1, jpj
328                  DO ji = 1, jpi
329                     DO jk = 1, jpkm1
330#else
331               DO jk = 1, jpkm1
332                  DO jj = 1, jpj
333                     DO ji = 1, jpi
334#endif
335                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
336                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
337                        !
338                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
339                        vb(ji,jj,jk) = zvf
340                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
341                        vn(ji,jj,jk) = va(ji,jj,jk)
342                     END DO
343                  END DO
344               END DO
345               !
346            ELSE
347               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
348               !-----------------------------------------------
349               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
350#if defined key_z_first
351               DO jj = 1, jpjm1
352                  DO ji = 1, jpim1
353                     DO jk = 1, jpkm1
354#else
355               DO jk = 1, jpkm1
356                  DO jj = 1, jpjm1
357                     DO ji = 1, jpim1
358#endif
359                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
360                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
361                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
362                        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) )
363                        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) )
364                     END DO
365                  END DO
366               END DO
367               ! lateral boundary conditions
368               CALL lbc_lnk( ze3u_f, 'U', 1. )
369               CALL lbc_lnk( ze3v_f, 'V', 1. )
370               ! Add initial scale factor to scale factor anomaly
371               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
372               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
373               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
374               ! -----------------------------------             ===========================
375#if defined key_z_first
376               DO jj = 1, jpj
377                  DO ji = 1, jpim1
378                     DO jk = 1, jpkm1
379#else
380               DO jk = 1, jpkm1
381                  DO jj = 1, jpj
382                     DO ji = 1, jpim1
383#endif
384                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
385                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
386                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
387                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
388                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
389                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
390                        !
391                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
392                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
393                        !
394                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
395                        vb(ji,jj,jk) = zvf
396                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
397                        vn(ji,jj,jk) = va(ji,jj,jk)
398                     END DO
399                  END DO
400               END DO
401               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
402               fse3v_b(:,:,:) = ze3v_f(:,:,:)
403               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
404               CALL lbc_lnk( vb, 'V', -1. )
405            ENDIF
406            !
407         ENDIF
408         !
409      ENDIF
410
411#if defined ARPDBGSUM
412      WRITE(*,*) 'ARPDBG: dyn_nxt end SUM(un)=',SUM(un),'at step=',kt
413#endif
414
415      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
416         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
417      !
418      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
419      !
420   END SUBROUTINE dyn_nxt
421
422   !!=========================================================================
423END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.