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

Last change on this file since 2219 was 2208, checked in by rblod, 14 years ago

Put FCM NEMO code changes in DEV_r2191_3partymerge2010 branch
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.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   !!-------------------------------------------------------------------------
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 in_out_manager  ! I/O manager
36   USE lbclnk          ! lateral boundary condition (or mpp link)
37   USE prtctl          ! Print control
38#if defined key_agrif
39   USE agrif_opa_update
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.2 , LOCEAN-IPSL (2009)
52   !! $Id$
53   !! Software is governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt)
54   !!-------------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE dyn_nxt ( kt )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE dyn_nxt  ***
61      !!                   
62      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
63      !!             condition on the after velocity, achieved the time stepping
64      !!             by applying the Asselin filter on now fields and swapping
65      !!             the fields.
66      !!
67      !! ** Method  : * After velocity is compute using a leap-frog scheme:
68      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
69      !!             Note that with flux form advection and variable volume layer
70      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
71      !!             velocity.
72      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
73      !!             the time stepping has already been done in dynspg module
74      !!
75      !!              * Apply lateral boundary conditions on after velocity
76      !!             at the local domain boundaries through lbc_lnk call,
77      !!             at the radiative open boundaries (lk_obc=T),
78      !!             at the relaxed   open boundaries (lk_bdy=T), and
79      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
80      !!
81      !!              * Apply the time filter applied and swap of the dynamics
82      !!             arrays to start the next time step:
83      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
84      !!                (un,vn) = (ua,va).
85      !!             Note that with flux form advection and variable volume layer
86      !!             (lk_vvl=T), the time filter is applied on thickness weighted
87      !!             velocity.
88      !!
89      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
90      !!               un,vn   now horizontal velocity of next time-step
91      !!----------------------------------------------------------------------
92      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
93      !!
94      INTEGER  ::   ji, jj, jk   ! dummy loop indices
95#if ! defined key_dynspg_flt
96      REAL(wp) ::   z2dt         ! temporary scalar
97#endif
98      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
99      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
100      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
101      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -
102      REAL(wp) ::   zuf   , zvf              !    -         -
103      !!----------------------------------------------------------------------
104
105      IF( kt == nit000 ) THEN
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
108         IF(lwp) WRITE(numout,*) '~~~~~~~'
109      ENDIF
110
111#if defined key_dynspg_flt
112      !
113      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
114      ! -------------
115
116      ! Update after velocity on domain lateral boundaries      (only local domain required)
117      ! --------------------------------------------------
118      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
119      CALL lbc_lnk( va, 'V', -1. ) 
120      !
121#else
122      ! Next velocity :   Leap-frog time stepping
123      ! -------------
124      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
125      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
126      !
127      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
128         DO jk = 1, jpkm1
129            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
130            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
131         END DO
132      ELSE                                            ! applied on thickness weighted velocity
133         DO jk = 1, jpkm1
134            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
135               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
136               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
137            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
138               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
139               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
140         END DO
141      ENDIF
142
143
144      ! Update after velocity on domain lateral boundaries
145      ! --------------------------------------------------     
146      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
147      CALL lbc_lnk( va, 'V', -1. ) 
148      !
149# if defined key_obc
150      !                                !* OBC open boundaries
151      IF( lk_obc )   CALL obc_dyn( kt )
152      !
153      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
154         ! Flather boundary condition : - Update sea surface height on each open boundary
155         !                                       sshn   (= after ssh   ) for explicit case
156         !                                       sshn_b (= after ssha_b) for time-splitting case
157         !                              - Correct the barotropic velocities
158         CALL obc_dyn_bt( kt )
159         !
160!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
161         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
162         !
163         IF( ln_vol_cst )   CALL obc_vol( kt )
164         !
165         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
166      ENDIF
167      !
168# elif defined key_bdy 
169      !                                !* BDY open boundaries
170      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
171         CALL bdy_dyn_frs( kt )
172      ENDIF
173# endif
174      !
175# if defined key_agrif
176      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
177# endif
178#endif
179
180      ! Time filter and swap of dynamics arrays
181      ! ------------------------------------------
182      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
183         DO jk = 1, jpkm1
184            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
185            vn(:,:,jk) = va(:,:,jk)
186         END DO
187      ELSE                                             !* Leap-Frog : Asselin filter and swap
188         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
189            DO jk = 1, jpkm1                             
190               DO jj = 1, jpj
191                  DO ji = 1, jpi   
192                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
193                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(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         ELSE                                                ! applied on thickness weighted velocity
203            DO jk = 1, jpkm1
204               DO jj = 1, jpj
205                  DO ji = 1, jpi
206                     ze3u_a = fse3u_a(ji,jj,jk)
207                     ze3v_a = fse3v_a(ji,jj,jk)
208                     ze3u_n = fse3u_n(ji,jj,jk)
209                     ze3v_n = fse3v_n(ji,jj,jk)
210                     ze3u_b = fse3u_b(ji,jj,jk)
211                     ze3v_b = fse3v_b(ji,jj,jk)
212                     !
213                     zue3a = ua(ji,jj,jk) * ze3u_a
214                     zve3a = va(ji,jj,jk) * ze3v_a
215                     zue3n = un(ji,jj,jk) * ze3u_n
216                     zve3n = vn(ji,jj,jk) * ze3v_n
217                     zue3b = ub(ji,jj,jk) * ze3u_b
218                     zve3b = vb(ji,jj,jk) * ze3v_b
219                     !
220                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
221                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
222                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
223                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
224                     !
225                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
226                     vb(ji,jj,jk) = zvf
227                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
228                     vn(ji,jj,jk) = va(ji,jj,jk)
229                  END DO
230               END DO
231            END DO
232         ENDIF
233      ENDIF
234
235#if defined key_agrif
236      ! Update velocity at AGRIF zoom boundaries
237      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
238#endif     
239
240      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
241         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
242      !
243   END SUBROUTINE dyn_nxt
244
245   !!=========================================================================
246END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.