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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 16.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   !!-------------------------------------------------------------------------
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   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE dyn_nxt ( kt )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE dyn_nxt  ***
62      !!                   
63      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
64      !!             condition on the after velocity, achieved the time stepping
65      !!             by applying the Asselin filter on now fields and swapping
66      !!             the fields.
67      !!
68      !! ** Method  : * After velocity is compute using a leap-frog scheme:
69      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
70      !!             Note that with flux form advection and variable volume layer
71      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
72      !!             velocity.
73      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
74      !!             the time stepping has already been done in dynspg module
75      !!
76      !!              * Apply lateral boundary conditions on after velocity
77      !!             at the local domain boundaries through lbc_lnk call,
78      !!             at the radiative open boundaries (lk_obc=T),
79      !!             at the relaxed   open boundaries (lk_bdy=T), and
80      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
81      !!
82      !!              * Apply the time filter applied and swap of the dynamics
83      !!             arrays to start the next time step:
84      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
85      !!                (un,vn) = (ua,va).
86      !!             Note that with flux form advection and variable volume layer
87      !!             (lk_vvl=T), the time filter is applied on thickness weighted
88      !!             velocity.
89      !!
90      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
91      !!               un,vn   now horizontal velocity of next time-step
92      !!----------------------------------------------------------------------
93      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
94      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
95      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
96      !
97      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
98      !
99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
100#if ! defined key_dynspg_flt
101      REAL(wp) ::   z2dt         ! temporary scalar
102#endif
103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
105      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
106      !!----------------------------------------------------------------------
107
108      IF( wrk_in_use(2, 1,2,3) ) THEN
109         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
110      ENDIF
111
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
118#if defined key_dynspg_flt
119      !
120      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
121      ! -------------
122
123      ! Update after velocity on domain lateral boundaries      (only local domain required)
124      ! --------------------------------------------------
125      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
126      CALL lbc_lnk( va, 'V', -1. ) 
127      !
128#else
129      ! Next velocity :   Leap-frog time stepping
130      ! -------------
131      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
132      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
133      !
134      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
135         DO jk = 1, jpkm1
136            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
137            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
138         END DO
139      ELSE                                            ! applied on thickness weighted velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
142               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
143               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
144            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
145               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
146               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
147         END DO
148      ENDIF
149
150
151      ! Update after velocity on domain lateral boundaries
152      ! --------------------------------------------------     
153      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
154      CALL lbc_lnk( va, 'V', -1. ) 
155      !
156# if defined key_obc
157      !                                !* OBC open boundaries
158      CALL obc_dyn( kt )
159      !
160      IF( .NOT. lk_dynspg_flt ) THEN
161         ! Flather boundary condition : - Update sea surface height on each open boundary
162         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
163         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
164         !                              - Correct the barotropic velocities
165         CALL obc_dyn_bt( kt )
166         !
167!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
168         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
169         !
170         IF( ln_vol_cst )   CALL obc_vol( kt )
171         !
172         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
173      ENDIF
174      !
175# elif defined key_bdy 
176      !                                !* BDY open boundaries
177      IF( .NOT. lk_dynspg_flt )   CALL bdy_dyn_frs( kt )
178# endif
179      !
180# if defined key_agrif
181      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
182# endif
183#endif
184
185      ! Time filter and swap of dynamics arrays
186      ! ------------------------------------------
187      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
188         DO jk = 1, jpkm1
189            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
190            vn(:,:,jk) = va(:,:,jk)
191         END DO
192      ELSE                                             !* Leap-Frog : Asselin filter and swap
193         !                                ! =============!
194         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
195            !                             ! =============!
196            DO jk = 1, jpkm1                             
197               DO jj = 1, jpj
198                  DO ji = 1, jpi   
199                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
200                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
201                     !
202                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
203                     vb(ji,jj,jk) = zvf
204                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
205                     vn(ji,jj,jk) = va(ji,jj,jk)
206                  END DO
207               END DO
208            END DO
209            !                             ! ================!
210         ELSE                             ! Variable volume !
211            !                             ! ================!
212            ! Before scale factor at t-points
213            ! -------------------------------
214            DO jk = 1, jpkm1
215               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
216                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
217                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
218            ENDDO
219            ! Add volume filter correction only at the first level of t-point scale factors
220            zec = atfp * rdt / rau0
221            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
222            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
223            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
224            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
225            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
226            !
227            IF( ln_dynadv_vec ) THEN
228               ! Before scale factor at (u/v)-points
229               ! -----------------------------------
230               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
231               DO jk = 1, jpkm1
232                  DO jj = 1, jpjm1
233                     DO ji = 1, jpim1
234                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
235                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
236                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
237                        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) )
238                        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) )
239                     END DO
240                  END DO
241               END DO
242               ! lateral boundary conditions
243               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
244               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
245               ! Add initial scale factor to scale factor anomaly
246               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
247               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
248               ! Leap-Frog - Asselin filter and swap: applied on velocity
249               ! -----------------------------------
250               DO jk = 1, jpkm1
251                  DO jj = 1, jpj
252                     DO ji = 1, jpi
253                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
254                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
255                        !
256                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
257                        vb(ji,jj,jk) = zvf
258                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
259                        vn(ji,jj,jk) = va(ji,jj,jk)
260                     END DO
261                  END DO
262               END DO
263               !
264            ELSE
265               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
266               !-----------------------------------------------
267               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
268               DO jk = 1, jpkm1
269                  DO jj = 1, jpjm1
270                     DO ji = 1, jpim1
271                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
272                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
273                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
274                        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) )
275                        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) )
276                     END DO
277                  END DO
278               END DO
279               ! lateral boundary conditions
280               CALL lbc_lnk( ze3u_f, 'U', 1. )
281               CALL lbc_lnk( ze3v_f, 'V', 1. )
282               ! Add initial scale factor to scale factor anomaly
283               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
284               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
285               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
286               ! -----------------------------------             ===========================
287               DO jk = 1, jpkm1
288                  DO jj = 1, jpj
289                     DO ji = 1, jpim1
290                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
291                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
292                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
293                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
294                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
295                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
296                        !
297                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
298                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
299                        !
300                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
301                        vb(ji,jj,jk) = zvf
302                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
303                        vn(ji,jj,jk) = va(ji,jj,jk)
304                     END DO
305                  END DO
306               END DO
307               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
308               fse3v_b(:,:,:) = ze3v_f(:,:,:)
309               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
310               CALL lbc_lnk( vb, 'V', -1. )
311            ENDIF
312            !
313         ENDIF
314         !
315      ENDIF
316
317      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
318         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
319      !
320      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
321      !
322   END SUBROUTINE dyn_nxt
323
324   !!=========================================================================
325END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.