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

Last change on this file since 3837 was 3837, checked in by trackstand2, 11 years ago

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