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

Last change on this file since 4328 was 4328, checked in by davestorkey, 7 years ago

Remove OBC module at NEMO 3.6. See ticket #1189.

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