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

source: trunk/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 1438

Last change on this file since 1438 was 1438, checked in by rblod, 15 years ago

Merge VVL branch with the trunk (act II), see ticket #429
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.5 KB
Line 
1MODULE dynnxt
2   !!======================================================================
3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
5   !!======================================================================
6   !!======================================================================
7   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
8   !!                 !  1990-10  (C. Levy, G. Madec)
9   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
10   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
11   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
12   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
13   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
14   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
15   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
16   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
17   !!            3.2  !  2009-04  (G. Madec, R.Benshila))  re-introduce the vvl option
18   !!----------------------------------------------------------------------
19 
20   !!----------------------------------------------------------------------
21   !!   dyn_nxt      : update the horizontal velocity from the momentum trend
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE in_out_manager  ! I/O manager
26   USE obc_oce         ! ocean open boundary conditions
27   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
28   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
29   USE obcvol          ! ocean open boundary condition (obc_vol routines)
30   USE bdy_oce         ! unstructured open boundary conditions
31   USE bdydta          ! unstructured open boundary conditions
32   USE bdydyn          ! unstructured open boundary conditions
33   USE dynspg_oce      ! type of surface pressure gradient
34   USE lbclnk          ! lateral boundary condition (or mpp link)
35   USE prtctl          ! Print control
36   USE agrif_opa_update
37   USE agrif_opa_interp
38   USE domvvl          ! variable volume
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC    dyn_nxt   ! routine called by step.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47   !!-------------------------------------------------------------------------
48   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
49   !! $Id$
50   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!-------------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE dyn_nxt ( kt )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dyn_nxt  ***
58      !!                   
59      !! ** Purpose :   Compute the after horizontal velocity from the
60      !!      momentum trend.
61      !!
62      !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va)
63      !!      through calls to routine lbc_lnk.
64      !!      After velocity is compute using a leap-frog scheme environment:
65      !!         (ua,va) = (ub,vb) + 2 rdt (ua,va)
66      !!      Note that if lk_dynspg_flt=T, the time stepping has already been
67      !!      performed in dynspg module
68      !!      Time filter applied on now horizontal velocity to avoid the
69      !!      divergence of two consecutive time-steps and swap of dynamics
70      !!      arrays to start the next time step:
71      !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
72      !!         (un,vn) = (ua,va)
73      !!
74      !! ** Action : - Update ub,vb arrays, the before horizontal velocity
75      !!             - Update un,vn arrays, the now horizontal velocity
76      !!
77      !!----------------------------------------------------------------------
78      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
79      !!
80      INTEGER  ::   ji, jj, jk   ! dummy loop indices
81      REAL(wp) ::   z2dt         ! temporary scalar
82      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
83      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
84      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
85      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -
86      REAL(wp) ::   zuf   , zvf              !    -         -
87     !!----------------------------------------------------------------------
88
89      IF( kt == nit000 ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
92         IF(lwp) WRITE(numout,*) '~~~~~~~'
93      ENDIF
94
95      ! Local constant initialization
96      z2dt = 2. * rdt
97      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
98
99      !! Explicit physics with thickness weighted updates
100
101      ! Lateral boundary conditions on ( ua, va )
102      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. ) 
103
104      ! Next velocity
105      ! -------------
106#if defined key_dynspg_flt
107      ! Leap-frog time stepping already done in dynspg_flt.F routine
108#else
109      IF( lk_vvl ) THEN                                  ! Varying levels
110         DO jk = 1, jpkm1
111            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     &
112               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   &
113               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
114            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     &
115               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   &
116               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
117         END DO
118      ELSE
119         DO jk = 1, jpkm1
120                  ! Leap-frog time stepping
121            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
122            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
123         END DO
124      ENDIF
125
126# if defined key_obc
127      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
128      CALL obc_dyn( kt )
129
130      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
131         ! Flather boundary condition :
132         !        - Update sea surface height on each open boundary
133         !                 sshn (= after ssh) for explicit case
134         !                 sshn_b (= after ssha_b) for time-splitting case
135         !        - Correct the barotropic velocities
136         CALL obc_dyn_bt( kt )
137         !
138         ! Boundary conditions on sshn ( after ssh)
139         CALL lbc_lnk( sshn, 'T', 1. )
140         !
141         IF( ln_vol_cst )   CALL obc_vol( kt )
142         !
143         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
144      ENDIF
145
146# elif defined key_bdy 
147      ! Update (ua,va) along open boundaries (for exp or ts options).
148      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
149         !
150         CALL bdy_dyn_frs( kt )
151         !
152         IF( ln_bdy_fla ) THEN
153            ua_e(:,:) = 0.e0
154            va_e(:,:) = 0.e0
155            ! Set these variables for use in bdy_dyn_fla
156            hu_e(:,:) = hu(:,:)
157            hv_e(:,:) = hv(:,:)
158            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
159               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
160               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
161            END DO
162            DO jk = 1 , jpkm1
163               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
164               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
165            END DO
166            CALL bdy_dta_bt( kt+1, 0)
167            CALL bdy_dyn_fla
168         ENDIF
169         !
170         DO jk = 1 , jpkm1
171            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
172            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
173         END DO
174      ENDIF
175# endif
176
177# if defined key_agrif
178      CALL Agrif_dyn( kt )
179# endif
180#endif
181
182      ! Time filter and swap of dynamics arrays
183      ! ------------------------------------------
184      IF( neuler == 0 .AND. kt == nit000 ) THEN
185         DO jk = 1, jpkm1                                 
186            un(:,:,jk) = ua(:,:,jk)
187            vn(:,:,jk) = va(:,:,jk)
188         END DO
189      ELSE
190         IF( lk_vvl ) THEN                                ! Varying levels
191            DO jk = 1, jpkm1                              ! filter applied on thickness weighted velocities
192               DO jj = 1, jpj
193                  DO ji = 1, jpi
194                     ze3u_a = fse3u_a(ji,jj,jk)
195                     ze3v_a = fse3v_a(ji,jj,jk)
196                     ze3u_n = fse3u_n(ji,jj,jk)
197                     ze3v_n = fse3v_n(ji,jj,jk)
198                     ze3u_b = fse3u_b(ji,jj,jk)
199                     ze3v_b = fse3v_b(ji,jj,jk)
200                     !
201                     zue3a = ua(ji,jj,jk) * ze3u_a
202                     zve3a = va(ji,jj,jk) * ze3v_a
203                     zue3n = un(ji,jj,jk) * ze3u_n
204                     zve3n = vn(ji,jj,jk) * ze3v_n
205                     zue3b = ub(ji,jj,jk) * ze3u_b
206                     zve3b = vb(ji,jj,jk) * ze3v_b
207                     !
208                     ub(ji,jj,jk) = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
209                        &         / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
210                     vb(ji,jj,jk) = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
211                        &         / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
212                     un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk)
213                     vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk)
214                  END DO
215               END DO
216            END DO
217         ELSE                                             ! Fixed levels
218            DO jk = 1, jpkm1                              ! filter applied on velocities
219               DO jj = 1, jpj
220                  DO ji = 1, jpi
221                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
222                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
223                     ub(ji,jj,jk) = zuf
224                     vb(ji,jj,jk) = zvf
225                     un(ji,jj,jk) = ua(ji,jj,jk)
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      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
235#endif     
236
237      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
238         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
239      !
240   END SUBROUTINE dyn_nxt
241
242   !!======================================================================
243END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.