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.
wzvmod.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/wzvmod.F90 @ 1438

Last change on this file since 1438 was 1438, checked in by rblod, 15 years ago

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