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

Last change on this file since 4427 was 4427, checked in by trackstand2, 7 years ago

First files changed on last FINISS work package. Stephen's work although
commited by Andy P.

  • Property svn:keywords set to Id
File size: 25.5 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#if defined key_z_first
97      INTEGER  ::   klim         ! upper bound on k loop
98#endif
99      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
100      !!----------------------------------------------------------------------
101
102      IF( wrk_in_use(2, 1,2) ) THEN
103         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
104      ENDIF
105
106      IF( kt == nit000 ) THEN
107         !
108         IF(lwp) WRITE(numout,*)
109         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
110         IF(lwp) WRITE(numout,*) '~~~~~~~ '
111         !
112#if defined key_z_first
113         DO jj=1,jpj
114            DO ji=1,jpi
115               DO jk=mbkmax(ji,jj), jpk
116                  wn(ji,jj,jk) = 0._wp
117               END DO
118            END DO
119         END DO
120#else
121         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
122#endif
123         !
124         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
125            DO jj = 1, jpjm1
126               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
127#if defined key_z_first
128                  zcoefu = 0.5  * umask_1(ji,jj) / ( e1u(ji,jj) * e2u(ji,jj) )
129                  zcoefv = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj) * e2v(ji,jj) )
130                  zcoeff = 0.25 * umask_1(ji,jj) * umask_1(ji,jj+1)
131#else
132                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
133                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
134                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
135#endif
136                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
137                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
138                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
139                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
140                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
141                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
142                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
143                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
144               END DO
145            END DO
146            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
147            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
148            DO jj = 1, jpjm1
149               DO ji = 1, jpim1      ! NO Vector Opt.
150#if defined key_z_first
151                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                   &
152                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
153                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
154                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
155#else
156                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
157                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
158                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
159                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
160#endif
161               END DO
162            END DO
163            CALL lbc_lnk( sshf_n, 'F', 1. )
164         ENDIF
165         !
166      ENDIF
167
168      !                                           !------------------------------------------!
169      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
170         !                                        !------------------------------------------!
171#if defined key_z_first
172            ! DCSE_NEMO: can't use implicit loop over k here because the domzgr_substitute.h90
173            ! file causes the line below to be expanded to:
174            ! gdept_1(1:jpkm1,:,:) = (gdept(1:jpkm1,:,:)*(1.+sshn(:,:)*mut(1:jpkm1,:,:)))
175            ! which contains non-conforming array expressions.
176         DO jj=1,jpj
177            DO ji=1,jpi
178               klim=mbkmax(ji,jj)
179               ! now local depths stored in fsdep. arrays
180               fsdept(ji,jj,1:klim) = fsdept_n(ji,jj,1:klim)
181               fsdepw(ji,jj,1:klim) = fsdepw_n(ji,jj,1:klim)
182               fsde3w(ji,jj,1:klim) = fsde3w_n(ji,jj,1:klim)
183               ! vertical scale factors stored in fse3. arrays
184               fse3t (ji,jj,1:klim) = fse3t_n (ji,jj,1:klim)
185               fse3u (ji,jj,1:klim) = fse3u_n (ji,jj,1:klim)
186               fse3v (ji,jj,1:klim) = fse3v_n (ji,jj,1:klim)
187               fse3f (ji,jj,1:klim) = fse3f_n (ji,jj,1:klim)
188               fse3w (ji,jj,1:klim) = fse3w_n (ji,jj,1:klim)
189               fse3uw(ji,jj,1:klim) = fse3uw_n(ji,jj,1:klim)
190               fse3vw(ji,jj,1:klim) = fse3vw_n(ji,jj,1:klim)
191            END DO
192         END DO
193#else
194         DO jk = 1, jpkm1
195            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
196            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
197            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
198            !
199            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
200            fse3u (:,:,jk) = fse3u_n (:,:,jk)
201            fse3v (:,:,jk) = fse3v_n (:,:,jk)
202            fse3f (:,:,jk) = fse3f_n (:,:,jk)
203            fse3w (:,:,jk) = fse3w_n (:,:,jk)
204            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
205            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
206         END DO
207#endif
208         !
209         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
210         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
211         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
212#if defined key_z_first
213         hur(:,:) = umask_1(:,:) / ( hu(:,:) + 1._wp - umask_1(:,:) )
214         hvr(:,:) = vmask_1(:,:) / ( hv(:,:) + 1._wp - vmask_1(:,:) )
215#else
216         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
217         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
218#endif
219         !
220      ENDIF
221      !
222      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
223      !
224      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
225      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
226
227      WRITE(*,*)'ARPDBG, ssh_wzv: sum WWW of hdivn=',SUM(hdivn),' at step=',kt
228      WRITE(*,*)'ARPDBG, ssh_wzv: sum WWW of fse3t=',SUM(fse3t(:,:,:)),' at step=',kt
229
230      !                                           !------------------------------!
231      !                                           !   After Sea Surface Height   !
232      !                                           !------------------------------!
233      zhdiv(:,:) = 0._wp
234#if defined key_z_first
235      DO jj = 1, jpj
236         DO ji = 1, jpi
237            DO jk = 1, mbkmax(ji,jj)-1                 ! Horizontal divergence of barotropic transports
238               zhdiv(ji,jj) = zhdiv(ji,jj) + fse3t(ji,jj,jk) * hdivn(ji,jj,jk)
239            END DO
240         END DO
241      END DO
242#else
243      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
244        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
245      END DO
246#endif
247      WRITE(*,*)'ARPDBG, ssh_wzv: sum XXX of zhdiv=',SUM(zhdiv),' at step=',kt
248      WRITE(*,*)'ARPDBG, ssh_wzv: sum XXX of ssha=',SUM(ssha),' at step=',kt
249      !                                                ! Sea surface elevation time stepping
250      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
251      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
252      z1_rau0 = 0.5 / rau0
253#if defined key_z_first
254      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask_1(:,:)
255#else
256      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
257#endif
258      WRITE(*,*)'ARPDBG, ssh_wzv: sum YYY of sshb=',SUM(sshb),' at step=',kt
259      WRITE(*,*)'ARPDBG, ssh_wzv: sum YYY of ssha=',SUM(ssha),' at step=',kt
260#if defined key_agrif
261      CALL agrif_ssh( kt )
262#endif
263#if defined key_obc
264      IF( Agrif_Root() ) THEN
265         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
266         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
267      ENDIF
268#endif
269#if defined key_bdy
270      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
271      CALL lbc_lnk( ssha, 'T', 1. ) 
272#endif
273
274      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
275      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
276         DO jj = 1, jpjm1
277            DO ji = 1, jpim1      ! NO Vector Opt.
278#if defined key_z_first
279               sshu_a(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
280                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
281                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
282               sshv_a(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
283                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
284                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
285#else
286               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
287                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
288                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
289               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
290                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
291                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
292#endif
293            END DO
294         END DO
295         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
296      ENDIF
297     
298#if defined key_asminc
299      !                                                ! Include the IAU weighted SSH increment
300      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
301         CALL ssh_asm_inc( kt )
302         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
303      ENDIF
304#endif
305
306      !                                           !------------------------------!
307      !                                           !     Now Vertical Velocity    !
308      !                                           !------------------------------!
309      z1_2dt = 1.e0 / z2dt
310#if defined key_z_first
311      DO jj = 1, jpj
312         DO ji = 1, jpi
313            DO jk = mbkmax(ji,jj)-1, 1, -1                      ! integrate from the bottom the hor. divergence
314                wn(ji,jj,jk) = wn(ji,jj,jk+1)                               &
315                   &         -   fse3t_n(ji,jj,jk) * hdivn(ji,jj,jk)        &
316                   &         - ( fse3t_a(ji,jj,jk) - fse3t_b(ji,jj,jk) )    &
317                   &            * tmask(ji,jj,jk) * z1_2dt
318#if defined key_bdy
319                wn(ji,jj,jk) = wn(ji,jj,jk) * bdytmask(ji,jj)
320#endif
321            END DO
322         END DO
323      END DO
324#else
325      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
326         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
327         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
328            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
329            &                         * tmask(:,:,jk) * z1_2dt
330#if defined key_bdy
331         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
332#endif
333      END DO
334#endif
335
336      !                                           !------------------------------!
337      !                                           !           outputs            !
338      !                                           !------------------------------!
339      CALL iom_put( "woce", wn                    )   ! vertical velocity
340      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
341      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
342      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
343         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
344         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
345#if defined key_z_first
346         DO jj = 1, jpj
347            DO ji = 1, jpi
348               DO jk = 1, mbkmax(ji,jj)
349                  z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj)
350               END DO
351            END DO
352         END DO
353#else
354         DO jk = 1, jpk
355            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
356         END DO
357#endif
358         CALL iom_put( "w_masstr" , z3d                     ) 
359         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
360      ENDIF
361      !
362      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
363      !
364      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
365      !
366   END SUBROUTINE ssh_wzv
367
368
369   SUBROUTINE ssh_nxt( kt )
370      !!----------------------------------------------------------------------
371      !!                    ***  ROUTINE ssh_nxt  ***
372      !!
373      !! ** Purpose :   achieve the sea surface  height time stepping by
374      !!              applying Asselin time filter and swapping the arrays
375      !!              ssha  already computed in ssh_wzv 
376      !!
377      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
378      !!              from the filter, see Leclair and Madec 2010) and swap :
379      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
380      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
381      !!                sshn = ssha
382      !!
383      !! ** action  : - sshb, sshn   : before & now sea surface height
384      !!                               ready for the next time step
385      !!
386      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
387      !!----------------------------------------------------------------------
388      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
389      !!
390      INTEGER  ::   ji, jj   ! dummy loop indices
391      REAL(wp) ::   zec      ! temporary scalar
392      !!----------------------------------------------------------------------
393
394      IF( kt == nit000 ) THEN
395         IF(lwp) WRITE(numout,*)
396         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
397         IF(lwp) WRITE(numout,*) '~~~~~~~ '
398      ENDIF
399
400      !                       !--------------------------!
401      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
402         !                    !--------------------------!
403         !
404         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
405            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
406            sshu_n(:,:) = sshu_a(:,:)
407            sshv_n(:,:) = sshv_a(:,:)
408            DO jj = 1, jpjm1                                ! ssh now at f-point
409               DO ji = 1, jpim1      ! NO Vector Opt.
410#if defined key_z_first
411                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
412                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
413                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
414                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
415#else
416                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
417                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
418                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
419                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
420#endif
421               END DO
422            END DO
423            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
424            !
425         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
426            WRITE(*,*) 'ARPDBG: ssh_nxt: SUM of sshb = ',SUM(sshb),' at step=',kt
427            WRITE(*,*) 'ARPDBG: ssh_nxt: SUM of sshn = ',SUM(sshn),' at step=',kt
428            zec = atfp * rdt / rau0
429            DO jj = 1, jpj
430               DO ji = 1, jpi                               ! before <-- now filtered
431#if defined key_z_first
432                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
433                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask_1(ji,jj)
434#else
435                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
436                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
437#endif
438                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
439                  sshu_n(ji,jj) = sshu_a(ji,jj)
440                  sshv_n(ji,jj) = sshv_a(ji,jj)
441               END DO
442            END DO
443            DO jj = 1, jpjm1                                ! ssh now at f-point
444               DO ji = 1, jpim1      ! NO Vector Opt.
445#if defined key_z_first
446                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
447                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
448                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
449                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
450#else
451                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
452                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
453                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
454                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
455#endif
456               END DO
457            END DO
458            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
459            !
460            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
461               DO ji = 1, jpim1      ! NO Vector Opt.
462#if defined key_z_first
463                  sshu_b(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
464                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
465                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
466                  sshv_b(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
467                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
468                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
469#else
470                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
471                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
472                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
473                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
474                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
475                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
476#endif
477               END DO
478            END DO
479            CALL lbc_lnk( sshu_b, 'U', 1. )
480            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
481            !
482         ENDIF
483         !                    !--------------------------!
484      ELSE                    !        fixed levels      !     (ssh at t-point only)
485         !                    !--------------------------!
486         !
487         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
488            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
489            !
490         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
491            DO jj = 1, jpj
492               DO ji = 1, jpi                               ! before <-- now filtered
493                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
494                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
495               END DO
496            END DO
497         ENDIF
498         !
499      ENDIF
500      !
501      ! Update velocity at AGRIF zoom boundaries
502#if defined key_agrif
503      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
504#endif
505      !
506      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
507      !
508   END SUBROUTINE ssh_nxt
509
510   !!======================================================================
511END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.