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 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 10 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185
  • Property svn:keywords set to Id
File size: 14.3 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 obc_par         ! open boundary cond. parameter
31   USE obc_oce
32   USE bdy_oce
33   USE bdy_par         
34   USE bdydyn2d        ! bdy_ssh routine
35   USE diaar5, ONLY:   lk_diaar5
36   USE iom
37   USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff
38   USE dynspg_ts,   ONLY: ln_bt_fw
39   USE dynspg_oce, ONLY: lk_dynspg_ts
40#if defined key_agrif
41   USE agrif_opa_update
42   USE agrif_opa_interp
43#endif
44#if defined key_asminc   
45   USE asminc          ! Assimilation increment
46#endif
47   USE wrk_nemo        ! Memory Allocation
48   USE timing          ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   ssh_nxt    ! called by step.F90
54   PUBLIC   wzv        ! called by step.F90
55   PUBLIC   ssh_swp    ! called by step.F90
56
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE ssh_nxt( kt )
68      !!----------------------------------------------------------------------
69      !!                ***  ROUTINE ssh_nxt  ***
70      !!                   
71      !! ** Purpose :   compute the after ssh (ssha)
72      !!
73      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
74      !!      is computed by integrating the horizontal divergence and multiply by
75      !!      by the time step.
76      !!
77      !! ** action  :   ssha    : after sea surface height
78      !!
79      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
80      !!----------------------------------------------------------------------
81      !
82      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
83      INTEGER, INTENT(in) ::   kt                      ! time step
84      !
85      INTEGER             ::   jk                      ! dummy loop indice
86      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
87      !!----------------------------------------------------------------------
88      !
89      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
90      !
91      CALL wrk_alloc( jpi, jpj, zhdiv ) 
92      !
93      IF( kt == nit000 ) THEN
94         !
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
97         IF(lwp) WRITE(numout,*) '~~~~~~~ '
98         !
99      ENDIF
100      !
101      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
102      !
103      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
104      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
105
106      !                                           !------------------------------!
107      !                                           !   After Sea Surface Height   !
108      !                                           !------------------------------!
109      zhdiv(:,:) = 0._wp
110      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
111        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
112      END DO
113      !                                                ! Sea surface elevation time stepping
114      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
115      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
116      z1_rau0 = 0.5_wp * r1_rau0
117      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
118
119#if defined key_agrif
120      CALL agrif_ssh( kt )
121#endif
122#if defined key_obc
123      IF( Agrif_Root() ) THEN
124         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
125         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
126      ENDIF
127#endif
128#if defined key_bdy
129      ! bg jchanut tschanges
130      ! These lines are not necessary with time splitting since
131      ! boundary condition on sea level is set during ts loop
132      IF (lk_bdy) THEN
133         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
134         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
135      ENDIF
136#endif
137      ! end jchanut tschanges
138#if defined key_asminc
139      !                                                ! Include the IAU weighted SSH increment
140      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
141         CALL ssh_asm_inc( kt )
142         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
143      ENDIF
144#endif
145
146      !                                           !------------------------------!
147      !                                           !           outputs            !
148      !                                           !------------------------------!
149      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
150      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
151      !
152      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
153      !
154      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
155      !
156      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
157      !
158   END SUBROUTINE ssh_nxt
159
160   
161   SUBROUTINE wzv( kt )
162      !!----------------------------------------------------------------------
163      !!                ***  ROUTINE wzv  ***
164      !!                   
165      !! ** Purpose :   compute the now vertical velocity
166      !!
167      !! ** Method  : - Using the incompressibility hypothesis, the vertical
168      !!      velocity is computed by integrating the horizontal divergence 
169      !!      from the bottom to the surface minus the scale factor evolution.
170      !!        The boundary conditions are w=0 at the bottom (no flux) and.
171      !!
172      !! ** action  :   wn      : now vertical velocity
173      !!
174      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
175      !!----------------------------------------------------------------------
176      !
177      INTEGER, INTENT(in) ::   kt           ! time step
178      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
179      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
180      !
181      INTEGER             ::   ji, jj, jk   ! dummy loop indices
182      REAL(wp)            ::   z1_2dt       ! local scalars
183      !!----------------------------------------------------------------------
184     
185      IF( nn_timing == 1 )  CALL timing_start('wzv')
186      !
187      IF( kt == nit000 ) THEN
188         !
189         IF(lwp) WRITE(numout,*)
190         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
191         IF(lwp) WRITE(numout,*) '~~~~~ '
192         !
193         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
194         !
195      ENDIF
196      !                                           !------------------------------!
197      !                                           !     Now Vertical Velocity    !
198      !                                           !------------------------------!
199      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
200      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
201      !
202      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
203         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
204         !
205         DO jk = 1, jpkm1
206            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
207            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
208            DO jj = 2, jpjm1
209               DO ji = fs_2, fs_jpim1   ! vector opt.
210                  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) )
211               END DO
212            END DO
213         END DO
214         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
215         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
216         !                             ! Same question holds for hdivn. Perhaps just for security
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) + zhdiv(:,:,jk)                    &
220               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
221         END DO
222         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
223         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
224      ELSE   ! z_star and linear free surface cases
225         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
226            ! computation of w
227            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
228               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
229         END DO
230      ENDIF
231
232#if defined key_bdy
233         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
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.