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.
dynzad.F90 in branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

Last change on this file was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.
File size: 13.8 KB
Line 
1MODULE dynzad
2   !!======================================================================
3   !!                       ***  MODULE  dynzad  ***
4   !! Ocean dynamics : vertical advection trend
5   !!======================================================================
6   !! History :  OPA  ! 1991-01  (G. Madec) Original code
7   !!            7.0  ! 1991-11  (G. Madec)
8   !!            7.5  ! 1996-01  (G. Madec) statement function for e3
9   !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
10   !!----------------------------------------------------------------------
11   
12   !!----------------------------------------------------------------------
13   !!   dyn_zad       : vertical advection momentum trend
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE sbc_oce        ! surface boundary condition: ocean
18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
21   USE in_out_manager ! I/O manager
22   USE lib_mpp        ! MPP library
23   USE prtctl         ! Print control
24   USE wrk_nemo       ! Memory Allocation
25   USE timing         ! Timing
26
27   USE yomhook, ONLY: lhook, dr_hook
28   USE parkind1, ONLY: jprb, jpim
29
30   IMPLICIT NONE
31   PRIVATE
32   
33   PUBLIC   dyn_zad       ! routine called by dynadv.F90
34   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dyn_zad ( kt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dynzad  ***
49      !!
50      !! ** Purpose :   Compute the now vertical momentum advection trend and
51      !!      add it to the general trend of momentum equation.
52      !!
53      !! ** Method  :   The now vertical advection of momentum is given by:
54      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
55      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
56      !!      Add this trend to the general trend (ua,va):
57      !!         (ua,va) = (ua,va) + w dz(u,v)
58      !!
59      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
60      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
63      !
64      INTEGER  ::   ji, jj, jk      ! dummy loop indices
65      REAL(wp) ::   zua, zva        ! temporary scalars
66      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw
67      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww
68      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
69      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
70      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
71      REAL(KIND=jprb)               :: zhook_handle
72
73      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_ZAD'
74
75      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
76
77      !!----------------------------------------------------------------------
78      !
79      IF( nn_timing == 1 )  CALL timing_start('dyn_zad')
80      !
81      CALL wrk_alloc( jpi,jpj, zww ) 
82      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 
83      !
84      IF( kt == nit000 ) THEN
85         IF(lwp)WRITE(numout,*)
86         IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
87      ENDIF
88
89      IF( l_trddyn )   THEN         ! Save ua and va trends
90         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
91         ztrdu(:,:,:) = ua(:,:,:) 
92         ztrdv(:,:,:) = va(:,:,:) 
93      ENDIF
94     
95      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
96         DO jj = 2, jpj                   ! vertical fluxes
97            DO ji = fs_2, jpi             ! vector opt.
98               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
99            END DO
100         END DO
101         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
102            DO ji = fs_2, fs_jpim1        ! vector opt.
103               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
104               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
105            END DO 
106         END DO   
107      END DO
108      !
109      ! Surface and bottom advective fluxes set to zero
110      IF ( ln_isfcav ) THEN
111         DO jj = 2, jpjm1
112            DO ji = fs_2, fs_jpim1           ! vector opt.
113               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
114               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
115               zwuw(ji,jj,jpk) = 0._wp
116               zwvw(ji,jj,jpk) = 0._wp
117            END DO
118         END DO
119      ELSE
120         DO jj = 2, jpjm1       
121            DO ji = fs_2, fs_jpim1           ! vector opt.
122               zwuw(ji,jj, 1 ) = 0._wp
123               zwvw(ji,jj, 1 ) = 0._wp
124               zwuw(ji,jj,jpk) = 0._wp
125               zwvw(ji,jj,jpk) = 0._wp
126            END DO 
127         END DO
128      END IF
129
130      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
131         DO jj = 2, jpjm1
132            DO ji = fs_2, fs_jpim1       ! vector opt.
133               !                         ! vertical momentum advective trends
134               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
135               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
136               !                         ! add the trends to the general momentum trends
137               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
138               va(ji,jj,jk) = va(ji,jj,jk) + zva
139            END DO 
140         END DO 
141      END DO
142
143      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
144         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
145         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
146         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
147         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
148      ENDIF
149      !                             ! Control print
150      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
151         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
152      !
153      CALL wrk_dealloc( jpi,jpj, zww ) 
154      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
155      !
156      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
157      !
158      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
159   END SUBROUTINE dyn_zad
160
161   SUBROUTINE dyn_zad_zts ( kt )
162      !!----------------------------------------------------------------------
163      !!                  ***  ROUTINE dynzad_zts  ***
164      !!
165      !! ** Purpose :   Compute the now vertical momentum advection trend and
166      !!      add it to the general trend of momentum equation. This version
167      !!      uses sub-timesteps for improved numerical stability with small
168      !!      vertical grid sizes. This is especially relevant when using
169      !!      embedded ice with thin surface boxes.
170      !!
171      !! ** Method  :   The now vertical advection of momentum is given by:
172      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
173      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
174      !!      Add this trend to the general trend (ua,va):
175      !!         (ua,va) = (ua,va) + w dz(u,v)
176      !!
177      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
178      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
179      !!----------------------------------------------------------------------
180      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
181      !
182      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
183      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
184      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
185      REAL(wp) ::   zua, zva        ! temporary scalars
186      REAL(wp) ::   zr_rdt          ! temporary scalar
187      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
188      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
189      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
190      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
191      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
192      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
193      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
194      REAL(KIND=jprb)               :: zhook_handle
195
196      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_ZAD_ZTS'
197
198      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
199
200      !!----------------------------------------------------------------------
201      !
202      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
203      !
204      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
205      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
206      !
207      IF( kt == nit000 ) THEN
208         IF(lwp)WRITE(numout,*)
209         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
210      ENDIF
211
212      IF( l_trddyn )   THEN         ! Save ua and va trends
213         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
214         ztrdu(:,:,:) = ua(:,:,:) 
215         ztrdv(:,:,:) = va(:,:,:) 
216      ENDIF
217     
218      IF( neuler == 0 .AND. kt == nit000 ) THEN
219          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
220      ELSE
221          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
222      ENDIF
223     
224      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
225         DO jj = 2, jpj                   
226            DO ji = fs_2, jpi             ! vector opt.
227               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
228            END DO
229         END DO
230      END DO
231      !
232      ! Surface and bottom advective fluxes set to zero
233      DO jj = 2, jpjm1       
234         DO ji = fs_2, fs_jpim1           ! vector opt.
235            zwuw(ji,jj, 1 ) = 0._wp
236            zwvw(ji,jj, 1 ) = 0._wp
237            zwuw(ji,jj,jpk) = 0._wp
238            zwvw(ji,jj,jpk) = 0._wp
239         END DO 
240      END DO
241
242! Start with before values and use sub timestepping to reach after values
243
244      zus(:,:,:,1) = ub(:,:,:)
245      zvs(:,:,:,1) = vb(:,:,:)
246
247      DO jl = 1, jnzts                   ! Start of sub timestepping loop
248
249         IF( jl == 1 ) THEN              ! Euler forward to kick things off
250           jtb = 1   ;   jtn = 1   ;   jta = 2
251           zts = z2dtzts
252         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
253           jtb = 1   ;   jtn = 2   ;   jta = 3
254           zts = 2._wp * z2dtzts
255         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
256           jtb = MOD(jtb,3) + 1
257           jtn = MOD(jtn,3) + 1
258           jta = MOD(jta,3) + 1
259         ENDIF
260
261         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
262            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
263               DO ji = fs_2, fs_jpim1        ! vector opt.
264                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)
265                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)
266               END DO 
267            END DO   
268         END DO
269         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
270            DO jj = 2, jpjm1
271               DO ji = fs_2, fs_jpim1       ! vector opt.
272                  !                         ! vertical momentum advective trends
273                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
274                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
275                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
276                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
277               END DO 
278            END DO 
279         END DO
280
281      END DO      ! End of sub timestepping loop
282
283      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
284      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
285         DO jj = 2, jpjm1
286            DO ji = fs_2, fs_jpim1       ! vector opt.
287               !                         ! vertical momentum advective trends
288               !                         ! add the trends to the general momentum trends
289               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
290               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
291            END DO 
292         END DO 
293      END DO
294
295      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
296         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
297         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
298         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
299         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
300      ENDIF
301      !                             ! Control print
302      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
303         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
304      !
305      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
306      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
307      !
308      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
309      !
310      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
311   END SUBROUTINE dyn_zad_zts
312
313   !!======================================================================
314END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.