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

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

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 13.2 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   IMPLICIT NONE
28   PRIVATE
29   
30   PUBLIC   dyn_zad       ! routine called by dynadv.F90
31   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dyn_zad ( kt )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dynzad  ***
46      !!
47      !! ** Purpose :   Compute the now vertical momentum advection trend and
48      !!      add it to the general trend of momentum equation.
49      !!
50      !! ** Method  :   The now vertical advection of momentum is given by:
51      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
52      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
53      !!      Add this trend to the general trend (ua,va):
54      !!         (ua,va) = (ua,va) + w dz(u,v)
55      !!
56      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
57      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
58      !!----------------------------------------------------------------------
59      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
60      !
61      INTEGER  ::   ji, jj, jk      ! dummy loop indices
62      REAL(wp) ::   zua, zva        ! temporary scalars
63      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw
64      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww
65      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
66      !!----------------------------------------------------------------------
67      !
68      IF( nn_timing == 1 )  CALL timing_start('dyn_zad')
69      !
70      CALL wrk_alloc( jpi,jpj, zww ) 
71      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 
72      !
73      IF( kt == nit000 ) THEN
74         IF(lwp)WRITE(numout,*)
75         IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
76         IF(lwp .AND. lflush) CALL flush(numout)
77      ENDIF
78
79      IF( l_trddyn )   THEN         ! Save ua and va trends
80         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
81         ztrdu(:,:,:) = ua(:,:,:) 
82         ztrdv(:,:,:) = va(:,:,:) 
83      ENDIF
84     
85      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
86         DO jj = 2, jpj                   ! vertical fluxes
87            DO ji = fs_2, jpi             ! vector opt.
88               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
89            END DO
90         END DO
91         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
92            DO ji = fs_2, fs_jpim1        ! vector opt.
93               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
94               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
95            END DO 
96         END DO   
97      END DO
98      !
99      ! Surface and bottom advective fluxes set to zero
100      IF ( ln_isfcav ) THEN
101         DO jj = 2, jpjm1
102            DO ji = fs_2, fs_jpim1           ! vector opt.
103               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
104               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
105               zwuw(ji,jj,jpk) = 0._wp
106               zwvw(ji,jj,jpk) = 0._wp
107            END DO
108         END DO
109      ELSE
110         DO jj = 2, jpjm1       
111            DO ji = fs_2, fs_jpim1           ! vector opt.
112               zwuw(ji,jj, 1 ) = 0._wp
113               zwvw(ji,jj, 1 ) = 0._wp
114               zwuw(ji,jj,jpk) = 0._wp
115               zwvw(ji,jj,jpk) = 0._wp
116            END DO 
117         END DO
118      END IF
119
120      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
121         DO jj = 2, jpjm1
122            DO ji = fs_2, fs_jpim1       ! vector opt.
123               !                         ! vertical momentum advective trends
124               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
125               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
126               !                         ! add the trends to the general momentum trends
127               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
128               va(ji,jj,jk) = va(ji,jj,jk) + zva
129            END DO 
130         END DO 
131      END DO
132
133      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
134         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
135         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
136         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
137         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
138      ENDIF
139      !                             ! Control print
140      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
141         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
142      !
143      CALL wrk_dealloc( jpi,jpj, zww ) 
144      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
145      !
146      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
147      !
148   END SUBROUTINE dyn_zad
149
150   SUBROUTINE dyn_zad_zts ( kt )
151      !!----------------------------------------------------------------------
152      !!                  ***  ROUTINE dynzad_zts  ***
153      !!
154      !! ** Purpose :   Compute the now vertical momentum advection trend and
155      !!      add it to the general trend of momentum equation. This version
156      !!      uses sub-timesteps for improved numerical stability with small
157      !!      vertical grid sizes. This is especially relevant when using
158      !!      embedded ice with thin surface boxes.
159      !!
160      !! ** Method  :   The now vertical advection of momentum is given by:
161      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
162      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
163      !!      Add this trend to the general trend (ua,va):
164      !!         (ua,va) = (ua,va) + w dz(u,v)
165      !!
166      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
167      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
168      !!----------------------------------------------------------------------
169      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
170      !
171      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
172      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
173      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
174      REAL(wp) ::   zua, zva        ! temporary scalars
175      REAL(wp) ::   zr_rdt          ! temporary scalar
176      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
177      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
178      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
179      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
180      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
181      !!----------------------------------------------------------------------
182      !
183      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
184      !
185      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
186      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
187      !
188      IF( kt == nit000 ) THEN
189         IF(lwp)WRITE(numout,*)
190         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
191         IF(lwp .AND. lflush) CALL flush(numout)
192      ENDIF
193
194      IF( l_trddyn )   THEN         ! Save ua and va trends
195         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
196         ztrdu(:,:,:) = ua(:,:,:) 
197         ztrdv(:,:,:) = va(:,:,:) 
198      ENDIF
199     
200      IF( neuler == 0 .AND. kt == nit000 ) THEN
201          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
202      ELSE
203          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
204      ENDIF
205     
206      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
207         DO jj = 2, jpj                   
208            DO ji = fs_2, jpi             ! vector opt.
209               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
210            END DO
211         END DO
212      END DO
213      !
214      ! Surface and bottom advective fluxes set to zero
215      DO jj = 2, jpjm1       
216         DO ji = fs_2, fs_jpim1           ! vector opt.
217            zwuw(ji,jj, 1 ) = 0._wp
218            zwvw(ji,jj, 1 ) = 0._wp
219            zwuw(ji,jj,jpk) = 0._wp
220            zwvw(ji,jj,jpk) = 0._wp
221         END DO 
222      END DO
223
224! Start with before values and use sub timestepping to reach after values
225
226      zus(:,:,:,1) = ub(:,:,:)
227      zvs(:,:,:,1) = vb(:,:,:)
228
229      DO jl = 1, jnzts                   ! Start of sub timestepping loop
230
231         IF( jl == 1 ) THEN              ! Euler forward to kick things off
232           jtb = 1   ;   jtn = 1   ;   jta = 2
233           zts = z2dtzts
234         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
235           jtb = 1   ;   jtn = 2   ;   jta = 3
236           zts = 2._wp * z2dtzts
237         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
238           jtb = MOD(jtb,3) + 1
239           jtn = MOD(jtn,3) + 1
240           jta = MOD(jta,3) + 1
241         ENDIF
242
243         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
244            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
245               DO ji = fs_2, fs_jpim1        ! vector opt.
246                  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)
247                  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)
248               END DO 
249            END DO   
250         END DO
251         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
252            DO jj = 2, jpjm1
253               DO ji = fs_2, fs_jpim1       ! vector opt.
254                  !                         ! vertical momentum advective trends
255                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
256                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
257                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
258                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
259               END DO 
260            END DO 
261         END DO
262
263      END DO      ! End of sub timestepping loop
264
265      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
266      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
267         DO jj = 2, jpjm1
268            DO ji = fs_2, fs_jpim1       ! vector opt.
269               !                         ! vertical momentum advective trends
270               !                         ! add the trends to the general momentum trends
271               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
272               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
273            END DO 
274         END DO 
275      END DO
276
277      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
278         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
279         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
280         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
281         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
282      ENDIF
283      !                             ! Control print
284      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
285         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
286      !
287      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
288      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
289      !
290      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
291      !
292   END SUBROUTINE dyn_zad_zts
293
294   !!======================================================================
295END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.