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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 13 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections
  • Property svn:keywords set to Id
File size: 17.9 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 obc_par         ! open boundary cond. parameter
28   USE obc_oce
29   USE bdy_oce
30   USE diaar5, ONLY :   lk_diaar5
31   USE iom
32   USE sbcrnf, ONLY  : h_rnf, nk_rnf  ! River runoff
33#if defined key_asminc   
34   USE asminc          ! Assimilation increment
35#endif
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   ssh_wzv    ! called by step.F90
41   PUBLIC   ssh_nxt    ! called by step.F90
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE ssh_wzv( kt ) 
55      !!----------------------------------------------------------------------
56      !!                ***  ROUTINE ssh_wzv  ***
57      !!                   
58      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
59      !!              and update the now vertical coordinate (lk_vvl=T).
60      !!
61      !! ** Method  : - Using the incompressibility hypothesis, the vertical
62      !!      velocity is computed by integrating the horizontal divergence 
63      !!      from the bottom to the surface minus the scale factor evolution.
64      !!        The boundary conditions are w=0 at the bottom (no flux) and.
65      !!
66      !! ** action  :   ssha    : after sea surface height
67      !!                wn      : now vertical velocity
68      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
69      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
70      !!
71      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
72      !!----------------------------------------------------------------------
73      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
74      !!
75      INTEGER, INTENT(in) ::   kt   ! time step
76      !!
77      INTEGER  ::   ji, jj, jk      ! dummy loop indices
78      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
79      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
80      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
81      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
82      !!----------------------------------------------------------------------
83
84      IF( kt == nit000 ) THEN
85         !
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
88         IF(lwp) WRITE(numout,*) '~~~~~~~ '
89         !
90         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
91         !
92         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
93            DO jj = 1, jpjm1
94               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
95                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
96                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
97                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
98                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
99                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
100                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
101                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
102                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
103                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
104                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
105                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
106               END DO
107            END DO
108            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
109            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
110            DO jj = 1, jpjm1
111               DO ji = 1, jpim1      ! NO Vector Opt.
112                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
113                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
114                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
115                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
116               END DO
117            END DO
118            CALL lbc_lnk( sshf_n, 'F', 1. )
119         ENDIF
120         !
121      ENDIF
122
123      !                                           !------------------------------------------!
124      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
125         !                                        !------------------------------------------!
126         DO jk = 1, jpkm1
127            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
128            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
129            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
130            !
131            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
132            fse3u (:,:,jk) = fse3u_n (:,:,jk)
133            fse3v (:,:,jk) = fse3v_n (:,:,jk)
134            fse3f (:,:,jk) = fse3f_n (:,:,jk)
135            fse3w (:,:,jk) = fse3w_n (:,:,jk)
136            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
137            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
138         END DO
139         !
140         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
141         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
142         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
143         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
144         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
145         !
146      ENDIF
147      !
148      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
149      !
150      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
151      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
152
153      !                                           !------------------------------!
154      !                                           !   After Sea Surface Height   !
155      !                                           !------------------------------!
156      zhdiv(:,:) = 0.e0
157      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
158        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
159      END DO
160      !                                                ! Sea surface elevation time stepping
161      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
162      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
163      z1_rau0 = 0.5 / rau0
164      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) &
165      &                      * tmask(:,:,1)
166
167#if defined key_obc
168      IF( Agrif_Root() ) THEN
169         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
170         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
171      ENDIF
172#endif
173#if defined key_bdy
174      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
175      CALL lbc_lnk( ssha, 'T', 1. ) 
176#endif
177
178      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
179      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
180         DO jj = 1, jpjm1
181            DO ji = 1, jpim1      ! NO Vector Opt.
182               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
183                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
184                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
185               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
186                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
187                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
188            END DO
189         END DO
190         ! Boundaries conditions
191         CALL lbc_lnk( sshu_a, 'U', 1. )
192         CALL lbc_lnk( sshv_a, 'V', 1. )
193      ENDIF
194! Include the IAU weighted SSH increment
195#if defined key_asminc
196      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
197         CALL ssh_asm_inc( kt )
198         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
199      ENDIF
200#endif
201
202      !                                           !------------------------------!
203      !                                           !     Now Vertical Velocity    !
204      !                                           !------------------------------!
205      z1_2dt = 1.e0 / z2dt
206      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
207         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
208         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
209            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
210            &                         * tmask(:,:,jk) * z1_2dt
211#if defined key_bdy
212         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
213#endif
214      END DO
215
216      !                                           !------------------------------!
217      !                                           !           outputs            !
218      !                                           !------------------------------!
219      CALL iom_put( "woce", wn                    )   ! vertical velocity
220      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
221      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
222      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
223         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
224         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
225         DO jk = 1, jpk
226            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
227         END DO
228         CALL iom_put( "w_masstr" , z3d                     ) 
229         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
230      ENDIF
231      !
232      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
233      !
234   END SUBROUTINE ssh_wzv
235
236
237   SUBROUTINE ssh_nxt( kt )
238      !!----------------------------------------------------------------------
239      !!                    ***  ROUTINE ssh_nxt  ***
240      !!
241      !! ** Purpose :   achieve the sea surface  height time stepping by
242      !!              applying Asselin time filter and swapping the arrays
243      !!              ssha  already computed in ssh_wzv 
244      !!
245      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
246      !!              from the filter, see Leclair and Madec 2010) and swap :
247      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
248      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
249      !!                sshn = ssha
250      !!
251      !! ** action  : - sshb, sshn   : before & now sea surface height
252      !!                               ready for the next time step
253      !!
254      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
255      !!----------------------------------------------------------------------
256      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
257      !!
258      INTEGER  ::   ji, jj   ! dummy loop indices
259      REAL(wp) ::   zec      ! temporary scalar
260      !!----------------------------------------------------------------------
261
262      IF( kt == nit000 ) THEN
263         IF(lwp) WRITE(numout,*)
264         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
265         IF(lwp) WRITE(numout,*) '~~~~~~~ '
266      ENDIF
267
268      !                       !--------------------------!
269      IF( lk_vvl ) THEN       !  Variable volume levels  !
270         !                    !--------------------------!
271         !
272         ! ssh at t-, u-, v, f-points
273         !===========================
274         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
275            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
276            sshu_n(:,:) = sshu_a(:,:)
277            sshv_n(:,:) = sshv_a(:,:)
278            DO jj = 1, jpjm1
279               DO ji = 1, jpim1      ! NO Vector Opt.
280                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
281                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
282                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
283                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
284               END DO
285            END DO
286            ! Boundaries conditions
287            CALL lbc_lnk( sshf_n, 'F', 1. )
288         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
289            zec = atfp * rdt / rau0
290            DO jj = 1, jpj
291               DO ji = 1, jpi                                ! before <-- now filtered
292                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
293                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
294                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
295                  sshu_n(ji,jj) = sshu_a(ji,jj)
296                  sshv_n(ji,jj) = sshv_a(ji,jj)
297               END DO
298            END DO
299            DO jj = 1, jpjm1
300               DO ji = 1, jpim1      ! NO Vector Opt.
301                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
302                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
303                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
304                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
305               END DO
306            END DO
307            ! Boundaries conditions
308            CALL lbc_lnk( sshf_n, 'F', 1. )
309            DO jj = 1, jpjm1
310               DO ji = 1, jpim1      ! NO Vector Opt.
311                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
312                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
313                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
314                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
315                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
316                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
317               END DO
318            END DO
319            ! Boundaries conditions
320            CALL lbc_lnk( sshu_b, 'U', 1. )
321            CALL lbc_lnk( sshv_b, 'V', 1. )
322         ENDIF
323         !                    !--------------------------!
324      ELSE                    !        fixed levels      !
325         !                    !--------------------------!
326         !
327         ! ssh at t-point only
328         !====================
329         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
330            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
331            !
332         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
333            DO jj = 1, jpj
334               DO ji = 1, jpi                                ! before <-- now filtered
335                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
336                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
337               END DO
338            END DO
339         ENDIF
340         !
341      ENDIF
342      !
343      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
344      !
345   END SUBROUTINE ssh_nxt
346
347   !!======================================================================
348END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.