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

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 9 years ago

added current state of wetting and drying code to test...note it does not work
  • Property svn:keywords set to Id
File size: 14.2 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.?  !  2014-09  (H. Liu) add wetting and drying 
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 wadlmt          ! Wetting/Drying flux limting
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      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
99      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
100      !
101      z1_rau0 = 0.5_wp * r1_rau0
102      !
103      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt)
104
105      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
106      !
107      !                                           !------------------------------!
108      !                                           !   After Sea Surface Height   !
109      !                                           !------------------------------!
110      zhdiv(:,:) = 0._wp
111      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
112        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
113      END DO
114      !                                                ! Sea surface elevation time stepping
115      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
116      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
117      !
118      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
119
120#if ! defined key_dynspg_ts
121      ! These lines are not necessary with time splitting since
122      ! boundary condition on sea level is set during ts loop
123#if defined key_agrif
124      CALL agrif_ssh( kt )
125#endif
126#if defined key_bdy
127      IF (lk_bdy) THEN
128         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
129         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
130      ENDIF
131#endif
132#endif
133
134#if defined key_asminc
135      !                                                ! Include the IAU weighted SSH increment
136      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
137         CALL ssh_asm_inc( kt )
138         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
139      ENDIF
140#endif
141
142      !                                           !------------------------------!
143      !                                           !           outputs            !
144      !                                           !------------------------------!
145      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
146      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
147      !
148      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
149      !
150      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
151      !
152      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
153      !
154   END SUBROUTINE ssh_nxt
155
156   
157   SUBROUTINE wzv( kt )
158      !!----------------------------------------------------------------------
159      !!                ***  ROUTINE wzv  ***
160      !!                   
161      !! ** Purpose :   compute the now vertical velocity
162      !!
163      !! ** Method  : - Using the incompressibility hypothesis, the vertical
164      !!      velocity is computed by integrating the horizontal divergence 
165      !!      from the bottom to the surface minus the scale factor evolution.
166      !!        The boundary conditions are w=0 at the bottom (no flux) and.
167      !!
168      !! ** action  :   wn      : now vertical velocity
169      !!
170      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
171      !!----------------------------------------------------------------------
172      !
173      INTEGER, INTENT(in) ::   kt           ! time step
174      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
175      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
176      !
177      INTEGER             ::   ji, jj, jk   ! dummy loop indices
178      REAL(wp)            ::   z1_2dt       ! local scalars
179      !!----------------------------------------------------------------------
180     
181      IF( nn_timing == 1 )  CALL timing_start('wzv')
182      !
183      IF( kt == nit000 ) THEN
184         !
185         IF(lwp) WRITE(numout,*)
186         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
187         IF(lwp) WRITE(numout,*) '~~~~~ '
188         !
189         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
190         !
191      ENDIF
192      !                                           !------------------------------!
193      !                                           !     Now Vertical Velocity    !
194      !                                           !------------------------------!
195      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
196      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
197      !
198      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
199         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
200         !
201         DO jk = 1, jpkm1
202            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
203            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
204            DO jj = 2, jpjm1
205               DO ji = fs_2, fs_jpim1   ! vector opt.
206                  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) )
207               END DO
208            END DO
209         END DO
210         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
211         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
212         !                             ! Same question holds for hdivn. Perhaps just for security
213         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
214            ! computation of w
215            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
216               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
217         END DO
218         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
219         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
220      ELSE   ! z_star and linear free surface cases
221         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
222            ! computation of w
223            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
224               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
225         END DO
226      ENDIF
227
228#if defined key_bdy
229      IF (lk_bdy) THEN
230         DO jk = 1, jpkm1
231            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
232         END DO
233      ENDIF
234#endif
235      !
236      !                                           !------------------------------!
237      !                                           !           outputs            !
238      !                                           !------------------------------!
239      CALL iom_put( "woce", wn )   ! vertical velocity
240      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
241         CALL wrk_alloc( jpi, jpj, z2d ) 
242         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
243         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
244         z2d(:,:) = rau0 * e12t(:,:)
245         DO jk = 1, jpk
246            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
247         END DO
248         CALL iom_put( "w_masstr" , z3d                     ) 
249         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
250         CALL wrk_dealloc( jpi, jpj, z2d  ) 
251         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
252      ENDIF
253      !
254      IF( nn_timing == 1 )  CALL timing_stop('wzv')
255
256
257   END SUBROUTINE wzv
258
259   SUBROUTINE ssh_swp( kt )
260      !!----------------------------------------------------------------------
261      !!                    ***  ROUTINE ssh_nxt  ***
262      !!
263      !! ** Purpose :   achieve the sea surface  height time stepping by
264      !!              applying Asselin time filter and swapping the arrays
265      !!              ssha  already computed in ssh_nxt 
266      !!
267      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
268      !!              from the filter, see Leclair and Madec 2010) and swap :
269      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
270      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
271      !!                sshn = ssha
272      !!
273      !! ** action  : - sshb, sshn   : before & now sea surface height
274      !!                               ready for the next time step
275      !!
276      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
277      !!----------------------------------------------------------------------
278      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
279      !!----------------------------------------------------------------------
280      !
281      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
282      !
283      IF( kt == nit000 ) THEN
284         IF(lwp) WRITE(numout,*)
285         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
286         IF(lwp) WRITE(numout,*) '~~~~~~~ '
287      ENDIF
288
289# if defined key_dynspg_ts
290      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
291# else
292      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
293#endif
294         sshb(:,:) = sshn(:,:)                           ! before <-- now
295         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
296      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
297         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
298         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
299         sshn(:,:) = ssha(:,:)                           ! now <-- after
300      ENDIF
301      !
302      ! Update velocity at AGRIF zoom boundaries
303#if defined key_agrif
304      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
305#endif
306      !
307      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
308      !
309      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
310      !
311   END SUBROUTINE ssh_swp
312
313   !!======================================================================
314END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.