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

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2005

Last change on this file since 2005 was 2005, checked in by mlelod, 14 years ago

ticket: #663 MLF: second part (local compatibility essentially)

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