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 tags/nemo_v3_2_2/NEMO/OPA_SRC/DYN – NEMO

source: tags/nemo_v3_2_2/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 4228

Last change on this file since 4228 was 2486, checked in by rblod, 14 years ago

Correct Agrif inconstency for ssh, nemo_v3_2 version, see ticket 669
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.7 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   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   ssh_wzv        : after ssh & now vertical velocity
11   !!   ssh_nxt        : filter ans swap the ssh arrays
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and tracers variables
14   USE dom_oce         ! ocean space and time domain variables
15   USE sbc_oce         ! surface boundary condition: ocean
16   USE domvvl          ! Variable volume
17   USE divcur          ! hor. divergence and curl      (div & cur routines)
18   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
19   USE iom             ! I/O library
20   USE restart         ! only for lrst_oce
21   USE in_out_manager  ! I/O manager
22   USE prtctl          ! Print control
23   USE phycst
24   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
25   USE obc_par         ! open boundary cond. parameter
26   USE obc_oce
27   USE diaar5, ONLY :   lk_diaar5
28   USE iom
29   USE agrif_opa_interp
30   USE agrif_opa_update
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ssh_wzv    ! called by step.F90
36   PUBLIC   ssh_nxt    ! called by step.F90
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
44   !! $Id$
45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE ssh_wzv( kt ) 
51      !!----------------------------------------------------------------------
52      !!                ***  ROUTINE ssh_wzv  ***
53      !!                   
54      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
55      !!              and update the now vertical coordinate (lk_vvl=T).
56      !!
57      !! ** Method  : -
58      !!
59      !!              - Using the incompressibility hypothesis, the vertical
60      !!      velocity is computed by integrating the horizontal divergence 
61      !!      from the bottom to the surface minus the scale factor evolution.
62      !!        The boundary conditions are w=0 at the bottom (no flux) and.
63      !!
64      !! ** action  :   ssha    : after sea surface height
65      !!                wn      : now vertical velocity
66      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
67      !!                          at u-, v-, f-point s
68      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
69      !!----------------------------------------------------------------------
70      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
71      !!
72      INTEGER, INTENT(in) ::   kt   ! time step
73      !!
74      INTEGER  ::   ji, jj, jk      ! dummy loop indices
75      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
76      REAL(wp) ::   z2dt, zraur     ! temporary scalars
77      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
78      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
79      !!----------------------------------------------------------------------
80
81      IF( kt == nit000 ) THEN
82         IF(lwp) WRITE(numout,*)
83         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
84         IF(lwp) WRITE(numout,*) '~~~~~~~ '
85         !
86         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
87         !
88         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
89            DO jj = 1, jpjm1
90               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
91                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
92                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
93                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
94                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
95                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
96                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
97                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
98                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
99                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
100                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
101                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
102                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
103                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
104                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
105                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
106               END DO
107            END DO
108            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
109            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
110            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
111         ENDIF
112         !
113      ENDIF
114
115      !                                           !------------------------------!
116      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
117      !                                           !------------------------------!
118         DO jk = 1, jpkm1
119            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
120            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
121            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
122            !
123            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
124            fse3u (:,:,jk) = fse3u_n (:,:,jk)
125            fse3v (:,:,jk) = fse3v_n (:,:,jk)
126            fse3f (:,:,jk) = fse3f_n (:,:,jk)
127            fse3w (:,:,jk) = fse3w_n (:,:,jk)
128            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
129            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
130         END DO
131         !                                             ! now ocean depth (at u- and v-points)
132         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
133         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
134         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
135         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
136         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
137         !
138      ENDIF
139
140                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
141      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
142
143      ! set time step size (Euler/Leapfrog)
144      z2dt = 2. * rdt 
145      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
146
147      zraur = 1. / rau0
148
149      !                                           !------------------------------!
150      !                                           !   After Sea Surface Height   !
151      !                                           !------------------------------!
152      zhdiv(:,:) = 0.e0
153      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
154        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
155      END DO
156
157      !                                                ! Sea surface elevation time stepping
158      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
159
160#if defined key_agrif
161      CALL agrif_ssh(kt)
162#endif
163#if defined key_obc
164      IF ( Agrif_Root() ) THEN
165         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
166         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
167      ENDIF
168#endif
169
170      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
171      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
172         DO jj = 1, jpjm1
173            DO ji = 1, jpim1      ! NO Vector Opt.
174               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
175                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
176                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
177               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
178                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
179                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
180               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
181                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
182                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
183            END DO
184         END DO
185         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
186         CALL lbc_lnk( sshv_a, 'V', 1. )
187         CALL lbc_lnk( sshf_a, 'F', 1. )
188      ENDIF
189
190      !                                           !------------------------------!
191      !                                           !     Now Vertical Velocity    !
192      !                                           !------------------------------!
193      !                                                ! integrate from the bottom the hor. divergence
194      DO jk = jpkm1, 1, -1
195         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
196              &                    - (  fse3t_a(:,:,jk)                   &
197              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
198      END DO
199      !
200      CALL iom_put( "woce", wn                    )   ! vertical velocity
201      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
202      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
203      IF( lk_diaar5 ) THEN
204         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
205         DO jk = 1, jpk
206            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
207         END DO
208         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
209         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
210      ENDIF
211      !
212   END SUBROUTINE ssh_wzv
213
214
215   SUBROUTINE ssh_nxt( kt )
216      !!----------------------------------------------------------------------
217      !!                    ***  ROUTINE ssh_nxt  ***
218      !!
219      !! ** Purpose :   achieve the sea surface  height time stepping by
220      !!              applying Asselin time filter and swapping the arrays
221      !!              ssha  already computed in ssh_wzv 
222      !!
223      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
224      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
225      !!             sshn = ssha
226      !!
227      !! ** action  : - sshb, sshn   : before & now sea surface height
228      !!                               ready for the next time step
229      !!----------------------------------------------------------------------
230      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
231      !!
232      INTEGER  ::   ji, jj               ! dummy loop indices
233      !!----------------------------------------------------------------------
234
235      IF( kt == nit000 ) THEN
236         IF(lwp) WRITE(numout,*)
237         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
238         IF(lwp) WRITE(numout,*) '~~~~~~~ '
239      ENDIF
240
241      ! Time filter and swap of the ssh
242      ! -------------------------------
243
244      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
245         !                   ! ---------------------- !
246         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
247            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
248            sshu_n(:,:) = sshu_a(:,:)
249            sshv_n(:,:) = sshv_a(:,:)
250            sshf_n(:,:) = sshf_a(:,:)
251         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
252            DO jj = 1, jpj
253               DO ji = 1, jpi                                ! before <-- now filtered
254                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
255                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
256                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
257                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
258                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
259                  sshu_n(ji,jj) = sshu_a(ji,jj)
260                  sshv_n(ji,jj) = sshv_a(ji,jj)
261                  sshf_n(ji,jj) = sshf_a(ji,jj)
262               END DO
263            END DO
264         ENDIF
265         !
266      ELSE                   ! fixed levels :   ssh at t-point only
267         !                   ! ------------ !
268         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
269            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
270            !
271         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
272            DO jj = 1, jpj
273               DO ji = 1, jpi                                ! before <-- now filtered
274                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   
275                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
276               END DO
277            END DO
278         ENDIF
279         !
280      ENDIF
281      !
282#if defined key_agrif
283      ! Update velocity at AGRIF zoom boundaries
284      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
285#endif
286
287      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
288      !
289   END SUBROUTINE ssh_nxt
290
291   !!======================================================================
292END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.