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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by acc, 11 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Nearly complete port of 2011/dev_r2739_LOCEAN8_ZTC development branch into v3.5aplha base. Compiles and runs but currently unstable after 8 timesteps with ORCA2_LIM reference configuration.
  • Property svn:keywords set to Id
File size: 13.5 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 diaar5, ONLY:   lk_diaar5
34   USE iom
35   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff
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
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   ssh_nxt    ! called by step.F90
50   PUBLIC   wzv        ! called by step.F90
51   PUBLIC   ssh_swp    ! called by step.F90
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE ssh_nxt( kt )
64      !!----------------------------------------------------------------------
65      !!                ***  ROUTINE ssh_nxt  ***
66      !!                   
67      !! ** Purpose :   compute the after ssh (ssha)
68      !!
69      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
70      !!      is computed by integrating the horizontal divergence and multiply by
71      !!      by the time step.
72      !!
73      !! ** action  :   ssha    : after sea surface height
74      !!
75      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
76      !!----------------------------------------------------------------------
77      !
78      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
79      INTEGER, INTENT(in) ::   kt                      ! time step
80      !
81      INTEGER             ::   jk                      ! dummy loop indice
82      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
83      !!----------------------------------------------------------------------
84      !
85      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
86      !
87      CALL wrk_alloc( jpi, jpj, zhdiv ) 
88      !
89      IF( kt == nit000 ) THEN
90         !
91         IF(lwp) WRITE(numout,*)
92         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
93         IF(lwp) WRITE(numout,*) '~~~~~~~ '
94         !
95      ENDIF
96      !
97      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
98      !
99      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
100      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
101
102      !                                           !------------------------------!
103      !                                           !   After Sea Surface Height   !
104      !                                           !------------------------------!
105      zhdiv(:,:) = 0._wp
106      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
107        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
108      END DO
109      !                                                ! Sea surface elevation time stepping
110      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
111      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
112      z1_rau0 = 0.5 / rau0
113      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
114
115#if defined key_agrif
116      CALL agrif_ssh( kt )
117#endif
118#if defined key_obc
119      IF( Agrif_Root() ) THEN
120         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
121         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
122      ENDIF
123#endif
124#if defined key_bdy
125      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
126      CALL lbc_lnk( ssha, 'T', 1. ) 
127#endif
128#if defined key_asminc
129      !                                                ! Include the IAU weighted SSH increment
130      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
131         CALL ssh_asm_inc( kt )
132         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
133      ENDIF
134#endif
135
136      !                                           !------------------------------!
137      !                                           !           outputs            !
138      !                                           !------------------------------!
139      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
140      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
141      !
142      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
143      !
144      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
145      !
146      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
147      !
148   END SUBROUTINE ssh_nxt
149
150   
151   SUBROUTINE wzv( kt )
152      !!----------------------------------------------------------------------
153      !!                ***  ROUTINE wzv  ***
154      !!                   
155      !! ** Purpose :   compute the now vertical velocity
156      !!
157      !! ** Method  : - Using the incompressibility hypothesis, the vertical
158      !!      velocity is computed by integrating the horizontal divergence 
159      !!      from the bottom to the surface minus the scale factor evolution.
160      !!        The boundary conditions are w=0 at the bottom (no flux) and.
161      !!
162      !! ** action  :   wn      : now vertical velocity
163      !!
164      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
165      !!----------------------------------------------------------------------
166      !
167      INTEGER, INTENT(in) ::   kt           ! time step
168      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
169      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
170      !
171      INTEGER             ::   ji, jj, jk   ! dummy loop indices
172      REAL(wp)            ::   z1_2dt       ! local scalars
173      !!----------------------------------------------------------------------
174     
175      IF( nn_timing == 1 )  CALL timing_start('wzv')
176      !
177      CALL wrk_alloc( jpi, jpj, z2d ) 
178      CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 
179      !
180      IF( kt == nit000 ) THEN
181         !
182         IF(lwp) WRITE(numout,*)
183         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
184         IF(lwp) WRITE(numout,*) '~~~ '
185         !
186         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
187         !
188      ENDIF
189      !                                           !------------------------------!
190      !                                           !     Now Vertical Velocity    !
191      !                                           !------------------------------!
192      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
193      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
194      !
195      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
196         DO jk = 1, jpkm1
197            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
198            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
199            DO jj = 2, jpjm1
200               DO ji = fs_2, fs_jpim1   ! vector opt.
201                  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) )
202               END DO
203            END DO
204         END DO
205         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
206         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
207         !                             ! Same question holds for hdivn. Perhaps just for security
208         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
209            ! computation of w
210            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
211               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
212         END DO
213         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
214      ELSE   ! z_star and linear free surface cases
215         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
216            ! computation of w
217            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
218               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
219         END DO
220      ENDIF
221
222#if defined key_bdy
223         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
224#endif
225
226      !                                           !------------------------------!
227      !                                           !           outputs            !
228      !                                           !------------------------------!
229      CALL iom_put( "woce", wn )   ! vertical velocity
230      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
231         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
232         z2d(:,:) = rau0 * e12t(:,:)
233         DO jk = 1, jpk
234            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
235         END DO
236         CALL iom_put( "w_masstr" , z3d                     ) 
237         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
238      ENDIF
239      !
240      CALL wrk_dealloc( jpi, jpj, z2d  ) 
241      CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 
242      !
243      IF( nn_timing == 1 )  CALL timing_stop('wzv')
244
245
246   END SUBROUTINE wzv
247
248
249   SUBROUTINE ssh_swp( kt )
250      !!----------------------------------------------------------------------
251      !!                    ***  ROUTINE ssh_nxt  ***
252      !!
253      !! ** Purpose :   achieve the sea surface  height time stepping by
254      !!              applying Asselin time filter and swapping the arrays
255      !!              ssha  already computed in ssh_nxt 
256      !!
257      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
258      !!              from the filter, see Leclair and Madec 2010) and swap :
259      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
260      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
261      !!                sshn = ssha
262      !!
263      !! ** action  : - sshb, sshn   : before & now sea surface height
264      !!                               ready for the next time step
265      !!
266      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
267      !!----------------------------------------------------------------------
268      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
269      !!----------------------------------------------------------------------
270      !
271      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
272      !
273      IF( kt == nit000 ) THEN
274         IF(lwp) WRITE(numout,*)
275         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
276         IF(lwp) WRITE(numout,*) '~~~~~~~ '
277      ENDIF
278
279      IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
280         sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
281      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
282         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
283         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
284         sshn(:,:) = ssha(:,:)                           ! now <-- after
285      ENDIF
286      !
287      ! Update velocity at AGRIF zoom boundaries
288#if defined key_agrif
289      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
290#endif
291      !
292      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
293      !
294      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
295      !
296   END SUBROUTINE ssh_swp
297
298   !!======================================================================
299END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.