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

Last change on this file since 911 was 911, checked in by ctlod, 16 years ago

Implementation of the BDY package, see ticket: #126

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.4 KB
Line 
1MODULE dynnxt
2   !!======================================================================
3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
5   !!======================================================================
6   
7   !!----------------------------------------------------------------------
8   !!   dyn_nxt      : update the horizontal velocity from the momentum trend
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers
12   USE dom_oce         ! ocean space and time domain
13   USE in_out_manager  ! I/O manager
14   USE obc_oce         ! ocean open boundary conditions
15   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
16   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
17   USE obcvol          ! ocean open boundary condition (obc_vol routines)
18   USE bdy_oce         ! unstructured open boundary conditions
19   USE bdydta          ! unstructured open boundary conditions
20   USE bdydyn          ! unstructured open boundary conditions
21   USE dynspg_oce      ! type of surface pressure gradient
22   USE lbclnk          ! lateral boundary condition (or mpp link)
23   USE prtctl          ! Print control
24   USE agrif_opa_update
25   USE agrif_opa_interp
26   USE domvvl          ! variable volume
27
28   IMPLICIT NONE
29   PRIVATE
30
31   !! * Accessibility
32   PUBLIC dyn_nxt                ! routine called by step.F90
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE dyn_nxt ( kt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dyn_nxt  ***
42      !!                   
43      !! ** Purpose :   Compute the after horizontal velocity from the
44      !!      momentum trend.
45      !!
46      !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va)
47      !!      through calls to routine lbc_lnk.
48      !!      After velocity is compute using a leap-frog scheme environment:
49      !!         (ua,va) = (ub,vb) + 2 rdt (ua,va)
50      !!      Note that if lk_dynspg_flt=T, the time stepping has already been
51      !!      performed in dynspg module
52      !!      Time filter applied on now horizontal velocity to avoid the
53      !!      divergence of two consecutive time-steps and swap of dynamics
54      !!      arrays to start the next time step:
55      !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
56      !!         (un,vn) = (ua,va)
57      !!
58      !! ** Action : - Update ub,vb arrays, the before horizontal velocity
59      !!             - Update un,vn arrays, the now horizontal velocity
60      !!
61      !! History :
62      !!        !  87-02  (P. Andrich, D. L Hostis)  Original code
63      !!        !  90-10  (C. Levy, G. Madec)
64      !!        !  93-03  (M. Guyon)  symetrical conditions
65      !!        !  97-02  (G. Madec & M. Imbard)  opa, release 8.0
66      !!        !  97-04  (A. Weaver)  Euler forward step
67      !!        !  97-06  (G. Madec)  lateral boudary cond., lbc routine
68      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
69      !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
70      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
71      !!    "   !  07-07  (D. Storkey) Calls to BDY routines.
72      !!----------------------------------------------------------------------
73      !! * Arguments
74      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
75
76      !! * Local declarations
77      INTEGER  ::   ji, jj, jk   ! dummy loop indices
78      REAL(wp) ::   z2dt         ! temporary scalar
79      REAL(wp) ::   zsshun1, zsshvn1         ! temporary scalar
80      !! Variable volume
81      REAL(wp), DIMENSION(jpi,jpj)     ::  & ! 2D workspace
82         zsshub, zsshun, zsshua,           &
83         zsshvb, zsshvn, zsshva
84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
85         zfse3ub, zfse3un, zfse3ua, &        ! 3D workspace
86         zfse3vb, zfse3vn, zfse3va
87      !!----------------------------------------------------------------------
88      !!  OPA 9.0 , LOCEAN-IPSL (2005)
89      !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DYN/dynnxt.F90,v 1.13 2007/05/25 15:51:50 opalod Exp $
90      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
91      !!----------------------------------------------------------------------
92
93      IF( kt == nit000 ) THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
96         IF(lwp) WRITE(numout,*) '~~~~~~~'
97      ENDIF
98
99      ! Local constant initialization
100      z2dt = 2. * rdt
101      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
102
103      !! Explicit physics with thickness weighted updates
104      IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
105
106         ! Sea surface elevation time stepping
107         ! -----------------------------------
108         !
109         DO jj = 1, jpjm1
110            DO ji = 1,jpim1
111
112               ! Sea Surface Height at u-point before
113               zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
114                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
115                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * sshbb(ji+1,jj  ) )
116
117               ! Sea Surface Height at v-point before
118               zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
119                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
120                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * sshbb(ji  ,jj+1) )
121
122               ! Sea Surface Height at u-point after
123               zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
124                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
125                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * ssha(ji+1,jj  ) )
126
127               ! Sea Surface Height at v-point after
128               zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
129                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
130                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * ssha(ji  ,jj+1) )
131
132            END DO
133         END DO
134
135         ! Boundaries conditions
136         CALL lbc_lnk( zsshua, 'U', 1. )   ;    CALL lbc_lnk( zsshva, 'V', 1. )
137         CALL lbc_lnk( zsshub, 'U', 1. )   ;    CALL lbc_lnk( zsshvb, 'V', 1. )
138
139         ! Scale factors at before and after time step
140         ! -------------------------------------------
141         CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 
142         CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 
143
144         ! Asselin filtered scale factor at now time step
145         ! ----------------------------------------------
146         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
147            CALL dom_vvl_sf_ini( 'U', zfse3un ) ;   CALL dom_vvl_sf_ini( 'V', zfse3vn ) 
148         ELSE
149            zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)
150            zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)
151            CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ;   CALL dom_vvl_sf( zsshvn, 'V', zfse3vn ) 
152         ENDIF
153
154         ! Thickness weighting
155         ! -------------------
156         DO jk = 1, jpkm1
157            DO jj = 1, jpj
158               DO ji = 1, jpi
159                  ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
160                  va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
161
162                  un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)
163                  vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)
164
165                  ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
166                  vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
167               END DO
168            END DO
169         END DO
170
171      ENDIF
172
173      ! Lateral boundary conditions on ( ua, va )
174      CALL lbc_lnk( ua, 'U', -1. )
175      CALL lbc_lnk( va, 'V', -1. )
176
177      !                                                ! ===============
178      DO jk = 1, jpkm1                                 ! Horizontal slab
179         !                                             ! ===============
180         ! Next velocity
181         ! -------------
182#if defined key_dynspg_flt
183         ! Leap-frog time stepping already done in dynspg.F routine
184#else
185         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
186            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
187               ! Leap-frog time stepping
188               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
189               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
190            END DO
191         END DO
192
193         IF( lk_vvl ) THEN
194            ! Unweight velocities prior to updating open boundaries.
195
196            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
197               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
198                  ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk)
199                  va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk)
200
201                  un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
202                  vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
203
204                  ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk)
205                  vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk)
206               END DO
207            END DO
208
209         ENDIF
210
211# if defined key_obc
212         !                                             ! ===============
213      END DO                                           !   End of slab
214      !                                                ! ===============
215      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
216      CALL obc_dyn( kt )
217
218      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
219         !Flather boundary condition :
220         !        - Update sea surface height on each open boundary
221         !                 sshn (= after ssh) for explicit case
222         !                 sshn_b (= after ssha_b) for time-splitting case
223         !        - Correct the barotropic velocities
224         CALL obc_dyn_bt( kt )
225
226         !Boundary conditions on sshn ( after ssh)
227         CALL lbc_lnk( sshn, 'T', 1. )
228
229         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
230            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
231         ENDIF
232
233         IF ( ln_vol_cst ) CALL obc_vol( kt )
234
235      ENDIF
236
237      !                                                ! ===============
238      DO jk = 1, jpkm1                                 ! Horizontal slab
239         !                                             ! ===============
240# elif defined key_bdy || key_bdy_tides
241         !                                             ! ===============
242      END DO                                           !   End of slab
243      !                                                ! ===============
244      ! Update (ua,va) along open boundaries (for exp or ts options).
245      IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
246 
247         CALL bdy_dyn_frs( kt )
248
249         IF ( ln_bdy_fla ) THEN
250
251            ua_e(:,:)=0.0
252            va_e(:,:)=0.0
253
254            ! Set these variables for use in bdy_dyn_fla
255            hu_e(:,:) = hu(:,:)
256            hv_e(:,:) = hv(:,:)
257
258            DO jk = 1, jpkm1
259               !! Vertically integrated momentum trends
260               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
261               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
262            END DO
263
264            DO jk = 1 , jpkm1
265               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
266               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
267            END DO
268
269            CALL bdy_dta_bt( kt+1, 0)
270            CALL bdy_dyn_fla
271
272         ENDIF
273
274         DO jk = 1 , jpkm1
275            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
276            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
277         END DO
278
279      ENDIF
280
281      !                                                ! ===============
282      DO jk = 1, jpkm1                                 ! Horizontal slab
283         !                                             ! ===============
284# endif
285# if defined key_agrif
286         !                                             ! ===============
287      END DO                                           !   End of slab
288      !                                                ! ===============
289      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
290      CALL Agrif_dyn( kt )
291      !                                                ! ===============
292      DO jk = 1, jpkm1                                 ! Horizontal slab
293         !                                             ! ===============
294# endif
295#endif
296
297         ! Time filter and swap of dynamics arrays
298         ! ------------------------------------------
299         IF( neuler == 0 .AND. kt == nit000 ) THEN
300            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
301               ! caution: don't use (:,:) for this loop
302               ! it causes optimization problems on NEC in auto-tasking
303               DO jj = 1, jpj
304                  DO ji = 1, jpi
305                     zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk)
306                     zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
307                     ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
308                     vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
309                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
310                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
311                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
312                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
313                  END DO
314               END DO
315            ELSE                                               ! Fixed levels
316               DO jj = 1, jpj
317                  DO ji = 1, jpi
318                     ! Euler (forward) time stepping
319                     ub(ji,jj,jk) = un(ji,jj,jk)
320                     vb(ji,jj,jk) = vn(ji,jj,jk)
321                     un(ji,jj,jk) = ua(ji,jj,jk)
322                     vn(ji,jj,jk) = va(ji,jj,jk)
323                  END DO
324               END DO
325            ENDIF
326         ELSE
327            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
328               ! caution: don't use (:,:) for this loop
329               ! it causes optimization problems on NEC in auto-tasking
330               DO jj = 1, jpj
331                  DO ji = 1, jpi
332                     zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
333                     zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
334                     ub(ji,jj,jk) = ( atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
335                       &            + atfp1 * un(ji,jj,jk) ) * zsshun1
336                     vb(ji,jj,jk) = ( atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
337                       &            + atfp1 * vn(ji,jj,jk) ) * zsshvn1
338                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
339                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
340                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1
341                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1
342                  END DO
343               END DO
344            ELSE                                               ! Fixed levels
345               DO jj = 1, jpj
346                  DO ji = 1, jpi
347                     ! Leap-frog time stepping
348                     ub(ji,jj,jk) = atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
349                     vb(ji,jj,jk) = atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
350                     un(ji,jj,jk) = ua(ji,jj,jk)
351                     vn(ji,jj,jk) = va(ji,jj,jk)
352                  END DO
353               END DO
354            ENDIF
355         ENDIF
356         !                                             ! ===============
357      END DO                                           !   End of slab
358      !                                                ! ===============
359
360      IF(ln_ctl)   THEN
361         CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, &
362            &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
363      ENDIF
364
365#if defined key_agrif
366      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
367#endif     
368
369   END SUBROUTINE dyn_nxt
370
371   !!======================================================================
372END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.