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

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

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