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

Last change on this file since 4438 was 4438, checked in by trackstand2, 10 years ago

Protect debug SUM output in dynnxt and sshwzv with ARPDBGSUM
  • Property svn:keywords set to Id
File size: 25.6 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#if defined ARPDBGSUM
228      WRITE(*,*)'ARPDBG, ssh_wzv: sum WWW of hdivn=',SUM(hdivn),' at step=',kt
229      WRITE(*,*)'ARPDBG, ssh_wzv: sum WWW of fse3t=',SUM(fse3t(:,:,:)),' at step=',kt
230#endif
231      !                                           !------------------------------!
232      !                                           !   After Sea Surface Height   !
233      !                                           !------------------------------!
234      zhdiv(:,:) = 0._wp
235#if defined key_z_first
236      DO jj = 1, jpj
237         DO ji = 1, jpi
238            DO jk = 1, mbkmax(ji,jj)-1                 ! Horizontal divergence of barotropic transports
239               zhdiv(ji,jj) = zhdiv(ji,jj) + fse3t(ji,jj,jk) * hdivn(ji,jj,jk)
240            END DO
241         END DO
242      END DO
243#else
244      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
245        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
246      END DO
247#endif
248
249#if defined ARPDBGSUM
250      WRITE(*,*)'ARPDBG, ssh_wzv: sum XXX of zhdiv=',SUM(zhdiv),' at step=',kt
251      WRITE(*,*)'ARPDBG, ssh_wzv: sum XXX of ssha=',SUM(ssha),' at step=',kt
252#endif
253      !                                                ! Sea surface elevation time stepping
254      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
255      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
256      z1_rau0 = 0.5 / rau0
257#if defined key_z_first
258      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask_1(:,:)
259#else
260      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
261#endif
262
263#if defined ARPDBGSUM
264      WRITE(*,*)'ARPDBG, ssh_wzv: sum YYY of sshb=',SUM(sshb),' at step=',kt
265      WRITE(*,*)'ARPDBG, ssh_wzv: sum YYY of ssha=',SUM(ssha),' at step=',kt
266#endif
267
268#if defined key_agrif
269      CALL agrif_ssh( kt )
270#endif
271#if defined key_obc
272      IF( Agrif_Root() ) THEN
273         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
274         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
275      ENDIF
276#endif
277#if defined key_bdy
278      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
279      CALL lbc_lnk( ssha, 'T', 1. ) 
280#endif
281
282      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
283      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
284         DO jj = 1, jpjm1
285            DO ji = 1, jpim1      ! NO Vector Opt.
286#if defined key_z_first
287               sshu_a(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
288                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
289                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
290               sshv_a(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
291                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
292                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
293#else
294               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
295                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
296                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
297               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
298                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
299                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
300#endif
301            END DO
302         END DO
303         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
304      ENDIF
305     
306#if defined key_asminc
307      !                                                ! Include the IAU weighted SSH increment
308      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
309         CALL ssh_asm_inc( kt )
310         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
311      ENDIF
312#endif
313
314      !                                           !------------------------------!
315      !                                           !     Now Vertical Velocity    !
316      !                                           !------------------------------!
317      z1_2dt = 1.e0 / z2dt
318#if defined key_z_first
319      DO jj = 1, jpj
320         DO ji = 1, jpi
321            DO jk = mbkmax(ji,jj)-1, 1, -1                      ! integrate from the bottom the hor. divergence
322                wn(ji,jj,jk) = wn(ji,jj,jk+1)                               &
323                   &         -   fse3t_n(ji,jj,jk) * hdivn(ji,jj,jk)        &
324                   &         - ( fse3t_a(ji,jj,jk) - fse3t_b(ji,jj,jk) )    &
325                   &            * tmask(ji,jj,jk) * z1_2dt
326#if defined key_bdy
327                wn(ji,jj,jk) = wn(ji,jj,jk) * bdytmask(ji,jj)
328#endif
329            END DO
330         END DO
331      END DO
332#else
333      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
334         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
335         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
336            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
337            &                         * tmask(:,:,jk) * z1_2dt
338#if defined key_bdy
339         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
340#endif
341      END DO
342#endif
343
344      !                                           !------------------------------!
345      !                                           !           outputs            !
346      !                                           !------------------------------!
347      CALL iom_put( "woce", wn                    )   ! vertical velocity
348      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
349      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
350      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
351         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
352         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
353#if defined key_z_first
354         DO jj = 1, jpj
355            DO ji = 1, jpi
356               DO jk = 1, mbkmax(ji,jj)
357                  z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj)
358               END DO
359            END DO
360         END DO
361#else
362         DO jk = 1, jpk
363            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
364         END DO
365#endif
366         CALL iom_put( "w_masstr" , z3d                     ) 
367         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
368      ENDIF
369      !
370      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
371      !
372      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
373      !
374   END SUBROUTINE ssh_wzv
375
376
377   SUBROUTINE ssh_nxt( kt )
378      !!----------------------------------------------------------------------
379      !!                    ***  ROUTINE ssh_nxt  ***
380      !!
381      !! ** Purpose :   achieve the sea surface  height time stepping by
382      !!              applying Asselin time filter and swapping the arrays
383      !!              ssha  already computed in ssh_wzv 
384      !!
385      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
386      !!              from the filter, see Leclair and Madec 2010) and swap :
387      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
388      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
389      !!                sshn = ssha
390      !!
391      !! ** action  : - sshb, sshn   : before & now sea surface height
392      !!                               ready for the next time step
393      !!
394      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
395      !!----------------------------------------------------------------------
396      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
397      !!
398      INTEGER  ::   ji, jj   ! dummy loop indices
399      REAL(wp) ::   zec      ! temporary scalar
400      !!----------------------------------------------------------------------
401
402      IF( kt == nit000 ) THEN
403         IF(lwp) WRITE(numout,*)
404         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
405         IF(lwp) WRITE(numout,*) '~~~~~~~ '
406      ENDIF
407
408      !                       !--------------------------!
409      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
410         !                    !--------------------------!
411         !
412         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
413            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
414            sshu_n(:,:) = sshu_a(:,:)
415            sshv_n(:,:) = sshv_a(:,:)
416            DO jj = 1, jpjm1                                ! ssh now at f-point
417               DO ji = 1, jpim1      ! NO Vector Opt.
418#if defined key_z_first
419                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
420                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
421                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
422                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
423#else
424                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
425                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
426                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
427                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
428#endif
429               END DO
430            END DO
431            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
432            !
433         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
434
435#if defined ARPDBGSUM
436            WRITE(*,*) 'ARPDBG: ssh_nxt: SUM of sshb = ',SUM(sshb),' at step=',kt
437            WRITE(*,*) 'ARPDBG: ssh_nxt: SUM of sshn = ',SUM(sshn),' at step=',kt
438#endif
439
440            zec = atfp * rdt / rau0
441            DO jj = 1, jpj
442               DO ji = 1, jpi                               ! before <-- now filtered
443#if defined key_z_first
444                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
445                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask_1(ji,jj)
446#else
447                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
448                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
449#endif
450                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
451                  sshu_n(ji,jj) = sshu_a(ji,jj)
452                  sshv_n(ji,jj) = sshv_a(ji,jj)
453               END DO
454            END DO
455            DO jj = 1, jpjm1                                ! ssh now at f-point
456               DO ji = 1, jpim1      ! NO Vector Opt.
457#if defined key_z_first
458                  sshf_n(ji,jj) = 0.5  * umask_1(ji,jj) * umask_1(ji,jj+1)                 &
459                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
460                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
461                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
462#else
463                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
464                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
465                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
466                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
467#endif
468               END DO
469            END DO
470            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
471            !
472            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
473               DO ji = 1, jpim1      ! NO Vector Opt.
474#if defined key_z_first
475                  sshu_b(ji,jj) = 0.5  * umask_1(ji,jj) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
476                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
477                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
478                  sshv_b(ji,jj) = 0.5  * vmask_1(ji,jj) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
479                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
480                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
481#else
482                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
483                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
484                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
485                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
486                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
487                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
488#endif
489               END DO
490            END DO
491            CALL lbc_lnk( sshu_b, 'U', 1. )
492            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
493            !
494         ENDIF
495         !                    !--------------------------!
496      ELSE                    !        fixed levels      !     (ssh at t-point only)
497         !                    !--------------------------!
498         !
499         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
500            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
501            !
502         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
503            DO jj = 1, jpj
504               DO ji = 1, jpi                               ! before <-- now filtered
505                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
506                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
507               END DO
508            END DO
509         ENDIF
510         !
511      ENDIF
512      !
513      ! Update velocity at AGRIF zoom boundaries
514#if defined key_agrif
515      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
516#endif
517      !
518      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
519      !
520   END SUBROUTINE ssh_nxt
521
522   !!======================================================================
523END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.