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

source: branches/DEV_r1784_mid_year_merge_2010/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2000

Last change on this file since 2000 was 2000, checked in by acc, 14 years ago

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