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

Last change on this file since 4416 was 4416, checked in by trackstand2, 7 years ago

Added debug SUMs to dyn_nxt

  • Property svn:keywords set to Id
File size: 19.3 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      WRITE (*,*) 'ARPDBG: dyn_nxt: after lbc update sum(ua)=',SUM(ua),' at step=',kt
171# if defined key_obc
172      !                                !* OBC open boundaries
173      CALL obc_dyn( kt )
174      !
175      IF( .NOT. lk_dynspg_flt ) THEN
176         ! Flather boundary condition : - Update sea surface height on each open boundary
177         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
178         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
179         !                              - Correct the barotropic velocities
180         CALL obc_dyn_bt( kt )
181         !
182!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
183         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
184         !
185         IF( ln_vol_cst )   CALL obc_vol( kt )
186         !
187         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
188      ENDIF
189      !
190# elif defined key_bdy 
191      !                                !* BDY open boundaries
192      IF( .NOT. lk_dynspg_flt ) THEN
193         CALL bdy_dyn_frs( kt )
194#  if ! defined key_vvl
195         ua_e(:,:) = 0.e0
196         va_e(:,:) = 0.e0
197         ! Set these variables for use in bdy_dyn_fla
198         hur_e(:,:) = hur(:,:)
199         hvr_e(:,:) = hvr(:,:)
200         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
201            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
202            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
203         END DO
204         ua_e(:,:) = ua_e(:,:) * hur(:,:)
205         va_e(:,:) = va_e(:,:) * hvr(:,:)
206         DO jk = 1 , jpkm1
207            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
208            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
209         END DO
210         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
211         CALL bdy_dyn_fla( sshn_b )
212         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
213         CALL lbc_lnk( va_e, 'V', -1. )   !
214         DO jk = 1 , jpkm1
215            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
216            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
217         END DO
218#  endif
219      ENDIF
220# endif
221      !
222# if defined key_agrif
223      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
224# endif
225#endif
226
227      WRITE(*,*) 'ARPDBG: dyn_nxt b4 time filter SUM(ua)=',SUM(ua),'at step=',kt
228
229      ! Time filter and swap of dynamics arrays
230      ! ------------------------------------------
231      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
232#if defined key_z_first
233         DO jj = 1, jpj
234            DO ji = 1, jpi
235               DO jk = 1, jpkm1
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#else
242         DO jk = 1, jpkm1
243            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
244            vn(:,:,jk) = va(:,:,jk)
245         END DO
246#endif
247      ELSE                                             !* Leap-Frog : Asselin filter and swap
248         !                                ! =============!
249         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
250            !                             ! =============!
251#if defined key_z_first
252            DO jj = 1, jpj
253               DO ji = 1, jpi   
254                  DO jk = 1, jpkm1                             
255#else
256            DO jk = 1, jpkm1                             
257               DO jj = 1, jpj
258                  DO ji = 1, jpi   
259#endif
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                             ! Variable volume !
272            !                             ! ================!
273            ! Before scale factor at t-points
274            ! -------------------------------
275            DO jk = 1, jpkm1
276               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
277                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
278                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
279            ENDDO
280            ! Add volume filter correction only at the first level of t-point scale factors
281            zec = atfp * rdt / rau0
282#if defined key_z_first
283            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask_1(:,:)
284#else
285            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
286#endif
287            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
288            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
289            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
290            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
291            !
292            IF( ln_dynadv_vec ) THEN
293               ! Before scale factor at (u/v)-points
294               ! -----------------------------------
295               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
296#if defined key_z_first
297               DO jj = 1, jpjm1
298                  DO ji = 1, jpim1
299                     DO jk = 1, jpkm1
300#else
301               DO jk = 1, jpkm1
302                  DO jj = 1, jpjm1
303                     DO ji = 1, jpim1
304#endif
305                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
306                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
307                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
308                        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) )
309                        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) )
310                     END DO
311                  END DO
312               END DO
313               ! lateral boundary conditions
314               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
315               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
316               ! Add initial scale factor to scale factor anomaly
317               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
318               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
319               ! Leap-Frog - Asselin filter and swap: applied on velocity
320               ! -----------------------------------
321#if defined key_z_first
322               DO jj = 1, jpj
323                  DO ji = 1, jpi
324                     DO jk = 1, jpkm1
325#else
326               DO jk = 1, jpkm1
327                  DO jj = 1, jpj
328                     DO ji = 1, jpi
329#endif
330                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
331                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
332                        !
333                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
334                        vb(ji,jj,jk) = zvf
335                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
336                        vn(ji,jj,jk) = va(ji,jj,jk)
337                     END DO
338                  END DO
339               END DO
340               !
341            ELSE
342               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
343               !-----------------------------------------------
344               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
345#if defined key_z_first
346               DO jj = 1, jpjm1
347                  DO ji = 1, jpim1
348                     DO jk = 1, jpkm1
349#else
350               DO jk = 1, jpkm1
351                  DO jj = 1, jpjm1
352                     DO ji = 1, jpim1
353#endif
354                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
355                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
356                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
357                        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) )
358                        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) )
359                     END DO
360                  END DO
361               END DO
362               ! lateral boundary conditions
363               CALL lbc_lnk( ze3u_f, 'U', 1. )
364               CALL lbc_lnk( ze3v_f, 'V', 1. )
365               ! Add initial scale factor to scale factor anomaly
366               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
367               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
368               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
369               ! -----------------------------------             ===========================
370#if defined key_z_first
371               DO jj = 1, jpj
372                  DO ji = 1, jpim1
373                     DO jk = 1, jpkm1
374#else
375               DO jk = 1, jpkm1
376                  DO jj = 1, jpj
377                     DO ji = 1, jpim1
378#endif
379                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
380                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
381                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
382                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
383                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
384                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
385                        !
386                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
387                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
388                        !
389                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
390                        vb(ji,jj,jk) = zvf
391                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
392                        vn(ji,jj,jk) = va(ji,jj,jk)
393                     END DO
394                  END DO
395               END DO
396               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
397               fse3v_b(:,:,:) = ze3v_f(:,:,:)
398               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
399               CALL lbc_lnk( vb, 'V', -1. )
400            ENDIF
401            !
402         ENDIF
403         !
404      ENDIF
405
406      WRITE(*,*) 'ARPDBG: dyn_nxt end SUM(un)=',SUM(un),'at step=',kt
407
408      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
409         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
410      !
411      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
412      !
413   END SUBROUTINE dyn_nxt
414
415   !!=========================================================================
416END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.