source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 13 months ago

The Dr Hook changes from my perl code.

File size: 20.0 KB
Line 
1MODULE dynnxt
2   !!=========================================================================
3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
5   !!=========================================================================
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.
16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
18   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
20   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
21   !!-------------------------------------------------------------------------
22 
23   !!-------------------------------------------------------------------------
24   !!   dyn_nxt      : obtain the next (after) horizontal velocity
25   !!-------------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE phycst          ! physical constants
30   USE dynspg_oce      ! type of surface pressure gradient
31   USE dynadv          ! dynamics: vector invariant versus flux form
32   USE domvvl          ! variable volume
33   USE bdy_oce         ! ocean open boundary conditions
34   USE bdydta          ! ocean open boundary conditions
35   USE bdydyn          ! ocean open boundary conditions
36   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
37   USE trd_oce         ! trends: ocean variables
38   USE trddyn          ! trend manager: dynamics
39   USE trdken          ! trend manager: kinetic energy
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! I/O manager library
43   USE lbclnk          ! lateral boundary condition (or mpp link)
44   USE lib_mpp         ! MPP library
45   USE wrk_nemo        ! Memory Allocation
46   USE prtctl          ! Print control
47   USE timing          ! Timing
48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
51
52   USE yomhook, ONLY: lhook, dr_hook
53   USE parkind1, ONLY: jprb, jpim
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC    dyn_nxt   ! routine called by step.F90
59
60   !! * Substitutions
61#  include "domzgr_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
64   !! $Id$
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE dyn_nxt ( kt )
70      !!----------------------------------------------------------------------
71      !!                  ***  ROUTINE dyn_nxt  ***
72      !!                   
73      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
74      !!             condition on the after velocity, achieved the time stepping
75      !!             by applying the Asselin filter on now fields and swapping
76      !!             the fields.
77      !!
78      !! ** Method  : * After velocity is compute using a leap-frog scheme:
79      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
80      !!             Note that with flux form advection and variable volume layer
81      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
82      !!             velocity.
83      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
84      !!             the time stepping has already been done in dynspg module
85      !!
86      !!              * Apply lateral boundary conditions on after velocity
87      !!             at the local domain boundaries through lbc_lnk call,
88      !!             at the one-way open boundaries (lk_bdy=T),
89      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
90      !!
91      !!              * Apply the time filter applied and swap of the dynamics
92      !!             arrays to start the next time step:
93      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
94      !!                (un,vn) = (ua,va).
95      !!             Note that with flux form advection and variable volume layer
96      !!             (lk_vvl=T), the time filter is applied on thickness weighted
97      !!             velocity.
98      !!
99      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
100      !!               un,vn   now horizontal velocity of next time-step
101      !!----------------------------------------------------------------------
102      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
103      !
104      INTEGER  ::   ji, jj, jk   ! dummy loop indices
105      INTEGER  ::   iku, ikv     ! local integers
106#if ! defined key_dynspg_flt
107      REAL(wp) ::   z2dt         ! temporary scalar
108#endif
109      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
110      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
111      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
112      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
113      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
114      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
115      REAL(KIND=jprb)               :: zhook_handle
116
117      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_NXT'
118
119      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
120
121      !!----------------------------------------------------------------------
122      !
123      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
124      !
125      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
126      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
127      !
128      IF( kt == nit000 ) THEN
129         IF(lwp) WRITE(numout,*)
130         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
131         IF(lwp) WRITE(numout,*) '~~~~~~~'
132      ENDIF
133
134#if defined key_dynspg_flt
135      !
136      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
137      ! -------------
138
139      ! Update after velocity on domain lateral boundaries      (only local domain required)
140      ! --------------------------------------------------
141      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
142      CALL lbc_lnk( va, 'V', -1. ) 
143      !
144#else
145
146# if defined key_dynspg_exp
147      ! Next velocity :   Leap-frog time stepping
148      ! -------------
149      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
150      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
151      !
152      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
153         DO jk = 1, jpkm1
154            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
155            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
156         END DO
157      ELSE                                            ! applied on thickness weighted velocity
158         DO jk = 1, jpkm1
159            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
160               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
161               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
162            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
163               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
164               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
165         END DO
166      ENDIF
167# endif
168
169# if defined key_dynspg_ts
170!!gm IF ( lk_dynspg_ts ) THEN ....
171      ! Ensure below that barotropic velocities match time splitting estimate
172      ! Compute actual transport and replace it with ts estimate at "after" time step
173      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
174      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
175      DO jk = 2, jpkm1
176         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
177         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
178      END DO
179      DO jk = 1, jpkm1
180         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
181         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
182      END DO
183
184      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
185         ! Remove advective velocity from "now velocities"
186         ! prior to asselin filtering     
187         ! In the forward case, this is done below after asselin filtering   
188         ! so that asselin contribution is removed at the same time
189         DO jk = 1, jpkm1
190            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
191            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
192         END DO 
193      ENDIF
194!!gm ENDIF
195# endif
196
197      ! Update after velocity on domain lateral boundaries
198      ! --------------------------------------------------     
199      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
200      CALL lbc_lnk( va, 'V', -1. ) 
201      !
202# if defined key_bdy
203      !                                !* BDY open boundaries
204      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
205      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
206
207!!$   Do we need a call to bdy_vol here??
208      !
209# endif
210      !
211# if defined key_agrif
212      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
213# endif
214#endif
215
216      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
217         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
218         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
219         !
220         !                                  ! Kinetic energy and Conversion
221         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
222         !
223         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
224            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
225            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
226            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
227            CALL iom_put( "vtrd_tot", zva )
228         ENDIF
229         !
230         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
231         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
232         !                                  !  computation of the asselin filter trends)
233      ENDIF
234
235      ! Time filter and swap of dynamics arrays
236      ! ------------------------------------------
237      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
238         DO jk = 1, jpkm1
239            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
240            vn(:,:,jk) = va(:,:,jk)
241         END DO
242         IF (lk_vvl) THEN
243            DO jk = 1, jpkm1
244               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
245               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
246               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
247            ENDDO
248         ENDIF
249      ELSE                                             !* Leap-Frog : Asselin filter and swap
250         !                                ! =============!
251         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
252            !                             ! =============!
253            DO jk = 1, jpkm1                             
254               DO jj = 1, jpj
255                  DO ji = 1, jpi   
256                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
257                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
258                     !
259                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
260                     vb(ji,jj,jk) = zvf
261                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
262                     vn(ji,jj,jk) = va(ji,jj,jk)
263                  END DO
264               END DO
265            END DO
266            !                             ! ================!
267         ELSE                             ! Variable volume !
268            !                             ! ================!
269            ! Before scale factor at t-points
270            ! (used as a now filtered scale factor until the swap)
271            ! ----------------------------------------------------
272            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
273               ! No asselin filtering on thicknesses if forward time splitting
274                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
275            ELSE
276               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
277               ! Add volume filter correction: compatibility with tracer advection scheme
278               ! => time filter + conservation correction (only at the first level)
279               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
280                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
281                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
282               ELSE                     ! if ice shelf melting
283                  DO jj = 1,jpj
284                     DO ji = 1,jpi
285                        jk = mikt(ji,jj)
286                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
287                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
288                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
289                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
290                     END DO
291                  END DO
292               END IF
293            ENDIF
294            !
295            IF( ln_dynadv_vec ) THEN
296               ! Before scale factor at (u/v)-points
297               ! -----------------------------------
298               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
299               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
300               ! Leap-Frog - Asselin filter and swap: applied on velocity
301               ! -----------------------------------
302               DO jk = 1, jpkm1
303                  DO jj = 1, jpj
304                     DO ji = 1, jpi
305                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
306                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
307                        !
308                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
309                        vb(ji,jj,jk) = zvf
310                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
311                        vn(ji,jj,jk) = va(ji,jj,jk)
312                     END DO
313                  END DO
314               END DO
315               !
316            ELSE
317               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
318               !------------------------------------------------
319               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
320               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
321               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
322               ! -----------------------------------             ===========================
323               DO jk = 1, jpkm1
324                  DO jj = 1, jpj
325                     DO ji = 1, jpi                 
326                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
327                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
328                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
329                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
330                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
331                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
332                        !
333                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
334                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
335                        !
336                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
337                        vb(ji,jj,jk) = zvf
338                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
339                        vn(ji,jj,jk) = va(ji,jj,jk)
340                     END DO
341                  END DO
342               END DO
343               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
344               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
345            ENDIF
346            !
347         ENDIF
348         !
349         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
350            ! Revert "before" velocities to time split estimate
351            ! Doing it here also means that asselin filter contribution is removed 
352            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
353            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
354            DO jk = 2, jpkm1
355               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
356               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
357            END DO
358            DO jk = 1, jpkm1
359               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
360               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
361            END DO
362         ENDIF
363         !
364      ENDIF ! neuler =/0
365      !
366      ! Set "now" and "before" barotropic velocities for next time step:
367      ! JC: Would be more clever to swap variables than to make a full vertical
368      ! integration
369      !
370      !
371      IF (lk_vvl) THEN
372         hu_b(:,:) = 0.
373         hv_b(:,:) = 0.
374         DO jk = 1, jpkm1
375            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
376            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
377         END DO
378         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
379         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
380      ENDIF
381      !
382      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
383      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
384      !
385      DO jk = 1, jpkm1
386         DO jj = 1, jpj
387            DO ji = 1, jpi
388               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
389               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
390               !
391               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
392               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
393            END DO
394         END DO
395      END DO
396      !
397      !
398      un_b(:,:) = un_b(:,:) * hur_a(:,:)
399      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
400      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
401      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
402      !
403      !
404
405      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
406         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
407         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
408         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
409      ENDIF
410      !
411      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
412         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
413      !
414      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
415      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
416      !
417      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
418      !
419      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
420   END SUBROUTINE dyn_nxt
421
422   !!=========================================================================
423END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.