source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 5944

Last change on this file since 5944 was 5944, checked in by mathiot, 5 years ago

ISF: change related to reviewers comments

  • Property svn:keywords set to Id
File size: 12.7 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      ENDIF
77
78      IF( l_trddyn )   THEN         ! Save ua and va trends
79         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
80         ztrdu(:,:,:) = ua(:,:,:) 
81         ztrdv(:,:,:) = va(:,:,:) 
82      ENDIF
83     
84      !                             ! Vertical momentum advection at uw and vw-points
85      DO jk = 2, jpk               
86         DO jj = 2, jpj                   ! vertical fluxes
87            DO ji = fs_2, jpi
88               zww(ji,jj) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk)
89            END DO
90         END DO
91         DO jj = 2, jpjm1                 ! interior vertical momentum advection at w-point
92            DO ji = fs_2, fs_jpim1
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      DO jj = 2, jpjm1                    ! First/last level advective fluxes
99         DO ji = fs_2, fs_jpim1           ! w is zero at k=1 and k=jpk
100            zwuw(ji,jj, 1 ) = 0._wp
101            zwvw(ji,jj, 1 ) = 0._wp
102            zwuw(ji,jj,jpk) = 0._wp
103            zwvw(ji,jj,jpk) = 0._wp
104         END DO
105      END DO
106
107      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
108         DO jj = 2, jpjm1
109            DO ji = fs_2, fs_jpim1       ! vector opt.
110               !                         ! vertical momentum advective trends
111               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
112               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
113               !                         ! add the trends to the general momentum trends
114               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
115               va(ji,jj,jk) = va(ji,jj,jk) + zva
116            END DO 
117         END DO 
118      END DO
119
120      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
121         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
122         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
123         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
124         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
125      ENDIF
126      !                             ! Control print
127      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
128         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
129      !
130      CALL wrk_dealloc( jpi,jpj, zww ) 
131      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
132      !
133      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
134      !
135   END SUBROUTINE dyn_zad
136
137   SUBROUTINE dyn_zad_zts ( kt )
138      !!----------------------------------------------------------------------
139      !!                  ***  ROUTINE dynzad_zts  ***
140      !!
141      !! ** Purpose :   Compute the now vertical momentum advection trend and
142      !!      add it to the general trend of momentum equation. This version
143      !!      uses sub-timesteps for improved numerical stability with small
144      !!      vertical grid sizes. This is especially relevant when using
145      !!      embedded ice with thin surface boxes.
146      !!
147      !! ** Method  :   The now vertical advection of momentum is given by:
148      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
149      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
150      !!      Add this trend to the general trend (ua,va):
151      !!         (ua,va) = (ua,va) + w dz(u,v)
152      !!
153      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
154      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
155      !!----------------------------------------------------------------------
156      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
157      !
158      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
159      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
160      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
161      REAL(wp) ::   zua, zva        ! temporary scalars
162      REAL(wp) ::   zr_rdt          ! temporary scalar
163      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
164      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
165      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
166      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
167      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
168      !!----------------------------------------------------------------------
169      !
170      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
171      !
172      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
173      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
174      !
175      IF( kt == nit000 ) THEN
176         IF(lwp)WRITE(numout,*)
177         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
178      ENDIF
179
180      IF( l_trddyn )   THEN         ! Save ua and va trends
181         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
182         ztrdu(:,:,:) = ua(:,:,:) 
183         ztrdv(:,:,:) = va(:,:,:) 
184      ENDIF
185     
186      IF( neuler == 0 .AND. kt == nit000 ) THEN
187          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
188      ELSE
189          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
190      ENDIF
191     
192      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
193         DO jj = 2, jpj                   
194            DO ji = fs_2, jpi             ! vector opt.
195               zww(ji,jj,jk) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk)
196            END DO
197         END DO
198      END DO
199      !
200      ! Surface and bottom advective fluxes set to zero
201      DO jj = 2, jpjm1       
202         DO ji = fs_2, fs_jpim1           ! vector opt.
203            zwuw(ji,jj, 1 ) = 0._wp
204            zwvw(ji,jj, 1 ) = 0._wp
205            zwuw(ji,jj,jpk) = 0._wp
206            zwvw(ji,jj,jpk) = 0._wp
207         END DO 
208      END DO
209
210! Start with before values and use sub timestepping to reach after values
211
212      zus(:,:,:,1) = ub(:,:,:)
213      zvs(:,:,:,1) = vb(:,:,:)
214
215      DO jl = 1, jnzts                   ! Start of sub timestepping loop
216
217         IF( jl == 1 ) THEN              ! Euler forward to kick things off
218           jtb = 1   ;   jtn = 1   ;   jta = 2
219           zts = z2dtzts
220         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
221           jtb = 1   ;   jtn = 2   ;   jta = 3
222           zts = 2._wp * z2dtzts
223         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
224           jtb = MOD(jtb,3) + 1
225           jtn = MOD(jtn,3) + 1
226           jta = MOD(jta,3) + 1
227         ENDIF
228
229         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
230            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
231               DO ji = fs_2, fs_jpim1        ! vector opt.
232                  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)
233                  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)
234               END DO 
235            END DO   
236         END DO
237         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
238            DO jj = 2, jpjm1
239               DO ji = fs_2, fs_jpim1       ! vector opt.
240                  !                         ! vertical momentum advective trends
241                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
242                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
243                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
244                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
245               END DO 
246            END DO 
247         END DO
248
249      END DO      ! End of sub timestepping loop
250
251      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
252      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
253         DO jj = 2, jpjm1
254            DO ji = fs_2, fs_jpim1       ! vector opt.
255               !                         ! vertical momentum advective trends
256               !                         ! add the trends to the general momentum trends
257               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
258               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
259            END DO 
260         END DO 
261      END DO
262
263      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
264         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
265         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
266         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
267         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
268      ENDIF
269      !                             ! Control print
270      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
271         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
272      !
273      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
274      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
275      !
276      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
277      !
278   END SUBROUTINE dyn_zad_zts
279
280   !!======================================================================
281END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.