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 @ 1528

Last change on this file since 1528 was 1482, checked in by smasson, 15 years ago

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