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

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

cosmetic changes again

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.8 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_n(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_n(ji,jj  )     &
109                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
110               END DO
111            END DO
112            CALL lbc_lnk( sshf_n, 'F', 1. )
113         ENDIF
114         !
115      ENDIF
116
117      !                                           !------------------------------------------!
118      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
119         !                                        !------------------------------------------!
120         DO jk = 1, jpkm1
121            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
122            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
123            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
124            !
125            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
126            fse3u (:,:,jk) = fse3u_n (:,:,jk)
127            fse3v (:,:,jk) = fse3v_n (:,:,jk)
128            fse3f (:,:,jk) = fse3f_n (:,:,jk)
129            fse3w (:,:,jk) = fse3w_n (:,:,jk)
130            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
131            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
132         END DO
133         !
134         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
135         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
136         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
137         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
138         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
139         !
140      ENDIF
141
142                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
143      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
144
145      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
146      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
147
148      !                                           !------------------------------!
149      !                                           !   After Sea Surface Height   !
150      !                                           !------------------------------!
151      zhdiv(:,:) = 0.e0
152      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
153        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
154      END DO
155      !                                                ! Sea surface elevation time stepping
156      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
157      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
158      z1_rau0 = 0.5 / rau0
159      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
160
161#if defined key_obc
162      IF( Agrif_Root() ) THEN
163         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
164         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
165      ENDIF
166#endif
167      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
168      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
169         DO jj = 1, jpjm1
170            DO ji = 1, jpim1      ! NO Vector Opt.
171               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
172                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
173                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
174               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
175                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
176                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
177            END DO
178         END DO
179         ! Boundaries conditions
180         CALL lbc_lnk( sshu_a, 'U', 1. )
181         CALL lbc_lnk( sshv_a, 'V', 1. )
182      ENDIF
183      !                                           !----------------------------------------!
184      !                                           !     vertical scale factor laplacian    !
185      !                                           !----------------------------------------!
186      ! Needed for Robert-Asselin time filter and for Brown & Campana semi implicit hydrostatic presure gradient
187      DO jk = 1, jpk
188         fse3t_d(:,:,jk) =          fse3t_b(:,:,jk)   &
189            &              - 2.e0 * fse3t_n(:,:,jk)   &
190            &              +        fse3t_a(:,:,jk)
191      ENDDO
192      !                                           !------------------------------!
193      !                                           !     Now Vertical Velocity    !
194      !                                           !------------------------------!
195      z1_2dt = 1.e0 / z2dt
196      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
197         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
198         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
199            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
200            &                         * tmask(:,:,jk) * z1_2dt
201      END DO
202
203      !                                           !------------------------------!
204      !                                           !           outputs            !
205      !                                           !------------------------------!
206      CALL iom_put( "woce", wn                    )   ! vertical velocity
207      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
208      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
209      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
210         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
211         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
212         DO jk = 1, jpk
213            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
214         END DO
215         CALL iom_put( "w_masstr" , z3d                     ) 
216         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
217      ENDIF
218      !
219      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
220      !
221   END SUBROUTINE ssh_wzv
222
223
224   SUBROUTINE ssh_nxt( kt )
225      !!----------------------------------------------------------------------
226      !!                    ***  ROUTINE ssh_nxt  ***
227      !!
228      !! ** Purpose :   achieve the sea surface  height time stepping by
229      !!              applying Asselin time filter and swapping the arrays
230      !!              ssha  already computed in ssh_wzv 
231      !!
232      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
233      !!              from the filter, see Leclair and Madec 2010) and swap :
234      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
235      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
236      !!                sshn = ssha
237      !!
238      !! ** action  : - sshb, sshn   : before & now sea surface height
239      !!                               ready for the next time step
240      !!
241      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
242      !!----------------------------------------------------------------------
243      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
244      !!
245      INTEGER  ::   ji, jj   ! dummy loop indices
246      REAL(wp) ::   zec      ! temporary scalar
247      !!----------------------------------------------------------------------
248
249      IF( kt == nit000 ) THEN
250         IF(lwp) WRITE(numout,*)
251         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
252         IF(lwp) WRITE(numout,*) '~~~~~~~ '
253      ENDIF
254
255      !                       !--------------------------!
256      IF( lk_vvl ) THEN       !  Variable volume levels  !
257         !                    !--------------------------!
258         !
259         ! ssh at t-, u-, v, f-points
260         !===========================
261         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
262            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
263            sshu_n(:,:) = sshu_a(:,:)
264            sshv_n(:,:) = sshv_a(:,:)
265            DO jj = 1, jpjm1
266               DO ji = 1, jpim1      ! NO Vector Opt.
267                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
268                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
269                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
270                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
271               END DO
272            END DO
273            ! Boundaries conditions
274            CALL lbc_lnk( sshf_n, 'F', 1. )
275         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
276            zec = atfp * rdt / rau0
277            DO jj = 1, jpj
278               DO ji = 1, jpi                                ! before <-- now filtered
279                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )   &
280                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
281                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
282                  sshu_n(ji,jj) = sshu_a(ji,jj)
283                  sshv_n(ji,jj) = sshv_a(ji,jj)
284               END DO
285            END DO
286            DO jj = 1, jpjm1
287               DO ji = 1, jpim1      ! NO Vector Opt.
288                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
289                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
290                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
291                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
292               END DO
293            END DO
294            ! Boundaries conditions
295            CALL lbc_lnk( sshf_n, 'F', 1. )
296            DO jj = 1, jpjm1
297               DO ji = 1, jpim1      ! NO Vector Opt.
298                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
299                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
300                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
301                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
302                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
303                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
304               END DO
305            END DO
306            ! Boundaries conditions
307            CALL lbc_lnk( sshu_b, 'U', 1. )
308            CALL lbc_lnk( sshv_b, 'V', 1. )
309         ENDIF
310         !                    !--------------------------!
311      ELSE                    !        fixed levels      !
312         !                    !--------------------------!
313         !
314         ! ssh at t-point only
315         !====================
316         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
317            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
318            !
319         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
320            DO jj = 1, jpj
321               DO ji = 1, jpi                                ! before <-- now filtered
322                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
323                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
324               END DO
325            END DO
326         ENDIF
327         !
328      ENDIF
329      !
330      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
331      !
332   END SUBROUTINE ssh_nxt
333
334   !!======================================================================
335END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.