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 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 19.0 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      !! DCSE_NEMO: need additional directives for renamed module variables
105!FTRANS ze3u_f :I :I :z
106!FTRANS ze3v_f :I :I :z
107      !
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
109      !
110      INTEGER  ::   ji, jj, jk   ! dummy loop indices
111#if ! defined key_dynspg_flt
112      REAL(wp) ::   z2dt         ! temporary scalar
113#endif
114      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
115      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
116      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
117      !!----------------------------------------------------------------------
118
119      IF( wrk_in_use(2, 1,2,3) ) THEN
120         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
121      ENDIF
122
123      IF( kt == nit000 ) THEN
124         IF(lwp) WRITE(numout,*)
125         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
126         IF(lwp) WRITE(numout,*) '~~~~~~~'
127      ENDIF
128
129#if defined key_dynspg_flt
130      !
131      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
132      ! -------------
133
134      ! Update after velocity on domain lateral boundaries      (only local domain required)
135      ! --------------------------------------------------
136      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
137      CALL lbc_lnk( va, 'V', -1. ) 
138      !
139#else
140      ! Next velocity :   Leap-frog time stepping
141      ! -------------
142      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
143      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
144      !
145      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
146         DO jk = 1, jpkm1
147            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
148            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
149         END DO
150      ELSE                                            ! applied on thickness weighted velocity
151         DO jk = 1, jpkm1
152            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
153               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
154               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
155            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
156               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
157               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
158         END DO
159      ENDIF
160
161
162      ! Update after velocity on domain lateral boundaries
163      ! --------------------------------------------------     
164      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
165      CALL lbc_lnk( va, 'V', -1. ) 
166      !
167# if defined key_obc
168      !                                !* OBC open boundaries
169      CALL obc_dyn( kt )
170      !
171      IF( .NOT. lk_dynspg_flt ) THEN
172         ! Flather boundary condition : - Update sea surface height on each open boundary
173         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
174         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
175         !                              - Correct the barotropic velocities
176         CALL obc_dyn_bt( kt )
177         !
178!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
179         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
180         !
181         IF( ln_vol_cst )   CALL obc_vol( kt )
182         !
183         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
184      ENDIF
185      !
186# elif defined key_bdy 
187      !                                !* BDY open boundaries
188      IF( .NOT. lk_dynspg_flt ) THEN
189         CALL bdy_dyn_frs( kt )
190#  if ! defined key_vvl
191         ua_e(:,:) = 0.e0
192         va_e(:,:) = 0.e0
193         ! Set these variables for use in bdy_dyn_fla
194         hur_e(:,:) = hur(:,:)
195         hvr_e(:,:) = hvr(:,:)
196         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
197            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
198            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
199         END DO
200         ua_e(:,:) = ua_e(:,:) * hur(:,:)
201         va_e(:,:) = va_e(:,:) * hvr(:,:)
202         DO jk = 1 , jpkm1
203            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
204            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
205         END DO
206         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
207         CALL bdy_dyn_fla( sshn_b )
208         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
209         CALL lbc_lnk( va_e, 'V', -1. )   !
210         DO jk = 1 , jpkm1
211            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
212            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
213         END DO
214#  endif
215      ENDIF
216# endif
217      !
218# if defined key_agrif
219      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
220# endif
221#endif
222
223      ! Time filter and swap of dynamics arrays
224      ! ------------------------------------------
225      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
226#if defined key_z_first
227         DO jj = 1, jpj
228            DO ji = 1, jpi
229               DO jk = 1, jpkm1
230                  un(ji,jj,jk) = ua(ji,jj,jk)                ! un <-- ua
231                  vn(ji,jj,jk) = va(ji,jj,jk)
232               END DO
233            END DO
234         END DO
235#else
236         DO jk = 1, jpkm1
237            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
238            vn(:,:,jk) = va(:,:,jk)
239         END DO
240#endif
241      ELSE                                             !* Leap-Frog : Asselin filter and swap
242         !                                ! =============!
243         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
244            !                             ! =============!
245#if defined key_z_first
246            DO jj = 1, jpj
247               DO ji = 1, jpi   
248                  DO jk = 1, jpkm1                             
249#else
250            DO jk = 1, jpkm1                             
251               DO jj = 1, jpj
252                  DO ji = 1, jpi   
253#endif
254                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
255                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
256                     !
257                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
258                     vb(ji,jj,jk) = zvf
259                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
260                     vn(ji,jj,jk) = va(ji,jj,jk)
261                  END DO
262               END DO
263            END DO
264            !                             ! ================!
265         ELSE                             ! Variable volume !
266            !                             ! ================!
267            ! Before scale factor at t-points
268            ! -------------------------------
269            DO jk = 1, jpkm1
270               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
271                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
272                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
273            ENDDO
274            ! Add volume filter correction only at the first level of t-point scale factors
275            zec = atfp * rdt / rau0
276#if defined key_z_first
277            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask_1(:,:)
278#else
279            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
280#endif
281            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
282            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
283            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
284            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
285            !
286            IF( ln_dynadv_vec ) THEN
287               ! Before scale factor at (u/v)-points
288               ! -----------------------------------
289               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
290#if defined key_z_first
291               DO jj = 1, jpjm1
292                  DO ji = 1, jpim1
293                     DO jk = 1, jpkm1
294#else
295               DO jk = 1, jpkm1
296                  DO jj = 1, jpjm1
297                     DO ji = 1, jpim1
298#endif
299                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
300                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
301                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
302                        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) )
303                        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) )
304                     END DO
305                  END DO
306               END DO
307               ! lateral boundary conditions
308               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
309               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
310               ! Add initial scale factor to scale factor anomaly
311               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
312               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
313               ! Leap-Frog - Asselin filter and swap: applied on velocity
314               ! -----------------------------------
315#if defined key_z_first
316               DO jj = 1, jpj
317                  DO ji = 1, jpi
318                     DO jk = 1, jpkm1
319#else
320               DO jk = 1, jpkm1
321                  DO jj = 1, jpj
322                     DO ji = 1, jpi
323#endif
324                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
325                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
326                        !
327                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
328                        vb(ji,jj,jk) = zvf
329                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
330                        vn(ji,jj,jk) = va(ji,jj,jk)
331                     END DO
332                  END DO
333               END DO
334               !
335            ELSE
336               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
337               !-----------------------------------------------
338               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
339#if defined key_z_first
340               DO jj = 1, jpjm1
341                  DO ji = 1, jpim1
342                     DO jk = 1, jpkm1
343#else
344               DO jk = 1, jpkm1
345                  DO jj = 1, jpjm1
346                     DO ji = 1, jpim1
347#endif
348                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
349                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
350                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
351                        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) )
352                        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) )
353                     END DO
354                  END DO
355               END DO
356               ! lateral boundary conditions
357               CALL lbc_lnk( ze3u_f, 'U', 1. )
358               CALL lbc_lnk( ze3v_f, 'V', 1. )
359               ! Add initial scale factor to scale factor anomaly
360               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
361               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
362               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
363               ! -----------------------------------             ===========================
364#if defined key_z_first
365               DO jj = 1, jpj
366                  DO ji = 1, jpim1
367                     DO jk = 1, jpkm1
368#else
369               DO jk = 1, jpkm1
370                  DO jj = 1, jpj
371                     DO ji = 1, jpim1
372#endif
373                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
374                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
375                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
376                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
377                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
378                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
379                        !
380                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
381                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
382                        !
383                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
384                        vb(ji,jj,jk) = zvf
385                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
386                        vn(ji,jj,jk) = va(ji,jj,jk)
387                     END DO
388                  END DO
389               END DO
390               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
391               fse3v_b(:,:,:) = ze3v_f(:,:,:)
392               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
393               CALL lbc_lnk( vb, 'V', -1. )
394            ENDIF
395            !
396         ENDIF
397         !
398      ENDIF
399
400      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
401         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
402      !
403      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
404      !
405   END SUBROUTINE dyn_nxt
406
407   !!=========================================================================
408END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.