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

source: branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2207

Last change on this file since 2207 was 2207, checked in by acc, 14 years ago

#733 DEV_r2191_3partymerge2010. Merged in changes from devukmo2010 branch
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.5 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 dynspg_oce      ! type of surface pressure gradient
26   USE dynadv          ! dynamics: vector invariant versus flux form
27   USE domvvl          ! variable volume
28   USE obc_oce         ! ocean open boundary conditions
29   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
30   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
31   USE obcvol          ! ocean open boundary condition (obc_vol routines)
32   USE bdy_oce         ! unstructured open boundary conditions
33   USE bdydta          ! unstructured open boundary conditions
34   USE bdydyn          ! unstructured open boundary conditions
35   USE agrif_opa_update
36   USE agrif_opa_interp
37   USE in_out_manager  ! I/O manager
38   USE lbclnk          ! lateral boundary condition (or mpp link)
39   USE prtctl          ! Print control
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC    dyn_nxt   ! routine called by step.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!-------------------------------------------------------------------------
49   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
50   !! $Id$
51   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
52   !!-------------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE dyn_nxt ( kt )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dyn_nxt  ***
59      !!                   
60      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
61      !!             condition on the after velocity, achieved the time stepping
62      !!             by applying the Asselin filter on now fields and swapping
63      !!             the fields.
64      !!
65      !! ** Method  : * After velocity is compute using a leap-frog scheme:
66      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
67      !!             Note that with flux form advection and variable volume layer
68      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
69      !!             velocity.
70      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
71      !!             the time stepping has already been done in dynspg module
72      !!
73      !!              * Apply lateral boundary conditions on after velocity
74      !!             at the local domain boundaries through lbc_lnk call,
75      !!             at the radiative open boundaries (lk_obc=T),
76      !!             at the relaxed   open boundaries (lk_bdy=T), and
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      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
91      !!
92      INTEGER  ::   ji, jj, jk   ! dummy loop indices
93#if ! defined key_dynspg_flt
94      REAL(wp) ::   z2dt         ! temporary scalar
95#endif
96      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
97      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
98      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
99      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -
100      REAL(wp) ::   zuf   , zvf              !    -         -
101      !!----------------------------------------------------------------------
102
103      IF( kt == nit000 ) THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
106         IF(lwp) WRITE(numout,*) '~~~~~~~'
107      ENDIF
108
109#if defined key_dynspg_flt
110      !
111      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
112      ! -------------
113
114      ! Update after velocity on domain lateral boundaries      (only local domain required)
115      ! --------------------------------------------------
116      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
117      CALL lbc_lnk( va, 'V', -1. ) 
118      !
119#else
120      ! Next velocity :   Leap-frog time stepping
121      ! -------------
122      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
123      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
124      !
125      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
126         DO jk = 1, jpkm1
127            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
128            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
129         END DO
130      ELSE                                            ! applied on thickness weighted velocity
131         DO jk = 1, jpkm1
132            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
133               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
134               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
135            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
136               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
137               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
138         END DO
139      ENDIF
140
141
142      ! Update after velocity on domain lateral boundaries
143      ! --------------------------------------------------     
144      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
145      CALL lbc_lnk( va, 'V', -1. ) 
146      !
147# if defined key_obc
148      !                                !* OBC open boundaries
149      IF( lk_obc )   CALL obc_dyn( kt )
150      !
151      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
152         ! Flather boundary condition : - Update sea surface height on each open boundary
153         !                                       sshn   (= after ssh   ) for explicit case
154         !                                       sshn_b (= after ssha_b) for time-splitting case
155         !                              - Correct the barotropic velocities
156         CALL obc_dyn_bt( kt )
157         !
158!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
159         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
160         !
161         IF( ln_vol_cst )   CALL obc_vol( kt )
162         !
163         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
164      ENDIF
165      !
166# elif defined key_bdy 
167      !                                !* BDY open boundaries
168      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
169         CALL bdy_dyn_frs( kt )
170      ENDIF
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         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
187            DO jk = 1, jpkm1                             
188               DO jj = 1, jpj
189                  DO ji = 1, jpi   
190                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
191                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
192                     !
193                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
194                     vb(ji,jj,jk) = zvf
195                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
196                     vn(ji,jj,jk) = va(ji,jj,jk)
197                  END DO
198               END DO
199            END DO
200         ELSE                                                ! applied on thickness weighted velocity
201            DO jk = 1, jpkm1
202               DO jj = 1, jpj
203                  DO ji = 1, jpi
204                     ze3u_a = fse3u_a(ji,jj,jk)
205                     ze3v_a = fse3v_a(ji,jj,jk)
206                     ze3u_n = fse3u_n(ji,jj,jk)
207                     ze3v_n = fse3v_n(ji,jj,jk)
208                     ze3u_b = fse3u_b(ji,jj,jk)
209                     ze3v_b = fse3v_b(ji,jj,jk)
210                     !
211                     zue3a = ua(ji,jj,jk) * ze3u_a
212                     zve3a = va(ji,jj,jk) * ze3v_a
213                     zue3n = un(ji,jj,jk) * ze3u_n
214                     zve3n = vn(ji,jj,jk) * ze3v_n
215                     zue3b = ub(ji,jj,jk) * ze3u_b
216                     zve3b = vb(ji,jj,jk) * ze3v_b
217                     !
218                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
219                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
220                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
221                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
222                     !
223                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
224                     vb(ji,jj,jk) = zvf
225                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
226                     vn(ji,jj,jk) = va(ji,jj,jk)
227                  END DO
228               END DO
229            END DO
230         ENDIF
231      ENDIF
232
233#if defined key_agrif
234      ! Update velocity at AGRIF zoom boundaries
235      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
236#endif     
237
238      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
239         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
240      !
241   END SUBROUTINE dyn_nxt
242
243   !!=========================================================================
244END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.