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

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

Bug fixes for the dynspg_ts case.
  • Property svn:keywords set to Id
File size: 15.6 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 obcdta          ! ocean open boundary conditions
33   USE obcdyn3d        ! ocean open boundary conditions
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!!$!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
160!!$         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
161!!$         !
162!!$         IF( ln_vol_cst )   CALL obc_vol( kt )
163!!$         !
164!!$         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
165
166      ENDIF
167      !
168# endif
169      !
170# if defined key_agrif
171      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
172# endif
173#endif
174
175      ! Time filter and swap of dynamics arrays
176      ! ------------------------------------------
177      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
178         DO jk = 1, jpkm1
179            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
180            vn(:,:,jk) = va(:,:,jk)
181         END DO
182      ELSE                                             !* Leap-Frog : Asselin filter and swap
183         !                                ! =============!
184         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
185            !                             ! =============!
186            DO jk = 1, jpkm1                             
187               DO jj = 1, jpj
188                  DO ji = 1, jpi   
189                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
190                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
191                     !
192                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
193                     vb(ji,jj,jk) = zvf
194                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
195                     vn(ji,jj,jk) = va(ji,jj,jk)
196                  END DO
197               END DO
198            END DO
199            !                             ! ================!
200         ELSE                             ! Variable volume !
201            !                             ! ================!
202            ! Before scale factor at t-points
203            ! -------------------------------
204            DO jk = 1, jpkm1
205               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
206                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
207                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
208            ENDDO
209            ! Add volume filter correction only at the first level of t-point scale factors
210            zec = atfp * rdt / rau0
211            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
212            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
213            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
214            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
215            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
216            !
217            IF( ln_dynadv_vec ) THEN
218               ! Before scale factor at (u/v)-points
219               ! -----------------------------------
220               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
221               DO jk = 1, jpkm1
222                  DO jj = 1, jpjm1
223                     DO ji = 1, jpim1
224                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
225                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
226                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
227                        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) )
228                        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) )
229                     END DO
230                  END DO
231               END DO
232               ! lateral boundary conditions
233               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
234               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
235               ! Add initial scale factor to scale factor anomaly
236               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
237               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
238               ! Leap-Frog - Asselin filter and swap: applied on velocity
239               ! -----------------------------------
240               DO jk = 1, jpkm1
241                  DO jj = 1, jpj
242                     DO ji = 1, jpi
243                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
244                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
245                        !
246                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
247                        vb(ji,jj,jk) = zvf
248                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
249                        vn(ji,jj,jk) = va(ji,jj,jk)
250                     END DO
251                  END DO
252               END DO
253               !
254            ELSE
255               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
256               !-----------------------------------------------
257               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
258               DO jk = 1, jpkm1
259                  DO jj = 1, jpjm1
260                     DO ji = 1, jpim1
261                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
262                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
263                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
264                        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) )
265                        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) )
266                     END DO
267                  END DO
268               END DO
269               ! lateral boundary conditions
270               CALL lbc_lnk( ze3u_f, 'U', 1. )
271               CALL lbc_lnk( ze3v_f, 'V', 1. )
272               ! Add initial scale factor to scale factor anomaly
273               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
274               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
275               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
276               ! -----------------------------------             ===========================
277               DO jk = 1, jpkm1
278                  DO jj = 1, jpj
279                     DO ji = 1, jpim1
280                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
281                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
282                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
283                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
284                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
285                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
286                        !
287                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
288                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
289                        !
290                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
291                        vb(ji,jj,jk) = zvf
292                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
293                        vn(ji,jj,jk) = va(ji,jj,jk)
294                     END DO
295                  END DO
296               END DO
297               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
298               fse3v_b(:,:,:) = ze3v_f(:,:,:)
299               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
300               CALL lbc_lnk( vb, 'V', -1. )
301            ENDIF
302            !
303         ENDIF
304         !
305      ENDIF
306
307      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
308         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
309      !
310      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
311      !
312   END SUBROUTINE dyn_nxt
313
314   !!=========================================================================
315END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.