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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 16.4 KB
RevLine 
[3]1MODULE dynnxt
[1502]2   !!=========================================================================
[3]3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
[1502]5   !!=========================================================================
[1438]6   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
7   !!                 !  1990-10  (C. Levy, G. Madec)
8   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
9   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
10   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
11   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
12   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
14   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
[1502]16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
[2528]17   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
[1502]18   !!-------------------------------------------------------------------------
[1438]19 
[1502]20   !!-------------------------------------------------------------------------
21   !!   dyn_nxt      : obtain the next (after) horizontal velocity
22   !!-------------------------------------------------------------------------
[3]23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
[2528]25   USE sbc_oce         ! Surface boundary condition: ocean fields
26   USE phycst          ! physical constants
[1502]27   USE dynspg_oce      ! type of surface pressure gradient
28   USE dynadv          ! dynamics: vector invariant versus flux form
29   USE domvvl          ! variable volume
[367]30   USE obc_oce         ! ocean open boundary conditions
[3]31   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
[367]32   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
33   USE obcvol          ! ocean open boundary condition (obc_vol routines)
[911]34   USE bdy_oce         ! unstructured open boundary conditions
35   USE bdydta          ! unstructured open boundary conditions
36   USE bdydyn          ! unstructured open boundary conditions
[1502]37   USE in_out_manager  ! I/O manager
[3]38   USE lbclnk          ! lateral boundary condition (or mpp link)
[2715]39   USE lib_mpp         ! MPP library
[258]40   USE prtctl          ! Print control
[2528]41#if defined key_agrif
42   USE agrif_opa_interp
43#endif
[3]44
45   IMPLICIT NONE
46   PRIVATE
47
[1438]48   PUBLIC    dyn_nxt   ! routine called by step.F90
49
[592]50   !! * Substitutions
51#  include "domzgr_substitute.h90"
[2715]52   !!----------------------------------------------------------------------
[2528]53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]54   !! $Id$
[2715]55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
[3]57CONTAINS
58
59   SUBROUTINE dyn_nxt ( kt )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE dyn_nxt  ***
62      !!                   
[1502]63      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
64      !!             condition on the after velocity, achieved the time stepping
65      !!             by applying the Asselin filter on now fields and swapping
66      !!             the fields.
[3]67      !!
[1502]68      !! ** Method  : * After velocity is compute using a leap-frog scheme:
69      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
70      !!             Note that with flux form advection and variable volume layer
71      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
72      !!             velocity.
73      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
74      !!             the time stepping has already been done in dynspg module
[3]75      !!
[1502]76      !!              * Apply lateral boundary conditions on after velocity
77      !!             at the local domain boundaries through lbc_lnk call,
78      !!             at the radiative open boundaries (lk_obc=T),
79      !!             at the relaxed   open boundaries (lk_bdy=T), and
80      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
[3]81      !!
[1502]82      !!              * Apply the time filter applied and swap of the dynamics
83      !!             arrays to start the next time step:
84      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
85      !!                (un,vn) = (ua,va).
86      !!             Note that with flux form advection and variable volume layer
87      !!             (lk_vvl=T), the time filter is applied on thickness weighted
88      !!             velocity.
89      !!
90      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
91      !!               un,vn   now horizontal velocity of next time-step
[3]92      !!----------------------------------------------------------------------
[2715]93      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
94      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
95      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
96      !
[3]97      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
[2715]98      !
[3]99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[1566]100#if ! defined key_dynspg_flt
[3]101      REAL(wp) ::   z2dt         ! temporary scalar
[1566]102#endif
[2715]103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
105      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
[1502]106      !!----------------------------------------------------------------------
[3]107
[2715]108      IF( wrk_in_use(2, 1,2,3) ) THEN
109         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
110      ENDIF
111
[3]112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
[1502]118#if defined key_dynspg_flt
119      !
120      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
121      ! -------------
[3]122
[1502]123      ! Update after velocity on domain lateral boundaries      (only local domain required)
124      ! --------------------------------------------------
125      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
126      CALL lbc_lnk( va, 'V', -1. ) 
127      !
128#else
129      ! Next velocity :   Leap-frog time stepping
[1438]130      ! -------------
[1502]131      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
132      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
133      !
134      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
[1438]135         DO jk = 1, jpkm1
[1502]136            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
137            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
138         END DO
139      ELSE                                            ! applied on thickness weighted velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
142               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
[1438]143               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
[1502]144            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
145               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
[1438]146               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
[592]147         END DO
148      ENDIF
149
[1502]150
151      ! Update after velocity on domain lateral boundaries
152      ! --------------------------------------------------     
153      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
154      CALL lbc_lnk( va, 'V', -1. ) 
155      !
[3]156# if defined key_obc
[1502]157      !                                !* OBC open boundaries
[2528]158      CALL obc_dyn( kt )
[1502]159      !
[2715]160      IF( .NOT. lk_dynspg_flt ) THEN
[1502]161         ! Flather boundary condition : - Update sea surface height on each open boundary
[2715]162         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
163         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
[1502]164         !                              - Correct the barotropic velocities
[367]165         CALL obc_dyn_bt( kt )
[1438]166         !
[1502]167!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
168         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
[1438]169         !
170         IF( ln_vol_cst )   CALL obc_vol( kt )
171         !
172         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
[367]173      ENDIF
[1502]174      !
[1125]175# elif defined key_bdy 
[1502]176      !                                !* BDY open boundaries
[2715]177      IF( .NOT. lk_dynspg_flt )   CALL bdy_dyn_frs( kt )
[1438]178# endif
[1502]179      !
[392]180# if defined key_agrif
[1502]181      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
[389]182# endif
[3]183#endif
[592]184
[1438]185      ! Time filter and swap of dynamics arrays
186      ! ------------------------------------------
[1502]187      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
188         DO jk = 1, jpkm1
189            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
[1438]190            vn(:,:,jk) = va(:,:,jk)
191         END DO
[1502]192      ELSE                                             !* Leap-Frog : Asselin filter and swap
[2528]193         !                                ! =============!
194         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
195            !                             ! =============!
[1502]196            DO jk = 1, jpkm1                             
[592]197               DO jj = 1, jpj
[1502]198                  DO ji = 1, jpi   
[2528]199                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
200                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]201                     !
202                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
203                     vb(ji,jj,jk) = zvf
204                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
205                     vn(ji,jj,jk) = va(ji,jj,jk)
206                  END DO
207               END DO
208            END DO
[2528]209            !                             ! ================!
210         ELSE                             ! Variable volume !
211            !                             ! ================!
212            ! Before scale factor at t-points
213            ! -------------------------------
[1502]214            DO jk = 1, jpkm1
[2528]215               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
216                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
217                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
218            ENDDO
219            ! Add volume filter correction only at the first level of t-point scale factors
220            zec = atfp * rdt / rau0
221            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
222            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
223            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
224            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
225            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
226            !
227            IF( ln_dynadv_vec ) THEN
228               ! Before scale factor at (u/v)-points
229               ! -----------------------------------
230               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
231               DO jk = 1, jpkm1
232                  DO jj = 1, jpjm1
233                     DO ji = 1, jpim1
234                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
235                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
236                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
237                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
238                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
239                     END DO
[592]240                  END DO
241               END DO
[2528]242               ! lateral boundary conditions
243               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
244               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
245               ! Add initial scale factor to scale factor anomaly
246               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
247               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
248               ! Leap-Frog - Asselin filter and swap: applied on velocity
249               ! -----------------------------------
250               DO jk = 1, jpkm1
251                  DO jj = 1, jpj
252                     DO ji = 1, jpi
253                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
254                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
255                        !
256                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
257                        vb(ji,jj,jk) = zvf
258                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
259                        vn(ji,jj,jk) = va(ji,jj,jk)
260                     END DO
261                  END DO
262               END DO
263               !
264            ELSE
265               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
266               !-----------------------------------------------
267               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
268               DO jk = 1, jpkm1
269                  DO jj = 1, jpjm1
270                     DO ji = 1, jpim1
271                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
272                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
273                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
274                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
275                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
276                     END DO
277                  END DO
278               END DO
279               ! lateral boundary conditions
280               CALL lbc_lnk( ze3u_f, 'U', 1. )
281               CALL lbc_lnk( ze3v_f, 'V', 1. )
282               ! Add initial scale factor to scale factor anomaly
283               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
284               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
285               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
286               ! -----------------------------------             ===========================
287               DO jk = 1, jpkm1
288                  DO jj = 1, jpj
289                     DO ji = 1, jpim1
290                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
291                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
292                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
293                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
294                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
295                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
296                        !
297                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
298                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
299                        !
300                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
301                        vb(ji,jj,jk) = zvf
302                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
303                        vn(ji,jj,jk) = va(ji,jj,jk)
304                     END DO
305                  END DO
306               END DO
307               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
308               fse3v_b(:,:,:) = ze3v_f(:,:,:)
309               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
310               CALL lbc_lnk( vb, 'V', -1. )
311            ENDIF
312            !
[3]313         ENDIF
[2528]314         !
[258]315      ENDIF
[3]316
[1438]317      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
318         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
319      !
[2715]320      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
321      !
[3]322   END SUBROUTINE dyn_nxt
323
[1502]324   !!=========================================================================
[3]325END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.