source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 4327

Last change on this file since 4327 was 4327, checked in by cbricaud, 7 years ago

Correct vertical open bdy condition on wn, see ticket #1197

  • Property svn:keywords set to Id
File size: 14.4 KB
Line 
1MODULE sshwzv   
2   !!==============================================================================
3   !!                       ***  MODULE  sshwzv  ***
4   !! Ocean dynamics : sea surface height and vertical velocity
5   !!==============================================================================
6   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ssh_nxt        : after ssh
15   !!   ssh_swp        : filter ans swap the ssh arrays
16   !!   wzv            : compute now vertical velocity
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE sbc_oce         ! surface boundary condition: ocean
21   USE domvvl          ! Variable volume
22   USE divcur          ! hor. divergence and curl      (div & cur routines)
23   USE iom             ! I/O library
24   USE restart         ! only for lrst_oce
25   USE in_out_manager  ! I/O manager
26   USE prtctl          ! Print control
27   USE phycst
28   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
29   USE lib_mpp         ! MPP library
30   USE obc_par         ! open boundary cond. parameter
31   USE obc_oce
32   USE bdy_oce
33   USE bdy_par         
34   USE bdydyn2d        ! bdy_ssh routine
35   USE diaar5, ONLY:   lk_diaar5
36   USE iom
37   USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff
38   USE dynspg_ts,   ONLY: ln_bt_fw
39   USE dynspg_oce, ONLY: lk_dynspg_ts
40#if defined key_agrif
41   USE agrif_opa_update
42   USE agrif_opa_interp
43#endif
44#if defined key_asminc   
45   USE asminc          ! Assimilation increment
46#endif
47   USE wrk_nemo        ! Memory Allocation
48   USE timing          ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   ssh_nxt    ! called by step.F90
54   PUBLIC   wzv        ! called by step.F90
55   PUBLIC   ssh_swp    ! called by step.F90
56
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE ssh_nxt( kt )
68      !!----------------------------------------------------------------------
69      !!                ***  ROUTINE ssh_nxt  ***
70      !!                   
71      !! ** Purpose :   compute the after ssh (ssha)
72      !!
73      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
74      !!      is computed by integrating the horizontal divergence and multiply by
75      !!      by the time step.
76      !!
77      !! ** action  :   ssha    : after sea surface height
78      !!
79      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
80      !!----------------------------------------------------------------------
81      !
82      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
83      INTEGER, INTENT(in) ::   kt                      ! time step
84      !
85      INTEGER             ::   jk                      ! dummy loop indice
86      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
87      !!----------------------------------------------------------------------
88      !
89      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
90      !
91      CALL wrk_alloc( jpi, jpj, zhdiv ) 
92      !
93      IF( kt == nit000 ) THEN
94         !
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
97         IF(lwp) WRITE(numout,*) '~~~~~~~ '
98         !
99      ENDIF
100      !
101      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
102      !
103      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
104      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
105
106      !                                           !------------------------------!
107      !                                           !   After Sea Surface Height   !
108      !                                           !------------------------------!
109      zhdiv(:,:) = 0._wp
110      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
111        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
112      END DO
113      !                                                ! Sea surface elevation time stepping
114      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
115      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
116      z1_rau0 = 0.5_wp * r1_rau0
117      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
118
119#if defined key_agrif
120      CALL agrif_ssh( kt )
121#endif
122#if defined key_obc
123      IF( Agrif_Root() ) THEN
124         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
125         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
126      ENDIF
127#endif
128#if defined key_bdy
129      ! bg jchanut tschanges
130      ! These lines are not necessary with time splitting since
131      ! boundary condition on sea level is set during ts loop
132      IF (lk_bdy) THEN
133         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
134         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
135      ENDIF
136#endif
137      ! end jchanut tschanges
138#if defined key_asminc
139      !                                                ! Include the IAU weighted SSH increment
140      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
141         CALL ssh_asm_inc( kt )
142         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
143      ENDIF
144#endif
145
146      !                                           !------------------------------!
147      !                                           !           outputs            !
148      !                                           !------------------------------!
149      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
150      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
151      !
152      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
153      !
154      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
155      !
156      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
157      !
158   END SUBROUTINE ssh_nxt
159
160   
161   SUBROUTINE wzv( kt )
162      !!----------------------------------------------------------------------
163      !!                ***  ROUTINE wzv  ***
164      !!                   
165      !! ** Purpose :   compute the now vertical velocity
166      !!
167      !! ** Method  : - Using the incompressibility hypothesis, the vertical
168      !!      velocity is computed by integrating the horizontal divergence 
169      !!      from the bottom to the surface minus the scale factor evolution.
170      !!        The boundary conditions are w=0 at the bottom (no flux) and.
171      !!
172      !! ** action  :   wn      : now vertical velocity
173      !!
174      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
175      !!----------------------------------------------------------------------
176      !
177      INTEGER, INTENT(in) ::   kt           ! time step
178      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
179      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
180      !
181      INTEGER             ::   ji, jj, jk   ! dummy loop indices
182      REAL(wp)            ::   z1_2dt       ! local scalars
183      !!----------------------------------------------------------------------
184     
185      IF( nn_timing == 1 )  CALL timing_start('wzv')
186      !
187      IF( kt == nit000 ) THEN
188         !
189         IF(lwp) WRITE(numout,*)
190         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
191         IF(lwp) WRITE(numout,*) '~~~~~ '
192         !
193         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
194         !
195      ENDIF
196      !                                           !------------------------------!
197      !                                           !     Now Vertical Velocity    !
198      !                                           !------------------------------!
199      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
200      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
201      !
202      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
203         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
204         !
205         DO jk = 1, jpkm1
206            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
207            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
208            DO jj = 2, jpjm1
209               DO ji = fs_2, fs_jpim1   ! vector opt.
210                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
211               END DO
212            END DO
213         END DO
214         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
215         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
216         !                             ! Same question holds for hdivn. Perhaps just for security
217         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
218            ! computation of w
219            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
220               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
221         END DO
222         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
223         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
224      ELSE   ! z_star and linear free surface cases
225         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
226            ! computation of w
227            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
228               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
229         END DO
230      ENDIF
231
232#if defined key_bdy
233      IF (lk_bdy) THEN
234         DO jk = 1, jpkm1
235            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
236         END DO
237      ENDIF
238#endif
239      !
240      !                                           !------------------------------!
241      !                                           !           outputs            !
242      !                                           !------------------------------!
243      CALL iom_put( "woce", wn )   ! vertical velocity
244      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
245         CALL wrk_alloc( jpi, jpj, z2d ) 
246         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
247         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
248         z2d(:,:) = rau0 * e12t(:,:)
249         DO jk = 1, jpk
250            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
251         END DO
252         CALL iom_put( "w_masstr" , z3d                     ) 
253         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
254         CALL wrk_dealloc( jpi, jpj, z2d  ) 
255         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
256      ENDIF
257      !
258      IF( nn_timing == 1 )  CALL timing_stop('wzv')
259
260
261   END SUBROUTINE wzv
262
263   SUBROUTINE ssh_swp( kt )
264      !!----------------------------------------------------------------------
265      !!                    ***  ROUTINE ssh_nxt  ***
266      !!
267      !! ** Purpose :   achieve the sea surface  height time stepping by
268      !!              applying Asselin time filter and swapping the arrays
269      !!              ssha  already computed in ssh_nxt 
270      !!
271      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
272      !!              from the filter, see Leclair and Madec 2010) and swap :
273      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
274      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
275      !!                sshn = ssha
276      !!
277      !! ** action  : - sshb, sshn   : before & now sea surface height
278      !!                               ready for the next time step
279      !!
280      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
281      !!----------------------------------------------------------------------
282      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
283      !!----------------------------------------------------------------------
284      !
285      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
286      !
287      IF( kt == nit000 ) THEN
288         IF(lwp) WRITE(numout,*)
289         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
290         IF(lwp) WRITE(numout,*) '~~~~~~~ '
291      ENDIF
292
293# if defined key_dynspg_ts
294      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
295# else
296      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
297#endif
298         sshb(:,:) = sshn(:,:)                           ! before <-- now
299         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
300      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
301         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
302         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
303         sshn(:,:) = ssha(:,:)                           ! now <-- after
304      ENDIF
305      !
306      ! Update velocity at AGRIF zoom boundaries
307#if defined key_agrif
308      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
309#endif
310      !
311      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
312      !
313      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
314      !
315   END SUBROUTINE ssh_swp
316
317   !!======================================================================
318END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.