New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sshwzv.F90 in branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN – NEMO

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2148

Last change on this file since 2148 was 2148, checked in by cetlod, 13 years ago

merge LOCEAN 2010 developments branches

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