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 on Ticket #1438 – Attachment – NEMO

Ticket #1438: sshwzv.F90

File sshwzv.F90, 14.9 KB (added by Robin_Waldman, 9 years ago)
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.2  !  2015-01 (R. Waldman) initialize sshb, sshn if key_istate_uvg
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   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
43   !! $Id: sshwzv.F90 1756 2009-11-25 14:15:20Z smasson $
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  : -
57      !!
58      !!              - Using the incompressibility hypothesis, the vertical
59      !!      velocity is computed by integrating the horizontal divergence 
60      !!      from the bottom to the surface minus the scale factor evolution.
61      !!        The boundary conditions are w=0 at the bottom (no flux) and.
62      !!
63      !! ** action  :   ssha    : after sea surface height
64      !!                wn      : now vertical velocity
65      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
66      !!                          at u-, v-, f-point s
67      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
68      !!----------------------------------------------------------------------
69      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
70      !!
71      INTEGER, INTENT(in) ::   kt   ! time step
72      !!
73      INTEGER  ::   ji, jj, jk      ! dummy loop indices
74      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
75      REAL(wp) ::   z2dt, zraur     ! temporary scalars
76      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
77      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
78      !!----------------------------------------------------------------------
79
80      IF( kt == nit000 ) THEN
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                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
98                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
99                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
100                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
101                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
102                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
103                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
104                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
105               END DO
106            END DO
107            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
108            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
109            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
110         ENDIF
111         !
112      ENDIF
113
114      !                                           !------------------------------!
115      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
116      !                                           !------------------------------!
117         DO jk = 1, jpkm1
118            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
119            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
120            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
121            !
122            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
123            fse3u (:,:,jk) = fse3u_n (:,:,jk)
124            fse3v (:,:,jk) = fse3v_n (:,:,jk)
125            fse3f (:,:,jk) = fse3f_n (:,:,jk)
126            fse3w (:,:,jk) = fse3w_n (:,:,jk)
127            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
128            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
129         END DO
130         !                                             ! now ocean depth (at u- and v-points)
131         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
132         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
133         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
134         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
135         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
136         !
137      ENDIF
138
139                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
140      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
141
142      ! set time step size (Euler/Leapfrog)
143      z2dt = 2. * rdt 
144      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
145
146      zraur = 1. / rau0
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
156      !                                                ! Sea surface elevation time stepping
157
158      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
159
160#if defined key_obc
161# if defined key_agrif
162      IF ( Agrif_Root() ) THEN 
163# endif
164         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
165         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
166# if defined key_agrif
167      ENDIF
168# endif
169#endif
170
171      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
172      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
173         DO jj = 1, jpjm1
174            DO ji = 1, jpim1      ! NO Vector Opt.
175               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
176                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
177                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
178               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
179                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
180                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
181               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
182                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
183                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
184            END DO
185         END DO
186         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
187         CALL lbc_lnk( sshv_a, 'V', 1. )
188         CALL lbc_lnk( sshf_a, 'F', 1. )
189      ENDIF
190
191      !                                           !------------------------------!
192      !                                           !     Now Vertical Velocity    !
193      !                                           !------------------------------!
194      !                                                ! integrate from the bottom the hor. divergence
195
196      DO jk = jpkm1, 1, -1
197         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
198              &                    - (  fse3t_a(:,:,jk)                   &
199              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
200      END DO
201
202      !
203         CALL iom_put( "woce", wn                    )   ! vertical velocity
204         CALL iom_put( "ssh" , sshn                  )   ! sea surface height
205         CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
206
207      IF( lk_diaar5 ) THEN
208         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
209         DO jk = 1, jpk
210            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
211         END DO
212         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
213         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
214      ENDIF
215      !
216
217#if defined key_istate_uvg
218      ! At 1st timestep put all ssh to ssha value in order to initialize ssh:
219      IF( kt == nit000 ) THEN
220       sshb(:,:)=ssha(:,:)
221       sshn(:,:)=ssha(:,:)
222      ENDIF
223#endif
224
225   END SUBROUTINE ssh_wzv
226
227
228   SUBROUTINE ssh_nxt( kt )
229      !!----------------------------------------------------------------------
230      !!                    ***  ROUTINE ssh_nxt  ***
231      !!
232      !! ** Purpose :   achieve the sea surface  height time stepping by
233      !!              applying Asselin time filter and swapping the arrays
234      !!              ssha  already computed in ssh_wzv 
235      !!
236      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
237      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
238      !!             sshn = ssha
239      !!
240      !! ** action  : - sshb, sshn   : before & now sea surface height
241      !!                               ready for the next time step
242      !!----------------------------------------------------------------------
243      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
244      !!
245      INTEGER  ::   ji, jj               ! dummy loop indices
246      !!----------------------------------------------------------------------
247
248      IF( kt == nit000 ) THEN
249         IF(lwp) WRITE(numout,*)
250         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
251         IF(lwp) WRITE(numout,*) '~~~~~~~ '
252      ENDIF
253
254      ! Time filter and swap of the ssh
255      ! -------------------------------
256
257      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
258         !                   ! ---------------------- !
259         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
260            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
261            sshu_n(:,:) = sshu_a(:,:)
262            sshv_n(:,:) = sshv_a(:,:)
263            sshf_n(:,:) = sshf_a(:,:)
264         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
265            DO jj = 1, jpj
266               DO ji = 1, jpi                                ! before <-- now filtered
267                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
268                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
269                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
270                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
271                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
272                  sshu_n(ji,jj) = sshu_a(ji,jj)
273                  sshv_n(ji,jj) = sshv_a(ji,jj)
274                  sshf_n(ji,jj) = sshf_a(ji,jj)
275               END DO
276            END DO
277         ENDIF
278         !
279      ELSE                   ! fixed levels :   ssh at t-point only
280         !                   ! ------------ !
281         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
282            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
283            !
284         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
285            DO jj = 1, jpj
286               DO ji = 1, jpi                                ! before <-- now filtered
287                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   
288                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
289               END DO
290            END DO
291         ENDIF
292         !
293      ENDIF
294      !
295      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
296      !
297   END SUBROUTINE ssh_nxt
298
299   !!======================================================================
300END MODULE sshwzv