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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 3432

Last change on this file since 3432 was 3432, checked in by trackstand2, 12 years ago

Merge branch 'ksection_partition'
  • Property svn:keywords set to Id
File size: 25.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   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ssh_wzv        : after ssh & now vertical velocity
14   !!   ssh_nxt        : filter ans swap the ssh arrays
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE domvvl          ! Variable volume
20   USE divcur          ! hor. divergence and curl      (div & cur routines)
21   USE iom             ! I/O library
22   USE restart         ! only for lrst_oce
23   USE in_out_manager  ! I/O manager
24   USE prtctl          ! Print control
25   USE phycst
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
27   USE lib_mpp         ! MPP library
28   USE obc_par         ! open boundary cond. parameter
29   USE obc_oce
30   USE bdy_oce
31   USE diaar5, ONLY:   lk_diaar5
32   USE iom
33   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff
34#if defined key_agrif
35   USE agrif_opa_update
36   USE agrif_opa_interp
37#endif
38#if defined key_asminc   
39   USE asminc          ! Assimilation increment
40#endif
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ssh_wzv    ! called by step.F90
46   PUBLIC   ssh_nxt    ! called by step.F90
47
48   !! * Control permutation of array indices
49#  include "oce_ftrans.h90"
50#  include "dom_oce_ftrans.h90"
51#  include "sbc_oce_ftrans.h90"
52#  include "domvvl_ftrans.h90"
53#  include "obc_oce_ftrans.h90"
54#if defined key_asminc 
55#  include "asminc_ftrans.h90"
56#endif
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE ssh_wzv( kt ) 
69      !!----------------------------------------------------------------------
70      !!                ***  ROUTINE ssh_wzv  ***
71      !!                   
72      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
73      !!              and update the now vertical coordinate (lk_vvl=T).
74      !!
75      !! ** Method  : - Using the incompressibility hypothesis, the vertical
76      !!      velocity is computed by integrating the horizontal divergence 
77      !!      from the bottom to the surface minus the scale factor evolution.
78      !!        The boundary conditions are w=0 at the bottom (no flux) and.
79      !!
80      !! ** action  :   ssha    : after sea surface height
81      !!                wn      : now vertical velocity
82      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
83      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
84      !!
85      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
86      !!----------------------------------------------------------------------
87      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
88      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
89      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
90      !! DCSE_NEMO: need additional directives for renamed module variables
91!FTRANS z3d :I :I :z
92      !
93      INTEGER, INTENT(in) ::   kt   ! time step
94      !
95      INTEGER  ::   ji, jj, jk   ! dummy loop indices
96      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
97      !!----------------------------------------------------------------------
98
99      IF( wrk_in_use(2, 1,2) ) THEN
100         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
101      ENDIF
102
103      IF( kt == nit000 ) THEN
104         !
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
107         IF(lwp) WRITE(numout,*) '~~~~~~~ '
108         !
109         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
110         !
111         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
112            DO jj = 1, jpjm1
113               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
114#if defined key_z_first
115                  zcoefu = 0.5  * umask_1(ji,jj) / ( e1u(ji,jj) * e2u(ji,jj) )
116                  zcoefv = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj) * e2v(ji,jj) )
117                  zcoeff = 0.25 * umask_1(ji,jj) * umask_1(ji,jj+1)
118#else
119                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
120                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
121                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
122#endif
123                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
124                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
125                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
126                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
127                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
128                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
129                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
130                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
131               END DO
132            END DO
133            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
134            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
135            DO jj = 1, jpjm1
136               DO ji = 1, jpim1      ! NO Vector Opt.
137#if defined key_z_first
138                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                   &
139                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
140                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
141                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
142#else
143                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
144                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
145                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
146                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
147#endif
148               END DO
149            END DO
150            CALL lbc_lnk( sshf_n, 'F', 1. )
151         ENDIF
152         !
153      ENDIF
154
155      !                                           !------------------------------------------!
156      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
157         !                                        !------------------------------------------!
158#if defined key_z_first
159            ! DCSE_NEMO: can't use implicit loop over k here because the domzgr_substitute.h90
160            ! file causes the line below to be expanded to:
161            ! gdept_1(1:jpkm1,:,:) = (gdept(1:jpkm1,:,:)*(1.+sshn(:,:)*mut(1:jpkm1,:,:)))
162            ! which contains non-conforming array expressions.
163         DO jj=1,jpj,1
164            DO ji=1,jpi,1
165               DO jk=1,jpk,1
166                  fsdept(ji,jj,jk) = fsdept_n(ji,jj,jk)   ! now local depths stored in fsdep. arrays
167               END DO
168            END DO
169         END DO
170         DO jj=1,jpj,1
171            DO ji=1,jpi,1
172               DO jk=1,jpk,1
173                  fsdepw(ji,jj,jk) = fsdepw_n(ji,jj,jk)
174               END DO
175            END DO
176         END DO
177         DO jj=1,jpj,1
178            DO ji=1,jpi,1
179               DO jk=1,jpk,1
180                  fsde3w(ji,jj,jk) = fsde3w_n(ji,jj,jk)
181               END DO
182            END DO
183         END DO
184            !
185         DO jj=1,jpj,1
186            DO ji=1,jpi,1
187               DO jk=1,jpk,1
188                  fse3t (ji,jj,jk) = fse3t_n (ji,jj,jk)   ! vertical scale factors stored in fse3. arrays
189               END DO
190            END DO
191         END DO
192         DO jj=1,jpj,1
193            DO ji=1,jpi,1
194               DO jk=1,jpk,1
195                  fse3u (ji,jj,jk) = fse3u_n (ji,jj,jk)
196               END DO
197            END DO
198         END DO
199         DO jj=1,jpj,1
200            DO ji=1,jpi,1
201               DO jk=1,jpk,1
202                  fse3v (ji,jj,jk) = fse3v_n (ji,jj,jk)
203               END DO
204            END DO
205         END DO
206         DO jj=1,jpj,1
207            DO ji=1,jpi,1
208               DO jk=1,jpk,1
209                  fse3f (ji,jj,jk) = fse3f_n (ji,jj,jk)
210               END DO
211            END DO
212         END DO
213         DO jj=1,jpj,1
214            DO ji=1,jpi,1
215               DO jk=1,jpk,1
216                  fse3w (ji,jj,jk) = fse3w_n (ji,jj,jk)
217               END DO
218            END DO
219         END DO
220
221
222         DO jj=1,jpj,1
223            DO ji=1,jpi,1
224               DO jk=1,jpk,1
225                  fse3uw(ji,jj,jk) = fse3uw_n(ji,jj,jk)
226               END DO
227            END DO
228         END DO
229
230         DO jj=1,jpj,1
231            DO ji=1,jpi,1
232               DO jk=1,jpk,1
233                  fse3vw(ji,jj,jk) = fse3vw_n(ji,jj,jk)
234               END DO
235            END DO
236         END DO
237#else
238         DO jk = 1, jpkm1
239            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
240            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
241            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
242            !
243            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
244            fse3u (:,:,jk) = fse3u_n (:,:,jk)
245            fse3v (:,:,jk) = fse3v_n (:,:,jk)
246            fse3f (:,:,jk) = fse3f_n (:,:,jk)
247            fse3w (:,:,jk) = fse3w_n (:,:,jk)
248            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
249            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
250         END DO
251#endif
252         !
253         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
254         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
255         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
256#if defined key_z_first
257         hur(:,:) = umask_1(:,:) / ( hu(:,:) + 1._wp - umask_1(:,:) )
258         hvr(:,:) = vmask_1(:,:) / ( hv(:,:) + 1._wp - vmask_1(:,:) )
259#else
260         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
261         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
262#endif
263         !
264      ENDIF
265      !
266      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
267      !
268      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
269      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
270
271      !                                           !------------------------------!
272      !                                           !   After Sea Surface Height   !
273      !                                           !------------------------------!
274      zhdiv(:,:) = 0._wp
275#if defined key_z_first
276      DO jj = 1, jpj
277         DO ji = 1, jpi
278            DO jk = 1, jpkm1                           ! Horizontal divergence of barotropic transports
279               zhdiv(ji,jj) = zhdiv(ji,jj) + fse3t(ji,jj,jk) * hdivn(ji,jj,jk)
280            END DO
281         END DO
282      END DO
283#else
284      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
285        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
286      END DO
287#endif
288      !                                                ! Sea surface elevation time stepping
289      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
290      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
291      z1_rau0 = 0.5 / rau0
292#if defined key_z_first
293      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask_1(:,:)
294#else
295      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
296#endif
297
298#if defined key_agrif
299      CALL agrif_ssh( kt )
300#endif
301#if defined key_obc
302      IF( Agrif_Root() ) THEN
303         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
304         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
305      ENDIF
306#endif
307#if defined key_bdy
308      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
309      CALL lbc_lnk( ssha, 'T', 1. ) 
310#endif
311
312      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
313      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
314         DO jj = 1, jpjm1
315            DO ji = 1, jpim1      ! NO Vector Opt.
316#if defined key_z_first
317               sshu_a(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
318                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
319                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
320               sshv_a(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
321                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
322                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
323#else
324               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
325                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
326                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
327               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
328                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
329                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
330#endif
331            END DO
332         END DO
333         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
334      ENDIF
335     
336#if defined key_asminc
337      !                                                ! Include the IAU weighted SSH increment
338      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
339         CALL ssh_asm_inc( kt )
340         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
341      ENDIF
342#endif
343
344      !                                           !------------------------------!
345      !                                           !     Now Vertical Velocity    !
346      !                                           !------------------------------!
347      z1_2dt = 1.e0 / z2dt
348#if defined key_z_first
349      DO jj = 1, jpj
350         DO ji = 1, jpi
351            DO jk = jpkm1, 1, -1                      ! integrate from the bottom the hor. divergence
352                wn(ji,jj,jk) = wn(ji,jj,jk+1)                               &
353                   &         -   fse3t_n(ji,jj,jk) * hdivn(ji,jj,jk)        &
354                   &         - ( fse3t_a(ji,jj,jk) - fse3t_b(ji,jj,jk) )    &
355                   &            * tmask(ji,jj,jk) * z1_2dt
356#if defined key_bdy
357                wn(ji,jj,jk) = wn(ji,jj,jk) * bdytmask(ji,jj)
358#endif
359            END DO
360         END DO
361      END DO
362#else
363      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
364         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
365         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
366            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
367            &                         * tmask(:,:,jk) * z1_2dt
368#if defined key_bdy
369         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
370#endif
371      END DO
372#endif
373
374      !                                           !------------------------------!
375      !                                           !           outputs            !
376      !                                           !------------------------------!
377      CALL iom_put( "woce", wn                    )   ! vertical velocity
378      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
379      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
380      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
381         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
382         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
383#if defined key_z_first
384         DO jj = 1, jpj
385            DO ji = 1, jpi
386               DO jk = 1, jpk
387                  z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj)
388               END DO
389            END DO
390         END DO
391#else
392         DO jk = 1, jpk
393            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
394         END DO
395#endif
396         CALL iom_put( "w_masstr" , z3d                     ) 
397         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
398      ENDIF
399      !
400      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
401      !
402      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
403      !
404   END SUBROUTINE ssh_wzv
405
406
407   SUBROUTINE ssh_nxt( kt )
408      !!----------------------------------------------------------------------
409      !!                    ***  ROUTINE ssh_nxt  ***
410      !!
411      !! ** Purpose :   achieve the sea surface  height time stepping by
412      !!              applying Asselin time filter and swapping the arrays
413      !!              ssha  already computed in ssh_wzv 
414      !!
415      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
416      !!              from the filter, see Leclair and Madec 2010) and swap :
417      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
418      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
419      !!                sshn = ssha
420      !!
421      !! ** action  : - sshb, sshn   : before & now sea surface height
422      !!                               ready for the next time step
423      !!
424      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
425      !!----------------------------------------------------------------------
426      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
427      !!
428      INTEGER  ::   ji, jj   ! dummy loop indices
429      REAL(wp) ::   zec      ! temporary scalar
430      !!----------------------------------------------------------------------
431
432      IF( kt == nit000 ) THEN
433         IF(lwp) WRITE(numout,*)
434         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
435         IF(lwp) WRITE(numout,*) '~~~~~~~ '
436      ENDIF
437
438      !                       !--------------------------!
439      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
440         !                    !--------------------------!
441         !
442         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
443            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
444            sshu_n(:,:) = sshu_a(:,:)
445            sshv_n(:,:) = sshv_a(:,:)
446            DO jj = 1, jpjm1                                ! ssh now at f-point
447               DO ji = 1, jpim1      ! NO Vector Opt.
448#if defined key_z_first
449                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
450                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
451                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
452                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
453#else
454                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
455                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
456                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
457                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
458#endif
459               END DO
460            END DO
461            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
462            !
463         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
464            zec = atfp * rdt / rau0
465            DO jj = 1, jpj
466               DO ji = 1, jpi                               ! before <-- now filtered
467#if defined key_z_first
468                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
469                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask_1(ji,jj)
470#else
471                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
472                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
473#endif
474                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
475                  sshu_n(ji,jj) = sshu_a(ji,jj)
476                  sshv_n(ji,jj) = sshv_a(ji,jj)
477               END DO
478            END DO
479            DO jj = 1, jpjm1                                ! ssh now at f-point
480               DO ji = 1, jpim1      ! NO Vector Opt.
481#if defined key_z_first
482                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
483                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
484                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
485                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
486#else
487                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
488                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
489                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
490                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
491#endif
492               END DO
493            END DO
494            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
495            !
496            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
497               DO ji = 1, jpim1      ! NO Vector Opt.
498#if defined key_z_first
499                  sshu_b(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
500                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
501                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
502                  sshv_b(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
503                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
504                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
505#else
506                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
507                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
508                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
509                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
510                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
511                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
512#endif
513               END DO
514            END DO
515            CALL lbc_lnk( sshu_b, 'U', 1. )
516            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
517            !
518         ENDIF
519         !                    !--------------------------!
520      ELSE                    !        fixed levels      !     (ssh at t-point only)
521         !                    !--------------------------!
522         !
523         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
524            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
525            !
526         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
527            DO jj = 1, jpj
528               DO ji = 1, jpi                               ! before <-- now filtered
529                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
530                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
531               END DO
532            END DO
533         ENDIF
534         !
535      ENDIF
536      !
537      ! Update velocity at AGRIF zoom boundaries
538#if defined key_agrif
539      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
540#endif
541      !
542      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
543      !
544   END SUBROUTINE ssh_nxt
545
546   !!======================================================================
547END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.