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

Last change on this file since 642 was 642, checked in by opalod, 17 years ago

nemo_v2_bugfix_028 : CT : correction to be able to compile without compilation stop when using key_zco cpp key
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.8 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$
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         zfse3ub(:,:,:) = sfe3( zsshub, 'U' )   ;    zfse3ua(:,:,:) = sfe3( zsshua, 'U' )
138         zfse3vb(:,:,:) = sfe3( zsshvb, 'V' )   ;    zfse3va(:,:,:) = sfe3( zsshva, 'V' )
139
140         ! Asselin filtered scale factor at now time step
141         ! ----------------------------------------------
142         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
143            zfse3un(:,:,:) = sfe3ini( 'U' )
144            zfse3vn(:,:,:) = sfe3ini( 'V' )
145         ELSE
146            zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)
147            zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)
148            zfse3un(:,:,:) = sfe3( zsshun, 'U' )   ;    zfse3vn(:,:,:) = sfe3( zsshvn, 'V' )
149         ENDIF
150
151         ! Thickness weighting
152         ! -------------------
153         DO jk = 1, jpkm1
154            DO jj = 1, jpj
155               DO ji = 1, jpi
156                  ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
157                  va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
158
159                  un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)
160                  vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)
161
162                  ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
163                  vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
164               END DO
165            END DO
166         END DO
167
168      ENDIF
169
170      ! Lateral boundary conditions on ( ua, va )
171      CALL lbc_lnk( ua, 'U', -1. )
172      CALL lbc_lnk( va, 'V', -1. )
173
174      !                                                ! ===============
175      DO jk = 1, jpkm1                                 ! Horizontal slab
176         !                                             ! ===============
177         ! Next velocity
178         ! -------------
179#if defined key_dynspg_flt
180         ! Leap-frog time stepping already done in dynspg.F routine
181#else
182         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
183            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
184               ! Leap-frog time stepping
185               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
186               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
187            END DO
188         END DO
189# if defined key_obc
190         !                                             ! ===============
191      END DO                                           !   End of slab
192      !                                                ! ===============
193      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
194      CALL obc_dyn( kt )
195
196      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
197         !Flather boundary condition :
198         !        - Update sea surface height on each open boundary
199         !                 sshn (= after ssh) for explicit case
200         !                 sshn_b (= after ssha_b) for time-splitting case
201         !        - Correct the barotropic velocities
202         CALL obc_dyn_bt( kt )
203
204         !Boundary conditions on sshn ( after ssh)
205         CALL lbc_lnk( sshn, 'T', 1. )
206
207         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
208            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
209         ENDIF
210
211         IF ( ln_vol_cst ) CALL obc_vol( kt )
212
213      ENDIF
214
215      !                                                ! ===============
216      DO jk = 1, jpkm1                                 ! Horizontal slab
217         !                                             ! ===============
218# endif
219# if defined key_agrif
220         !                                             ! ===============
221      END DO                                           !   End of slab
222      !                                                ! ===============
223      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
224      CALL Agrif_dyn( kt )
225      !                                                ! ===============
226      DO jk = 1, jpkm1                                 ! Horizontal slab
227         !                                             ! ===============
228# endif
229#endif
230
231         ! Time filter and swap of dynamics arrays
232         ! ------------------------------------------
233         IF( neuler == 0 .AND. kt == nit000 ) THEN
234            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
235               ! caution: don't use (:,:) for this loop
236               ! it causes optimization problems on NEC in auto-tasking
237               DO jj = 1, jpj
238                  DO ji = 1, jpi
239                     zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk)
240                     zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
241                     ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
242                     vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
243                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
244                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
245                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
246                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
247                  END DO
248               END DO
249            ELSE                                               ! Fixed levels
250               DO jj = 1, jpj
251                  DO ji = 1, jpi
252                     ! Euler (forward) time stepping
253                     ub(ji,jj,jk) = un(ji,jj,jk)
254                     vb(ji,jj,jk) = vn(ji,jj,jk)
255                     un(ji,jj,jk) = ua(ji,jj,jk)
256                     vn(ji,jj,jk) = va(ji,jj,jk)
257                  END DO
258               END DO
259            ENDIF
260         ELSE
261            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
262               ! caution: don't use (:,:) for this loop
263               ! it causes optimization problems on NEC in auto-tasking
264               DO jj = 1, jpj
265                  DO ji = 1, jpi
266                     zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
267                     zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
268                     ub(ji,jj,jk) = ( atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
269                       &            + atfp1 * un(ji,jj,jk) ) * zsshun1
270                     vb(ji,jj,jk) = ( atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
271                       &            + atfp1 * vn(ji,jj,jk) ) * zsshvn1
272                     zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
273                     zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
274                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1
275                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1
276                  END DO
277               END DO
278            ELSE                                               ! Fixed levels
279               DO jj = 1, jpj
280                  DO ji = 1, jpi
281                     ! Leap-frog time stepping
282                     ub(ji,jj,jk) = atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
283                     vb(ji,jj,jk) = atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
284                     un(ji,jj,jk) = ua(ji,jj,jk)
285                     vn(ji,jj,jk) = va(ji,jj,jk)
286                  END DO
287               END DO
288            ENDIF
289         ENDIF
290         !                                             ! ===============
291      END DO                                           !   End of slab
292      !                                                ! ===============
293
294      IF(ln_ctl)   THEN
295         CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, &
296            &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
297      ENDIF
298
299#if defined key_agrif
300      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
301#endif     
302
303   END SUBROUTINE dyn_nxt
304
305   !!======================================================================
306END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.