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

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

nemo_v2_update_001 : CT : - add non linear free surface (variable volume) with new cpp key key_vvl

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.1 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      !! Variable volume
76      REAL(wp) ::   zsshun, zsshvn           ! temporary scalars
77      REAL(wp), DIMENSION(jpi,jpj)     ::  &
78         zsshub, zsshua, zsshvb, zsshva      ! 2D workspace
79      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
80         zfse3ub, zfse3un, zfse3ua, &        ! 3D workspace
81         zfse3vb, zfse3vn, zfse3va
82      !!----------------------------------------------------------------------
83      !!  OPA 9.0 , LOCEAN-IPSL (2005)
84      !! $Header$
85      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
86      !!----------------------------------------------------------------------
87
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
91         IF(lwp) WRITE(numout,*) '~~~~~~~'
92      ENDIF
93
94      ! Local constant initialization
95      z2dt = 2. * rdt
96      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
97
98      !! Explicit physics with thickness weighted updates
99      IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
100
101         ! Sea surface elevation time stepping
102         ! -----------------------------------
103         !
104         DO jj = 1, jpjm1
105            DO ji = 1,jpim1
106
107               ! Sea Surface Height at u-point before
108               zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
109                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
110                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * sshbb(ji+1,jj  ) )
111
112               ! Sea Surface Height at v-point before
113               zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
114                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   &
115                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * sshbb(ji  ,jj+1) )
116
117               ! Sea Surface Height at u-point after
118               zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     &
119                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
120                  &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * ssha(ji+1,jj  ) )
121
122               ! Sea Surface Height at v-point after
123               zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     &
124                  &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    &
125                  &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * ssha(ji  ,jj+1) )
126
127            END DO
128         END DO
129
130         ! Boundaries conditions
131         CALL lbc_lnk( zsshua, 'U', 1. )   ;    CALL lbc_lnk( zsshva, 'V', 1. )
132         CALL lbc_lnk( zsshub, 'U', 1. )   ;    CALL lbc_lnk( zsshvb, 'V', 1. )
133
134         ! Scale factors at before and after time step
135         ! -------------------------------------------
136         DO jk = 1, jpkm1
137            zfse3ub(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zsshub(:,:) * muu(:,:,jk) )
138            zfse3ua(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zsshua(:,:) * muu(:,:,jk) )
139            zfse3vb(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zsshvb(:,:) * muv(:,:,jk) )
140            zfse3va(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zsshva(:,:) * muv(:,:,jk) )
141         END DO
142
143         ! Asselin filtered scale factor at now time step
144         ! ----------------------------------------------
145         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
146            zfse3un(:,:,:) = fse3u(:,:,:)
147            zfse3vn(:,:,:) = fse3v(:,:,:)
148         ELSE
149            DO jk = 1, jpkm1
150               DO jj = 1, jpj
151                  DO ji = 1, jpi
152                     zsshun = atfp * ( zsshub(ji,jj) + zsshua(ji,jj) ) + atfp1 * sshu(ji,jj)
153                     zsshvn = atfp * ( zsshvb(ji,jj) + zsshva(ji,jj) ) + atfp1 * sshv(ji,jj)
154                     zfse3un(ji,jj,jk) = fsve3u(ji,jj,jk) * ( 1 + zsshun * muu(ji,jj,jk) )
155                     zfse3vn(ji,jj,jk) = fsve3v(ji,jj,jk) * ( 1 + zsshvn * muv(ji,jj,jk) )
156                  END DO
157               END DO
158            END DO
159         ENDIF
160
161         ! Thickness weighting
162         ! -------------------
163         ua(:,:,1:jpkm1) = ua(:,:,1:jpkm1) * fse3u (:,:,1:jpkm1)
164         va(:,:,1:jpkm1) = va(:,:,1:jpkm1) * fse3v (:,:,1:jpkm1)
165
166         un(:,:,1:jpkm1) = un(:,:,1:jpkm1) * fse3u (:,:,1:jpkm1)
167         vn(:,:,1:jpkm1) = vn(:,:,1:jpkm1) * fse3v (:,:,1:jpkm1)
168
169         ub(:,:,1:jpkm1) = ub(:,:,1:jpkm1) * zfse3ub(:,:,1:jpkm1)
170         vb(:,:,1:jpkm1) = vb(:,:,1:jpkm1) * zfse3vb(:,:,1:jpkm1)
171
172      ENDIF
173
174      ! Lateral boundary conditions on ( ua, va )
175      CALL lbc_lnk( ua, 'U', -1. )
176      CALL lbc_lnk( va, 'V', -1. )
177
178      !                                                ! ===============
179      DO jk = 1, jpkm1                                 ! Horizontal slab
180         !                                             ! ===============
181         ! Next velocity
182         ! -------------
183#if defined key_dynspg_flt
184         ! Leap-frog time stepping already done in dynspg.F routine
185#else
186         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
187            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
188               ! Leap-frog time stepping
189               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
190               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
191            END DO
192         END DO
193# if defined key_obc
194         !                                             ! ===============
195      END DO                                           !   End of slab
196      !                                                ! ===============
197      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
198      CALL obc_dyn( kt )
199
200      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
201         !Flather boundary condition :
202         !        - Update sea surface height on each open boundary
203         !                 sshn (= after ssh) for explicit case
204         !                 sshn_b (= after ssha_b) for time-splitting case
205         !        - Correct the barotropic velocities
206         CALL obc_dyn_bt( kt )
207
208         !Boundary conditions on sshn ( after ssh)
209         CALL lbc_lnk( sshn, 'T', 1. )
210
211         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
212            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
213         ENDIF
214
215         IF ( ln_vol_cst ) CALL obc_vol( kt )
216
217      ENDIF
218
219      !                                                ! ===============
220      DO jk = 1, jpkm1                                 ! Horizontal slab
221         !                                             ! ===============
222# endif
223# if defined key_agrif
224         !                                             ! ===============
225      END DO                                           !   End of slab
226      !                                                ! ===============
227      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
228      CALL Agrif_dyn( kt )
229      !                                                ! ===============
230      DO jk = 1, jpkm1                                 ! Horizontal slab
231         !                                             ! ===============
232# endif
233#endif
234
235         ! Time filter and swap of dynamics arrays
236         ! ------------------------------------------
237         IF( neuler == 0 .AND. kt == nit000 ) THEN
238            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
239               ! caution: don't use (:,:) for this loop
240               ! it causes optimization problems on NEC in auto-tasking
241               DO jj = 1, jpj
242                  DO ji = 1, jpi
243                     zsshun = umask(ji,jj,jk) / fse3u(ji,jj,jk)
244                     zsshvn = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
245                     ub(ji,jj,jk) = un(ji,jj,jk) * zsshun * umask(ji,jj,jk)
246                     vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn * vmask(ji,jj,jk)
247                     zsshun = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
248                     zsshvn = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
249                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun * umask(ji,jj,jk)
250                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn * vmask(ji,jj,jk)
251                  END DO
252               END DO
253            ELSE                                               ! Fixed levels
254               DO jj = 1, jpj
255                  DO ji = 1, jpi
256                     ! Euler (forward) time stepping
257                     ub(ji,jj,jk) = un(ji,jj,jk)
258                     vb(ji,jj,jk) = vn(ji,jj,jk)
259                     un(ji,jj,jk) = ua(ji,jj,jk)
260                     vn(ji,jj,jk) = va(ji,jj,jk)
261                  END DO
262               END DO
263            ENDIF
264         ELSE
265            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
266               ! caution: don't use (:,:) for this loop
267               ! it causes optimization problems on NEC in auto-tasking
268               DO jj = 1, jpj
269                  DO ji = 1, jpi
270                     zsshun = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
271                     zsshvn = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
272                     ub(ji,jj,jk) = ( atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
273                       &            + atfp1 * un(ji,jj,jk) ) * zsshun
274                     vb(ji,jj,jk) = ( atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
275                       &            + atfp1 * vn(ji,jj,jk) ) * zsshvn
276                     zsshun = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
277                     zsshvn = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
278                     un(ji,jj,jk) = ua(ji,jj,jk) * zsshun
279                     vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn
280                  END DO
281               END DO
282            ELSE                                               ! Fixed levels
283               DO jj = 1, jpj
284                  DO ji = 1, jpi
285                     ! Leap-frog time stepping
286                     ub(ji,jj,jk) = atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
287                     vb(ji,jj,jk) = atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
288                     un(ji,jj,jk) = ua(ji,jj,jk)
289                     vn(ji,jj,jk) = va(ji,jj,jk)
290                  END DO
291               END DO
292            ENDIF
293         ENDIF
294         !                                             ! ===============
295      END DO                                           !   End of slab
296      !                                                ! ===============
297
298      IF(ln_ctl)   THEN
299         CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, &
300            &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
301      ENDIF
302
303#if defined key_agrif
304      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
305#endif     
306
307   END SUBROUTINE dyn_nxt
308
309   !!======================================================================
310END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.