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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 8 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

  • Property svn:keywords set to Id
File size: 12.6 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 divhor          ! horizontal divergence
23   USE phycst          ! physical constants
24   USE bdy_oce
25   USE bdy_par         
26   USE bdydyn2d        ! bdy_ssh routine
27#if defined key_agrif
28   USE agrif_opa_interp
29#endif
30#if defined key_asminc   
31   USE asminc          ! Assimilation increment
32#endif
33   USE in_out_manager  ! I/O manager
34   USE restart         ! only for lrst_oce
35   USE prtctl          ! Print control
36   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
37   USE lib_mpp         ! MPP library
38   USE wrk_nemo        ! Memory Allocation
39   USE timing          ! Timing
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   ssh_nxt    ! called by step.F90
45   PUBLIC   wzv        ! called by step.F90
46   PUBLIC   ssh_swp    ! called by step.F90
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ssh_nxt( kt )
58      !!----------------------------------------------------------------------
59      !!                ***  ROUTINE ssh_nxt  ***
60      !!                   
61      !! ** Purpose :   compute the after ssh (ssha)
62      !!
63      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
64      !!      is computed by integrating the horizontal divergence and multiply by
65      !!      by the time step.
66      !!
67      !! ** action  :   ssha, after sea surface height
68      !!
69      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
70      !!----------------------------------------------------------------------
71      INTEGER, INTENT(in) ::   kt   ! time step
72      !
73      INTEGER  ::   jk            ! dummy loop indice
74      REAL(wp) ::   z2dt, zcoef   ! local scalars
75      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
76      !!----------------------------------------------------------------------
77      !
78      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
79      !
80      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
81      !
82      IF( kt == nit000 ) THEN
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
85         IF(lwp) WRITE(numout,*) '~~~~~~~ '
86      ENDIF
87      !
88      CALL div_hor( kt )                              ! Horizontal divergence
89      !
90      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
91      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
92
93      !                                           !------------------------------!
94      !                                           !   After Sea Surface Height   !
95      !                                           !------------------------------!
96      zhdiv(:,:) = 0._wp
97      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
98        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
99      END DO
100      !                                                ! Sea surface elevation time stepping
101      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
102      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
103      !
104      zcoef = 0.5_wp * r1_rau0
105      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
106
107#if ! defined key_dynspg_ts
108      ! These lines are not necessary with time splitting since
109      ! boundary condition on sea level is set during ts loop
110# if defined key_agrif
111      CALL agrif_ssh( kt )
112# endif
113# if defined key_bdy
114      IF( lk_bdy ) THEN
115         CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
116         CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
117      ENDIF
118# endif
119#endif
120
121#if defined key_asminc
122      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment
123         CALL ssh_asm_inc( kt )
124         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
125      ENDIF
126#endif
127      !                                           !------------------------------!
128      !                                           !           outputs            !
129      !                                           !------------------------------!
130      !
131      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
132      !
133      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
134      !
135      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
136      !
137   END SUBROUTINE ssh_nxt
138
139   
140   SUBROUTINE wzv( kt )
141      !!----------------------------------------------------------------------
142      !!                ***  ROUTINE wzv  ***
143      !!                   
144      !! ** Purpose :   compute the now vertical velocity
145      !!
146      !! ** Method  : - Using the incompressibility hypothesis, the vertical
147      !!      velocity is computed by integrating the horizontal divergence 
148      !!      from the bottom to the surface minus the scale factor evolution.
149      !!        The boundary conditions are w=0 at the bottom (no flux) and.
150      !!
151      !! ** action  :   wn      : now vertical velocity
152      !!
153      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
154      !!----------------------------------------------------------------------
155      INTEGER, INTENT(in) ::   kt   ! time step
156      !
157      INTEGER  ::   ji, jj, jk   ! dummy loop indices
158      REAL(wp) ::   z1_2dt       ! local scalars
159      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
160      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
161      !!----------------------------------------------------------------------
162      !
163      IF( nn_timing == 1 )   CALL timing_start('wzv')
164      !
165      IF( kt == nit000 ) THEN
166         IF(lwp) WRITE(numout,*)
167         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
168         IF(lwp) WRITE(numout,*) '~~~~~ '
169         !
170         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
171      ENDIF
172      !                                           !------------------------------!
173      !                                           !     Now Vertical Velocity    !
174      !                                           !------------------------------!
175      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
176      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
177      !
178      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
179         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
180         !
181         DO jk = 1, jpkm1
182            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
183            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
184            DO jj = 2, jpjm1
185               DO ji = fs_2, fs_jpim1   ! vector opt.
186                  zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
187               END DO
188            END DO
189         END DO
190         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
191         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
192         !                             ! Same question holds for hdivn. Perhaps just for security
193         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
194            ! computation of w
195            wn(:,:,jk) = wn(:,:,jk+1) - (   e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
196               &                          + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
197         END DO
198         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
199         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
200      ELSE   ! z_star and linear free surface cases
201         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
202            ! computation of w
203            wn(:,:,jk) = wn(:,:,jk+1) - (   e3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
204               &                          + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
205         END DO
206      ENDIF
207
208#if defined key_bdy
209      IF( lk_bdy ) THEN
210         DO jk = 1, jpkm1
211            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
212         END DO
213      ENDIF
214#endif
215      !
216      IF( nn_timing == 1 )  CALL timing_stop('wzv')
217      !
218   END SUBROUTINE wzv
219
220
221   SUBROUTINE ssh_swp( kt )
222      !!----------------------------------------------------------------------
223      !!                    ***  ROUTINE ssh_nxt  ***
224      !!
225      !! ** Purpose :   achieve the sea surface  height time stepping by
226      !!              applying Asselin time filter and swapping the arrays
227      !!              ssha  already computed in ssh_nxt 
228      !!
229      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
230      !!              from the filter, see Leclair and Madec 2010) and swap :
231      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
232      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
233      !!                sshn = ssha
234      !!
235      !! ** action  : - sshb, sshn   : before & now sea surface height
236      !!                               ready for the next time step
237      !!
238      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
239      !!----------------------------------------------------------------------
240      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
241      !!----------------------------------------------------------------------
242      !
243      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
244      !
245      IF( kt == nit000 ) THEN
246         IF(lwp) WRITE(numout,*)
247         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
248         IF(lwp) WRITE(numout,*) '~~~~~~~ '
249      ENDIF
250
251# if defined key_dynspg_ts
252      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
253# else
254      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
255#endif
256         sshb(:,:) = sshn(:,:)                           ! before <-- now
257         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
258         !
259      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
260         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
261         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:)    - emp(:,:)    &
262                                &                                 - rnf_b(:,:)    + rnf(:,:)    &
263                                &                                 + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
264         sshn(:,:) = ssha(:,:)                           ! now <-- after
265      ENDIF
266      !
267      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
268      !
269      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
270      !
271   END SUBROUTINE ssh_swp
272
273   !!======================================================================
274END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.