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/2010_and_older/DEV_r1837_MLF/NEMO/OPA_SRC/DYN – NEMO

source: branches/2010_and_older/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 5776

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

ticket: #663 allocatable variables in diahsb and remove unnecessary fse3t_d field

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.3 KB
Line 
1
2MODULE sshwzv   
3   !!==============================================================================
4   !!                       ***  MODULE  sshwzv  ***
5   !! Ocean dynamics : sea surface height and vertical velocity
6   !!==============================================================================
7   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
8   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   ssh_wzv        : after ssh & now vertical velocity
13   !!   ssh_nxt        : filter ans swap the ssh arrays
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers variables
16   USE dom_oce         ! ocean space and time domain variables
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE domvvl          ! Variable volume
19   USE divcur          ! hor. divergence and curl      (div & cur routines)
20   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
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 diaar5, ONLY :   lk_diaar5
30   USE iom
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                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
144      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
145
146      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
147      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
148
149      !                                           !------------------------------!
150      !                                           !   After Sea Surface Height   !
151      !                                           !------------------------------!
152      zhdiv(:,:) = 0.e0
153      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
154        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
155      END DO
156      !                                                ! Sea surface elevation time stepping
157      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
158      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
159      z1_rau0 = 0.5 / rau0
160      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
161
162#if defined key_obc
163      IF( Agrif_Root() ) THEN
164         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
165         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
166      ENDIF
167#endif
168      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
169      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
170         DO jj = 1, jpjm1
171            DO ji = 1, jpim1      ! NO Vector Opt.
172               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
173                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
174                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
175               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
176                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
177                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
178            END DO
179         END DO
180         ! Boundaries conditions
181         CALL lbc_lnk( sshu_a, 'U', 1. )
182         CALL lbc_lnk( sshv_a, 'V', 1. )
183      ENDIF
184      !                                           !------------------------------!
185      !                                           !     Now Vertical Velocity    !
186      !                                           !------------------------------!
187      z1_2dt = 1.e0 / z2dt
188      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
189         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
190         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
191            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
192            &                         * tmask(:,:,jk) * z1_2dt
193      END DO
194
195      !                                           !------------------------------!
196      !                                           !           outputs            !
197      !                                           !------------------------------!
198      CALL iom_put( "woce", wn                    )   ! vertical velocity
199      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
200      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
201      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
202         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
203         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
204         DO jk = 1, jpk
205            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
206         END DO
207         CALL iom_put( "w_masstr" , z3d                     ) 
208         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
209      ENDIF
210      !
211      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
212      !
213   END SUBROUTINE ssh_wzv
214
215
216   SUBROUTINE ssh_nxt( kt )
217      !!----------------------------------------------------------------------
218      !!                    ***  ROUTINE ssh_nxt  ***
219      !!
220      !! ** Purpose :   achieve the sea surface  height time stepping by
221      !!              applying Asselin time filter and swapping the arrays
222      !!              ssha  already computed in ssh_wzv 
223      !!
224      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
225      !!              from the filter, see Leclair and Madec 2010) and swap :
226      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
227      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
228      !!                sshn = ssha
229      !!
230      !! ** action  : - sshb, sshn   : before & now sea surface height
231      !!                               ready for the next time step
232      !!
233      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
234      !!----------------------------------------------------------------------
235      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
236      !!
237      INTEGER  ::   ji, jj   ! dummy loop indices
238      REAL(wp) ::   zec      ! temporary scalar
239      !!----------------------------------------------------------------------
240
241      IF( kt == nit000 ) THEN
242         IF(lwp) WRITE(numout,*)
243         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
244         IF(lwp) WRITE(numout,*) '~~~~~~~ '
245      ENDIF
246
247      !                       !--------------------------!
248      IF( lk_vvl ) THEN       !  Variable volume levels  !
249         !                    !--------------------------!
250         !
251         ! ssh at t-, u-, v, f-points
252         !===========================
253         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
254            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
255            sshu_n(:,:) = sshu_a(:,:)
256            sshv_n(:,:) = sshv_a(:,:)
257            DO jj = 1, jpjm1
258               DO ji = 1, jpim1      ! NO Vector Opt.
259                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
260                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
261                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
262                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
263               END DO
264            END DO
265            ! Boundaries conditions
266            CALL lbc_lnk( sshf_n, 'F', 1. )
267         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
268            zec = atfp * rdt / rau0
269            DO jj = 1, jpj
270               DO ji = 1, jpi                                ! before <-- now filtered
271                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )   &
272                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
273                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
274                  sshu_n(ji,jj) = sshu_a(ji,jj)
275                  sshv_n(ji,jj) = sshv_a(ji,jj)
276               END DO
277            END DO
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            DO jj = 1, jpjm1
289               DO ji = 1, jpim1      ! NO Vector Opt.
290                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
291                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
292                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
293                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
294                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
295                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
296               END DO
297            END DO
298            ! Boundaries conditions
299            CALL lbc_lnk( sshu_b, 'U', 1. )
300            CALL lbc_lnk( sshv_b, 'V', 1. )
301         ENDIF
302         !                    !--------------------------!
303      ELSE                    !        fixed levels      !
304         !                    !--------------------------!
305         !
306         ! ssh at t-point only
307         !====================
308         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
309            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
310            !
311         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
312            DO jj = 1, jpj
313               DO ji = 1, jpi                                ! before <-- now filtered
314                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
315                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
316               END DO
317            END DO
318         ENDIF
319         !
320      ENDIF
321      !
322      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
323      !
324   END SUBROUTINE ssh_nxt
325
326   !!======================================================================
327END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.