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.
sshwzv.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

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

Last change on this file since 4338 was 4338, checked in by acc, 10 years ago

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