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

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2068

Last change on this file since 2068 was 2068, checked in by mlelod, 14 years ago

ticket: #663 ensuring restartability and conservation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 18.1 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   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   ssh_wzv        : after ssh & now vertical velocity
12   !!   ssh_nxt        : filter ans swap the ssh arrays
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 diaar5, ONLY :   lk_diaar5
29   USE iom
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   ssh_wzv    ! called by step.F90
35   PUBLIC   ssh_nxt    ! called by step.F90
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
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  : - 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      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
63      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
64      !!
65      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
66      !!----------------------------------------------------------------------
67      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
68      !!
69      INTEGER, INTENT(in) ::   kt   ! time step
70      !!
71      INTEGER  ::   ji, jj, jk      ! dummy loop indices
72      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
73      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
74      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
75      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
76      !!----------------------------------------------------------------------
77
78      IF( kt == nit000 ) THEN
79         !
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                  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               END DO
101            END DO
102            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
103            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
104            DO jj = 1, jpjm1
105               DO ji = 1, jpim1      ! NO Vector Opt.
106                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
107                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
108                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
109                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
110               END DO
111            END DO
112            CALL lbc_lnk( sshf_n, 'F', 1. )
113         ENDIF
114         !
115      ENDIF
116
117      !                                           !------------------------------------------!
118      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
119         !                                        !------------------------------------------!
120         DO jk = 1, jpkm1
121            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
122            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
123            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
124            !
125            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
126            fse3u (:,:,jk) = fse3u_n (:,:,jk)
127            fse3v (:,:,jk) = fse3v_n (:,:,jk)
128            fse3f (:,:,jk) = fse3f_n (:,:,jk)
129            fse3w (:,:,jk) = fse3w_n (:,:,jk)
130            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
131            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
132         END DO
133         !
134         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
135         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
136         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
137         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
138         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
139         !
140      ENDIF
141
142                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
143      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
144
145      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
146      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
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      !                                                ! Sea surface elevation time stepping
156      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
157      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
158      IF( neuler == 0 .AND. kt == nit000 ) THEN
159         z1_rau0 = 1.e0 / rau0
160         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 *   emp(:,:)                + zhdiv(:,:) )  ) * tmask(:,:,1)
161      ELSE
162         z1_rau0 = 0.5 / rau0
163         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
164      ENDIF
165
166#if defined key_obc
167      IF( Agrif_Root() ) THEN
168         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
169         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
170      ENDIF
171#endif
172      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
173      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
174         DO jj = 1, jpjm1
175            DO ji = 1, jpim1      ! NO Vector Opt.
176               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
177                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
178                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
179               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
180                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
181                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
182            END DO
183         END DO
184         ! Boundaries conditions
185         CALL lbc_lnk( sshu_a, 'U', 1. )
186         CALL lbc_lnk( sshv_a, 'V', 1. )
187      ENDIF
188      !                                           !----------------------------------------!
189      !                                           !     vertical scale factor laplacian    !
190      !                                           !----------------------------------------!
191      ! Needed for Robert-Asselin time filter and for Brown & Campana semi implicit hydrostatic presure gradient
192      DO jk = 1, jpk
193         fse3t_d(:,:,jk) =          fse3t_b(:,:,jk)   &
194            &              - 2.e0 * fse3t_n(:,:,jk)   &
195            &              +        fse3t_a(:,:,jk)
196      ENDDO
197      !                                           !------------------------------!
198      !                                           !     Now Vertical Velocity    !
199      !                                           !------------------------------!
200      z1_2dt = 1.e0 / z2dt
201      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
202         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
203         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
204            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
205            &                         * tmask(:,:,jk) * z1_2dt
206      END DO
207
208      !                                           !------------------------------!
209      !                                           !           outputs            !
210      !                                           !------------------------------!
211      CALL iom_put( "woce", wn                    )   ! vertical velocity
212      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
213      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
214      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
215         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
216         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
217         DO jk = 1, jpk
218            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
219         END DO
220         CALL iom_put( "w_masstr" , z3d                     ) 
221         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
222      ENDIF
223      !
224      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
225      !
226   END SUBROUTINE ssh_wzv
227
228
229   SUBROUTINE ssh_nxt( kt )
230      !!----------------------------------------------------------------------
231      !!                    ***  ROUTINE ssh_nxt  ***
232      !!
233      !! ** Purpose :   achieve the sea surface  height time stepping by
234      !!              applying Asselin time filter and swapping the arrays
235      !!              ssha  already computed in ssh_wzv 
236      !!
237      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
238      !!              from the filter, see Leclair and Madec 2010) and swap :
239      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
240      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
241      !!                sshn = ssha
242      !!
243      !! ** action  : - sshb, sshn   : before & now sea surface height
244      !!                               ready for the next time step
245      !!
246      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
247      !!----------------------------------------------------------------------
248      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
249      !!
250      INTEGER  ::   ji, jj   ! dummy loop indices
251      REAL(wp) ::   zec      ! temporary scalar
252      !!----------------------------------------------------------------------
253
254      IF( kt == nit000 ) THEN
255         IF(lwp) WRITE(numout,*)
256         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
257         IF(lwp) WRITE(numout,*) '~~~~~~~ '
258      ENDIF
259
260      !                       !--------------------------!
261      IF( lk_vvl ) THEN       !  Variable volume levels  !
262         !                    !--------------------------!
263         !
264         ! ssh at t-, u-, v, f-points
265         !===========================
266         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
267            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
268            sshu_n(:,:) = sshu_a(:,:)
269            sshv_n(:,:) = sshv_a(:,:)
270            DO jj = 1, jpjm1
271               DO ji = 1, jpim1      ! NO Vector Opt.
272                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
273                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
274                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
275                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
276               END DO
277            END DO
278            ! Boundaries conditions
279            CALL lbc_lnk( sshf_n, 'F', 1. )
280         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
281            zec = atfp * rdt / rau0
282            DO jj = 1, jpj
283               DO ji = 1, jpi                                ! before <-- now filtered
284                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )   &
285                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
286                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
287                  sshu_n(ji,jj) = sshu_a(ji,jj)
288                  sshv_n(ji,jj) = sshv_a(ji,jj)
289               END DO
290            END DO
291            DO jj = 1, jpjm1
292               DO ji = 1, jpim1      ! NO Vector Opt.
293                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
294                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
295                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
296                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
297               END DO
298            END DO
299            ! Boundaries conditions
300            CALL lbc_lnk( sshf_n, 'F', 1. )
301            DO jj = 1, jpjm1
302               DO ji = 1, jpim1      ! NO Vector Opt.
303                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
304                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
305                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
306                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
307                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
308                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
309               END DO
310            END DO
311            ! Boundaries conditions
312            CALL lbc_lnk( sshu_b, 'U', 1. )
313            CALL lbc_lnk( sshv_b, 'V', 1. )
314         ENDIF
315         !                    !--------------------------!
316      ELSE                    !        fixed levels      !
317         !                    !--------------------------!
318         !
319         ! ssh at t-point only
320         !====================
321         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
322            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
323            !
324         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
325            DO jj = 1, jpj
326               DO ji = 1, jpi                                ! before <-- now filtered
327                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
328                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
329               END DO
330            END DO
331         ENDIF
332         !
333      ENDIF
334      !
335      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
336      !
337   END SUBROUTINE ssh_nxt
338
339   !!======================================================================
340END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.