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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 13.1 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/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ]
52      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*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      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
85         DO jj = 2, jpj                   ! vertical fluxes
86            DO ji = fs_2, jpi             ! vector opt.
87               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
88            END DO
89         END DO
90         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
91            DO ji = fs_2, fs_jpim1        ! vector opt.
92               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
93               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
94            END DO 
95         END DO   
96      END DO
97      !
98      ! Surface and bottom advective fluxes set to zero
99      IF ( ln_isfcav ) THEN
100         DO jj = 2, jpjm1
101            DO ji = fs_2, fs_jpim1           ! vector opt.
102               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
103               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
104               zwuw(ji,jj,jpk) = 0._wp
105               zwvw(ji,jj,jpk) = 0._wp
106            END DO
107         END DO
108      ELSE
109         DO jj = 2, jpjm1       
110            DO ji = fs_2, fs_jpim1           ! vector opt.
111               zwuw(ji,jj, 1 ) = 0._wp
112               zwvw(ji,jj, 1 ) = 0._wp
113               zwuw(ji,jj,jpk) = 0._wp
114               zwvw(ji,jj,jpk) = 0._wp
115            END DO 
116         END DO
117      END IF
118
119      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
120         DO jj = 2, jpjm1
121            DO ji = fs_2, fs_jpim1       ! vector opt.
122               !                         ! vertical momentum advective trends
123               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
124               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
125               !                         ! add the trends to the general momentum trends
126               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
127               va(ji,jj,jk) = va(ji,jj,jk) + zva
128            END DO 
129         END DO 
130      END DO
131
132      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
133         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
134         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
135         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
136         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
137      ENDIF
138      !                             ! Control print
139      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
140         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
141      !
142      CALL wrk_dealloc( jpi,jpj, zww ) 
143      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
144      !
145      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
146      !
147   END SUBROUTINE dyn_zad
148
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      ENDIF
192
193      IF( l_trddyn )   THEN         ! Save ua and va trends
194         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
195         ztrdu(:,:,:) = ua(:,:,:) 
196         ztrdv(:,:,:) = va(:,:,:) 
197      ENDIF
198     
199      IF( neuler == 0 .AND. kt == nit000 ) THEN
200          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
201      ELSE
202          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
203      ENDIF
204     
205      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
206         DO jj = 2, jpj                   
207            DO ji = fs_2, jpi             ! vector opt.
208               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
209            END DO
210         END DO
211      END DO
212
213      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
214         DO ji = fs_2, fs_jpim1           ! vector opt.
215 !!gm missing ISF boundary condition
216            zwuw(ji,jj, 1 ) = 0._wp
217            zwvw(ji,jj, 1 ) = 0._wp
218            zwuw(ji,jj,jpk) = 0._wp
219            zwvw(ji,jj,jpk) = 0._wp
220         END DO 
221      END DO
222
223! Start with before values and use sub timestepping to reach after values
224
225      zus(:,:,:,1) = ub(:,:,:)
226      zvs(:,:,:,1) = vb(:,:,:)
227
228      DO jl = 1, jnzts                   ! Start of sub timestepping loop
229
230         IF( jl == 1 ) THEN              ! Euler forward to kick things off
231           jtb = 1   ;   jtn = 1   ;   jta = 2
232           zts = z2dtzts
233         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
234           jtb = 1   ;   jtn = 2   ;   jta = 3
235           zts = 2._wp * z2dtzts
236         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
237           jtb = MOD(jtb,3) + 1
238           jtn = MOD(jtn,3) + 1
239           jta = MOD(jta,3) + 1
240         ENDIF
241
242         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
243            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
244               DO ji = fs_2, fs_jpim1        ! vector opt.
245                  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)
246                  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)
247               END DO 
248            END DO   
249         END DO
250         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
251            DO jj = 2, jpjm1
252               DO ji = fs_2, fs_jpim1       ! vector opt.
253                  !                         ! vertical momentum advective trends
254                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
255                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
256                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
257                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
258               END DO 
259            END DO 
260         END DO
261
262      END DO      ! End of sub timestepping loop
263
264      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
265      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
266         DO jj = 2, jpjm1
267            DO ji = fs_2, fs_jpim1       ! vector opt.
268               !                         ! vertical momentum advective trends
269               !                         ! add the trends to the general momentum trends
270               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
271               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
272            END DO 
273         END DO 
274      END DO
275
276      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
277         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
278         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
279         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
280         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
281      ENDIF
282      !                             ! Control print
283      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
284         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
285      !
286      CALL wrk_dealloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
287      CALL wrk_dealloc( jpi,jpj,jpk,3,   zus , zvs ) 
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
290      !
291   END SUBROUTINE dyn_zad_zts
292
293   !!======================================================================
294END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.