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

source: trunk/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 1565

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

Use appropiate e3 for divergence computation and for sea-ice, see ticket #507
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 15.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   !!   ssh_rst        : read/write ssh restart fields in the ocean restart file
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers variables
15   USE dom_oce         ! ocean space and time domain variables
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE domvvl          ! Variable volume
18   USE divcur          ! hor. divergence and curl      (div & cur routines)
19   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
20   USE iom             ! I/O library
21   USE restart         ! only for lrst_oce
22   USE in_out_manager  ! I/O manager
23   USE prtctl          ! Print control
24   USE phycst
25   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
26   USE obc_par         ! open boundary cond. parameter
27   USE obc_oce
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      INTEGER, INTENT(in) ::   kt   ! time step
69      !!
70      INTEGER  ::   ji, jj, jk      ! dummy loop indices
71      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
72      REAL(wp) ::   z2dt, zraur     ! temporary scalars
73      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
74      !!----------------------------------------------------------------------
75
76      IF( kt == nit000 ) THEN
77         IF(lwp) WRITE(numout,*)
78         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
79         IF(lwp) WRITE(numout,*) '~~~~~~~ '
80         !
81         CALL ssh_rst( nit000, 'READ' )       ! read or initialize sshb and sshn
82         !
83         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
84         !
85         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
86            DO jj = 1, jpjm1
87               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
88                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
89                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
90                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
91                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
92                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
93                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
94                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
95                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
96                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
97                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
98                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
99                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
100                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
101                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
102                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
103               END DO
104            END DO
105            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
106            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
107            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
108         ENDIF
109         !
110      ENDIF
111
112      !                                           !------------------------------!
113      !                                           !  Update Now Vertical coord.  !
114      !                                           !------------------------------!
115      IF( lk_vvl ) THEN                           ! only in vvl case)
116         !                                             ! now local depth and scale factors (stored in fse3. arrays)
117         DO jk = 1, jpkm1
118            fsdept(:,:,jk) = fsdept_n(:,:,jk)               ! depths
119            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
120            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
121            !
122            fse3t (:,:,jk) = fse3t_n (:,:,jk)               ! vertical scale factors
123            fse3u (:,:,jk) = fse3u_n (:,:,jk)
124            fse3v (:,:,jk) = fse3v_n (:,:,jk)
125            fse3f (:,:,jk) = fse3f_n (:,:,jk)
126            fse3w (:,:,jk) = fse3w_n (:,:,jk)
127            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
128            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
129         END DO
130         !                                             ! ocean depth (at u- and v-points)
131         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
132         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
133         !                                             ! masked inverse of the ocean depth (at u- and v-points)
134         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
135         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
136         !
137      ENDIF
138
139                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
140      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
141
142      ! set time step size (Euler/Leapfrog)
143      z2dt = 2. * rdt 
144      IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt
145
146      zraur = 1. / rauw
147
148      !                                           !------------------------------!
149      !                                           !   After Sea Surface Height   !
150      !                                           !------------------------------!
151      zhdiv(:,:) = 0.e0
152      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
153        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
154      END DO
155
156      !                                                ! Sea surface elevation time stepping
157      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
158
159#if defined key_obc
160# if defined key_agrif
161      IF ( Agrif_Root() ) THEN 
162# endif
163         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
164         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
165# if defined key_agrif
166      ENDIF
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, fs_jpim1   ! 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)   &   ! Caution : fmask not used
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   )                     ! vert. current
201      CALL iom_put( "ssh" , sshn )                     ! sea surface height
202      !
203   END SUBROUTINE ssh_wzv
204
205
206   SUBROUTINE ssh_nxt( kt )
207      !!----------------------------------------------------------------------
208      !!                    ***  ROUTINE ssh_nxt  ***
209      !!
210      !! ** Purpose :   achieve the sea surface  height time stepping by
211      !!              applying Asselin time filter and swapping the arrays
212      !!              ssha  already computed in ssh_wzv 
213      !!
214      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
215      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
216      !!             sshn = ssha
217      !!
218      !! ** action  : - sshb, sshn   : before & now sea surface height
219      !!                               ready for the next time step
220      !!----------------------------------------------------------------------
221      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
222      !!
223      INTEGER  ::   ji, jj               ! dummy loop indices
224      !!----------------------------------------------------------------------
225
226      IF( kt == nit000 ) THEN
227         IF(lwp) WRITE(numout,*)
228         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
229         IF(lwp) WRITE(numout,*) '~~~~~~~ '
230      ENDIF
231
232      ! Time filter and swap of the ssh
233      ! -------------------------------
234
235      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
236         !                   ! ---------------------- !
237         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
238            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
239            sshu_n(:,:) = sshu_a(:,:)
240            sshv_n(:,:) = sshv_a(:,:)
241            sshf_n(:,:) = sshf_a(:,:)
242         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
243            DO jj = 1, jpj
244               DO ji = 1, jpi                                ! before <-- now filtered
245                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
246                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
247                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
248                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
249                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
250                  sshu_n(ji,jj) = sshu_a(ji,jj)
251                  sshv_n(ji,jj) = sshv_a(ji,jj)
252                  sshf_n(ji,jj) = sshf_a(ji,jj)
253               END DO
254            END DO
255         ENDIF
256         !
257      ELSE                   ! fixed levels :   ssh at t-point only
258         !                   ! ------------ !
259         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
260            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
261            !
262         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
263            DO jj = 1, jpj
264               DO ji = 1, jpi                                ! before <-- now filtered
265                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   
266                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
267               END DO
268            END DO
269         ENDIF
270         !
271      ENDIF
272
273      !                      ! write filtered free surface arrays in restart file
274      IF( lrst_oce )   CALL ssh_rst( kt, 'WRITE' )
275      !
276      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
277      !
278   END SUBROUTINE ssh_nxt
279
280
281   SUBROUTINE ssh_rst( kt, cdrw )
282      !!---------------------------------------------------------------------
283      !!                   ***  ROUTINE ssh_rst  ***
284      !!
285      !! ** Purpose :   Read or write Sea Surface Height arrays in restart file
286      !!
287      !! ** action  : - cdrw = READ  :  sshb, sshn  read    in ocean restart file
288      !!              - cdrw = WRITE :  sshb, sshn  written in ocean restart file
289      !!----------------------------------------------------------------------
290      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
291      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
292      !!----------------------------------------------------------------------
293      !
294      IF( TRIM(cdrw) == 'READ' ) THEN
295         IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
296            CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   )
297            CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   )
298            IF( neuler == 0 )   sshb(:,:) = sshn(:,:)
299         ELSE
300            IF( nn_rstssh == 1 ) THEN
301               sshb(:,:) = 0.e0
302               sshn(:,:) = 0.e0
303            ENDIF
304         ENDIF
305      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
306         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb(:,:) )
307         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn(:,:) )
308      ENDIF
309      !
310   END SUBROUTINE ssh_rst
311
312   !!======================================================================
313END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.