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

source: trunk/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 1756

Last change on this file since 1756 was 1756, checked in by smasson, 14 years ago

implement AR5 diagnostics, see ticket:610
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.5 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   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   ssh_wzv        : after ssh & now vertical velocity
11   !!   ssh_nxt        : filter ans swap the ssh arrays
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and tracers variables
14   USE dom_oce         ! ocean space and time domain variables
15   USE sbc_oce         ! surface boundary condition: ocean
16   USE domvvl          ! Variable volume
17   USE divcur          ! hor. divergence and curl      (div & cur routines)
18   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
19   USE iom             ! I/O library
20   USE restart         ! only for lrst_oce
21   USE in_out_manager  ! I/O manager
22   USE prtctl          ! Print control
23   USE phycst
24   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
25   USE obc_par         ! open boundary cond. parameter
26   USE obc_oce
27   USE diaar5, ONLY :   lk_diaar5
28   USE iom
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ssh_wzv    ! called by step.F90
34   PUBLIC   ssh_nxt    ! called by step.F90
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
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  : -
56      !!
57      !!              - Using the incompressibility hypothesis, the vertical
58      !!      velocity is computed by integrating the horizontal divergence 
59      !!      from the bottom to the surface minus the scale factor evolution.
60      !!        The boundary conditions are w=0 at the bottom (no flux) and.
61      !!
62      !! ** action  :   ssha    : after sea surface height
63      !!                wn      : now vertical velocity
64      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
65      !!                          at u-, v-, f-point s
66      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
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, zraur     ! 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         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                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
97                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
98                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
99                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
100                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
101                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
102                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
103                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
104               END DO
105            END DO
106            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
107            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
108            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
109         ENDIF
110         !
111      ENDIF
112
113      !                                           !------------------------------!
114      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
115      !                                           !------------------------------!
116         DO jk = 1, jpkm1
117            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
118            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
119            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
120            !
121            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
122            fse3u (:,:,jk) = fse3u_n (:,:,jk)
123            fse3v (:,:,jk) = fse3v_n (:,:,jk)
124            fse3f (:,:,jk) = fse3f_n (:,:,jk)
125            fse3w (:,:,jk) = fse3w_n (:,:,jk)
126            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
127            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
128         END DO
129         !                                             ! now ocean depth (at u- and v-points)
130         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
131         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
132         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
133         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
134         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
135         !
136      ENDIF
137
138                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
139      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
140
141      ! set time step size (Euler/Leapfrog)
142      z2dt = 2. * rdt 
143      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
144
145      zraur = 1. / rau0
146
147      !                                           !------------------------------!
148      !                                           !   After Sea Surface Height   !
149      !                                           !------------------------------!
150      zhdiv(:,:) = 0.e0
151      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
152        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
153      END DO
154
155      !                                                ! Sea surface elevation time stepping
156      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
157
158#if defined key_obc
159# if defined key_agrif
160      IF ( Agrif_Root() ) THEN 
161# endif
162         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
163         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
164# if defined key_agrif
165      ENDIF
166# endif
167#endif
168
169      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
170      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
171         DO jj = 1, jpjm1
172            DO ji = 1, jpim1      ! NO Vector Opt.
173               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
174                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
175                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
176               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
177                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
178                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
179               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
180                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
181                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
182            END DO
183         END DO
184         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
185         CALL lbc_lnk( sshv_a, 'V', 1. )
186         CALL lbc_lnk( sshf_a, 'F', 1. )
187      ENDIF
188
189      !                                           !------------------------------!
190      !                                           !     Now Vertical Velocity    !
191      !                                           !------------------------------!
192      !                                                ! integrate from the bottom the hor. divergence
193      DO jk = jpkm1, 1, -1
194         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
195              &                    - (  fse3t_a(:,:,jk)                   &
196              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
197      END DO
198      !
199      CALL iom_put( "woce", wn                    )   ! vertical velocity
200      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
201      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
202      IF( lk_diaar5 ) THEN
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                     )   !           vertical mass transport
208         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
209      ENDIF
210      !
211   END SUBROUTINE ssh_wzv
212
213
214   SUBROUTINE ssh_nxt( kt )
215      !!----------------------------------------------------------------------
216      !!                    ***  ROUTINE ssh_nxt  ***
217      !!
218      !! ** Purpose :   achieve the sea surface  height time stepping by
219      !!              applying Asselin time filter and swapping the arrays
220      !!              ssha  already computed in ssh_wzv 
221      !!
222      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
223      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
224      !!             sshn = ssha
225      !!
226      !! ** action  : - sshb, sshn   : before & now sea surface height
227      !!                               ready for the next time step
228      !!----------------------------------------------------------------------
229      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
230      !!
231      INTEGER  ::   ji, jj               ! dummy loop indices
232      !!----------------------------------------------------------------------
233
234      IF( kt == nit000 ) THEN
235         IF(lwp) WRITE(numout,*)
236         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
237         IF(lwp) WRITE(numout,*) '~~~~~~~ '
238      ENDIF
239
240      ! Time filter and swap of the ssh
241      ! -------------------------------
242
243      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
244         !                   ! ---------------------- !
245         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
246            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
247            sshu_n(:,:) = sshu_a(:,:)
248            sshv_n(:,:) = sshv_a(:,:)
249            sshf_n(:,:) = sshf_a(:,:)
250         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
251            DO jj = 1, jpj
252               DO ji = 1, jpi                                ! before <-- now filtered
253                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
254                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
255                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
256                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
257                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
258                  sshu_n(ji,jj) = sshu_a(ji,jj)
259                  sshv_n(ji,jj) = sshv_a(ji,jj)
260                  sshf_n(ji,jj) = sshf_a(ji,jj)
261               END DO
262            END DO
263         ENDIF
264         !
265      ELSE                   ! fixed levels :   ssh at t-point only
266         !                   ! ------------ !
267         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
268            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
269            !
270         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
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                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
275               END DO
276            END DO
277         ENDIF
278         !
279      ENDIF
280      !
281      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
282      !
283   END SUBROUTINE ssh_nxt
284
285   !!======================================================================
286END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.