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

Last change on this file since 782 was 782, checked in by rblod, 16 years ago

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