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

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 5014

Last change on this file since 5014 was 5014, checked in by hliu, 9 years ago

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