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

source: branches/2014/dev_r4822_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 4825

Last change on this file since 4825 was 4825, checked in by poddo, 10 years ago

#1404: NEMO-wave coupling WG
  • Property svn:keywords set to Id
File size: 15.8 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   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) add wave contribution to surface vertical velocity
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   ssh_nxt        : after ssh
16   !!   ssh_swp        : filter ans swap the ssh arrays
17   !!   wzv            : compute now vertical velocity
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
21   USE sbc_oce         ! surface boundary condition: ocean
22   USE domvvl          ! Variable volume
23   USE divcur          ! hor. divergence and curl      (div & cur routines)
24   USE iom             ! I/O library
25   USE restart         ! only for lrst_oce
26   USE in_out_manager  ! I/O manager
27   USE prtctl          ! Print control
28   USE phycst
29   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
30   USE lib_mpp         ! MPP library
31   USE bdy_oce
32   USE bdy_par         
33   USE bdydyn2d        ! bdy_ssh routine
34   USE diaar5, ONLY:   lk_diaar5
35   USE iom
36#if defined key_agrif
37   USE agrif_opa_update
38   USE agrif_opa_interp
39#endif
40#if defined key_asminc   
41   USE asminc          ! Assimilation increment
42#endif
43   USE wrk_nemo        ! Memory Allocation
44   USE timing          ! Timing
45   USE sbcwave,  ONLY: usd2dt, vsd2dt,wsd3d
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   ssh_nxt    ! called by step.F90
51   PUBLIC   wzv        ! called by step.F90
52   PUBLIC   ssh_swp    ! called by step.F90
53
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56#  include "vectopt_loop_substitute.h90"
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE ssh_nxt( kt )
65      !!----------------------------------------------------------------------
66      !!                ***  ROUTINE ssh_nxt  ***
67      !!                   
68      !! ** Purpose :   compute the after ssh (ssha)
69      !!
70      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
71      !!      is computed by integrating the horizontal divergence and multiply by
72      !!      by the time step.
73      !!
74      !! ** action  :   ssha    : after sea surface height
75      !!
76      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
77      !!----------------------------------------------------------------------
78      !
79      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
80      INTEGER, INTENT(in) ::   kt                      ! time step
81      !
82      INTEGER             ::   jk                      ! dummy loop indice
83      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
84      !!----------------------------------------------------------------------
85      !
86      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
87      !
88      CALL wrk_alloc( jpi, jpj, zhdiv ) 
89      !
90      IF( kt == nit000 ) THEN
91         !
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
94         IF(lwp) WRITE(numout,*) '~~~~~~~ '
95         !
96      ENDIF
97      !
98      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
99      !
100      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
101      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
102
103      !                                           !------------------------------!
104      !                                           !   After Sea Surface Height   !
105      !                                           !------------------------------!
106      zhdiv(:,:) = 0._wp
107      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
108        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
109      END DO
110      !                                                ! Sea surface elevation time stepping
111      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
112      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
113      !
114      z1_rau0 = 0.5_wp * r1_rau0
115      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
116
117#if ! defined key_dynspg_ts
118      ! These lines are not necessary with time splitting since
119      ! boundary condition on sea level is set during ts loop
120#if defined key_agrif
121      CALL agrif_ssh( kt )
122#endif
123#if defined key_bdy
124      IF (lk_bdy) THEN
125         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
126         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
127      ENDIF
128#endif
129#endif
130
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      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  sshnu, sshnv, dsshnu,dsshnv
178      !!----------------------------------------------------------------------
179     
180      IF( nn_timing == 1 )  CALL timing_start('wzv')
181      !
182      IF( kt == nit000 ) THEN
183         !
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
186         IF(lwp) WRITE(numout,*) '~~~~~ '
187         !
188         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
189         !
190      ENDIF
191      !                                           !------------------------------!
192      !                                           !     Now Vertical Velocity    !
193      !                                           !------------------------------!
194      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
195      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
196      !
197      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
198         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
199         !
200         DO jk = 1, jpkm1
201            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
202            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
203            DO jj = 2, jpjm1
204               DO ji = fs_2, fs_jpim1   ! vector opt.
205                  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) )
206               END DO
207            END DO
208         END DO
209         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
210         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
211         !                             ! Same question holds for hdivn. Perhaps just for security
212         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
213            ! computation of w
214            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
215               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
216         END DO
217         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
218         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
219      ELSE   ! z_star and linear free surface cases
220         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
221            ! computation of w
222            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
223               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
224         END DO
225      ENDIF
226!
227!     In case ln_wave and ln_sdw the surface vertical velocity is modified
228!     accounting for Sokes Drift velocity
229      IF ( ln_sdw )  THEN
230       ALLOCATE(sshnu(jpi,jpj),sshnv(jpi,jpj),dsshnu(jpi,jpj),dsshnv(jpi,jpj))
231       sshnu (:,:) = 0._wp
232       sshnv (:,:) = 0._wp
233       dsshnu(:,:) = 0._wp
234       dsshnv(:,:) = 0._wp
235       ! sshn interpolated on U-V grid
236       !---------------------------------
237       DO jj = 1 , jpjm1
238         DO ji = 1 , jpim1
239           sshnu(ji,jj) =  0.5 * ( 2. - umask(ji,jj,1) ) *               &
240                        &        ( sshn(ji  ,jj) * tmask(ji  ,jj,1)      &
241                        &        + sshn(ji+1,jj) * tmask(ji+1,jj,1) )
242           sshnv(ji,jj) =  0.5 * ( 2. - vmask(ji,jj,1) ) *               &
243                        &        ( sshn(ji,jj  ) * tmask(ji,jj  ,1)      &
244                        &        + sshn(ji,jj+1) * tmask(ji,jj+1,1) )
245         ENDDO
246       ENDDO
247       ! Compute d(ssh)/dx  and d(ssh)/dy 
248       !---------------------------------
249       DO jj = 1 , jpjm1
250         DO ji = 1 , jpim1
251           dsshnu(ji,jj) = ( sshnu(ji+1,jj) - sshnu(ji,jj) ) / e1u(ji,jj)
252           dsshnv(ji,jj) = ( sshnv(ji,jj+1) - sshnv(ji,jj) ) / e2v(ji,jj)
253         ENDDO
254       ENDDO
255       ! Compute the surface vertical velocity accounting for the Stokes Drift
256       !---------------------------------------------------------------------
257       wn(:,:,1) = wn(:,:,1) + usd2dt(:,:) * dsshnu(:,:)     &
258                 &           + vsd2dt(:,:) * dsshnv(:,:)      &
259                 &           - ( wsd3d (:,:,1) ) * tmask(:,:,1)
260      ENDIF
261!
262
263#if defined key_bdy
264      IF (lk_bdy) THEN
265         DO jk = 1, jpkm1
266            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
267         END DO
268      ENDIF
269#endif
270      !
271      !                                           !------------------------------!
272      !                                           !           outputs            !
273      !                                           !------------------------------!
274      CALL iom_put( "woce", wn )   ! vertical velocity
275      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
276         CALL wrk_alloc( jpi, jpj, z2d ) 
277         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
278         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
279         z2d(:,:) = rau0 * e12t(:,:)
280         DO jk = 1, jpk
281            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
282         END DO
283         CALL iom_put( "w_masstr" , z3d                     ) 
284         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
285         CALL wrk_dealloc( jpi, jpj, z2d  ) 
286         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
287      ENDIF
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('wzv')
290
291
292   END SUBROUTINE wzv
293
294   SUBROUTINE ssh_swp( kt )
295      !!----------------------------------------------------------------------
296      !!                    ***  ROUTINE ssh_nxt  ***
297      !!
298      !! ** Purpose :   achieve the sea surface  height time stepping by
299      !!              applying Asselin time filter and swapping the arrays
300      !!              ssha  already computed in ssh_nxt 
301      !!
302      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
303      !!              from the filter, see Leclair and Madec 2010) and swap :
304      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
305      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
306      !!                sshn = ssha
307      !!
308      !! ** action  : - sshb, sshn   : before & now sea surface height
309      !!                               ready for the next time step
310      !!
311      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
312      !!----------------------------------------------------------------------
313      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
314      !!----------------------------------------------------------------------
315      !
316      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
317      !
318      IF( kt == nit000 ) THEN
319         IF(lwp) WRITE(numout,*)
320         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
321         IF(lwp) WRITE(numout,*) '~~~~~~~ '
322      ENDIF
323
324# if defined key_dynspg_ts
325      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
326# else
327      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
328#endif
329         sshb(:,:) = sshn(:,:)                           ! before <-- now
330         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
331      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
332         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
333         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
334         sshn(:,:) = ssha(:,:)                           ! now <-- after
335      ENDIF
336      !
337      ! Update velocity at AGRIF zoom boundaries
338#if defined key_agrif
339      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
340#endif
341      !
342      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
343      !
344      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
345      !
346   END SUBROUTINE ssh_swp
347
348   !!======================================================================
349END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.