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

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2113

Last change on this file since 2113 was 2113, checked in by cetlod, 13 years ago

phasing DEV_r2106_LOCEAN2010 with DEV_r1821_Rivers branch

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.9 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   USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep  ! River runoff
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$
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         !
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               rnf_dep(ji,jj) = 0.
141               DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth
142                  rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box
143               ENDDO
144            ENDDO
145         ENDDO
146         !
147      ENDIF
148
149                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
150      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
151
152      ! set time step size (Euler/Leapfrog)
153      z2dt = 2. * rdt 
154      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
155
156      zraur = 1. / rau0
157
158      !                                           !------------------------------!
159      !                                           !   After Sea Surface Height   !
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
166      !                                                ! Sea surface elevation time stepping
167      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
168
169#if defined key_obc
170      IF ( Agrif_Root() ) THEN
171         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
172         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
173      ENDIF
174#endif
175
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               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
187                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
188                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
189            END DO
190         END DO
191         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
192         CALL lbc_lnk( sshv_a, 'V', 1. )
193         CALL lbc_lnk( sshf_a, 'F', 1. )
194      ENDIF
195
196      !                                           !------------------------------!
197      !                                           !     Now Vertical Velocity    !
198      !                                           !------------------------------!
199      !                                                ! integrate from the bottom the hor. divergence
200      DO jk = jpkm1, 1, -1
201         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
202              &                    - (  fse3t_a(:,:,jk)                   &
203              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
204      END DO
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
210         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
211         DO jk = 1, jpk
212            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
213         END DO
214         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
215         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
216      ENDIF
217      !
218   END SUBROUTINE ssh_wzv
219
220
221   SUBROUTINE ssh_nxt( kt )
222      !!----------------------------------------------------------------------
223      !!                    ***  ROUTINE ssh_nxt  ***
224      !!
225      !! ** Purpose :   achieve the sea surface  height time stepping by
226      !!              applying Asselin time filter and swapping the arrays
227      !!              ssha  already computed in ssh_wzv 
228      !!
229      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
230      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
231      !!             sshn = ssha
232      !!
233      !! ** action  : - sshb, sshn   : before & now sea surface height
234      !!                               ready for the next time step
235      !!----------------------------------------------------------------------
236      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
237      !!
238      INTEGER  ::   ji, jj               ! dummy loop indices
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      ! Time filter and swap of the ssh
248      ! -------------------------------
249
250      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
251         !                   ! ---------------------- !
252         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
253            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
254            sshu_n(:,:) = sshu_a(:,:)
255            sshv_n(:,:) = sshv_a(:,:)
256            sshf_n(:,:) = sshf_a(:,:)
257         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
258            DO jj = 1, jpj
259               DO ji = 1, jpi                                ! before <-- now filtered
260                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
261                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
262                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
263                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
264                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
265                  sshu_n(ji,jj) = sshu_a(ji,jj)
266                  sshv_n(ji,jj) = sshv_a(ji,jj)
267                  sshf_n(ji,jj) = sshf_a(ji,jj)
268               END DO
269            END DO
270         ENDIF
271         !
272      ELSE                   ! fixed levels :   ssh at t-point only
273         !                   ! ------------ !
274         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
275            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
276            !
277         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
278            DO jj = 1, jpj
279               DO ji = 1, jpi                                ! before <-- now filtered
280                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   
281                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
282               END DO
283            END DO
284         ENDIF
285         !
286      ENDIF
287      !
288      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
289      !
290   END SUBROUTINE ssh_nxt
291
292   !!======================================================================
293END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.