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

Last change on this file since 1870 was 1870, checked in by gm, 14 years ago

ticket: #663 step-1 : introduce the modified forcing term

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 15.6 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, zraur     ! 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                  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                           !  Regridding: 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         !
130         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
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      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
142      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
143
144      !                                           !------------------------------!
145      !                                           !   After Sea Surface Height   !
146      !                                           !------------------------------!
147      zraur = 0.5 / rau0
148      !
149      zhdiv(:,:) = 0.e0
150      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
151        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
152      END DO
153      !                                                ! Sea surface elevation time stepping
154      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
155
156#if defined key_obc
157      IF( Agrif_Root() ) THEN
158         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
159         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
160      ENDIF
161#endif
162      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
163      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
164         DO jj = 1, jpjm1
165            DO ji = 1, jpim1      ! NO Vector Opt.
166               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
167                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
168                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
169               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
170                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
171                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
172               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
173                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
174                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
175            END DO
176         END DO
177         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
178         CALL lbc_lnk( sshv_a, 'V', 1. )
179         CALL lbc_lnk( sshf_a, 'F', 1. )
180      ENDIF
181
182      !                                           !------------------------------!
183      !                                           !     Now Vertical Velocity    !
184      !                                           !------------------------------!
185      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
186         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)     &
187            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
188      END DO
189
190      !                                           !------------------------------!
191      !                                           !           outputs            !
192      !                                           !------------------------------!
193      CALL iom_put( "woce", wn                    )   ! vertical velocity
194      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
195      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
196      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
197         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
198         DO jk = 1, jpk
199            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
200         END DO
201         CALL iom_put( "w_masstr" , z3d                     ) 
202         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
203      ENDIF
204      !
205      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
206      !
207   END SUBROUTINE ssh_wzv
208
209
210   SUBROUTINE ssh_nxt( kt )
211      !!----------------------------------------------------------------------
212      !!                    ***  ROUTINE ssh_nxt  ***
213      !!
214      !! ** Purpose :   achieve the sea surface  height time stepping by
215      !!              applying Asselin time filter and swapping the arrays
216      !!              ssha  already computed in ssh_wzv 
217      !!
218      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
219      !!              from the filter, see Leclair and Madec 2010) and swap :
220      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
221      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
222      !!                sshn = ssha
223      !!
224      !! ** action  : - sshb, sshn   : before & now sea surface height
225      !!                               ready for the next time step
226      !!
227      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
228      !!----------------------------------------------------------------------
229      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
230      !!
231      INTEGER  ::   ji, jj   ! dummy loop indices
232      REAL(wp) ::   zec      ! temporary scalare
233      !!----------------------------------------------------------------------
234
235      IF( kt == nit000 ) THEN
236         IF(lwp) WRITE(numout,*)
237         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
238         IF(lwp) WRITE(numout,*) '~~~~~~~ '
239      ENDIF
240
241      !                       !--------------------------!
242      IF( lk_vvl ) THEN       !  Variable volume levels  !   ssh at t-, u-, v, f-points
243         !                    !--------------------------!
244         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
245            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
246            sshu_n(:,:) = sshu_a(:,:)
247            sshv_n(:,:) = sshv_a(:,:)
248            sshf_n(:,:) = sshf_a(:,:)
249         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
250            DO jj = 1, jpj
251               DO ji = 1, jpi                                ! before <-- now filtered
252                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
253                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
254                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
255                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
256                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
257                  sshu_n(ji,jj) = sshu_a(ji,jj)
258                  sshv_n(ji,jj) = sshv_a(ji,jj)
259                  sshf_n(ji,jj) = sshf_a(ji,jj)
260               END DO
261            END DO
262         ENDIF
263         !                    !--------------------------!
264      ELSE                    !        fixed levels      !   ssh at t-point only
265         !                    !--------------------------!
266         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
267            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
268            !
269         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
270            zec = atfp * rdt / rau0
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                     &                      - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
275                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
276               END DO
277            END DO
278         ENDIF
279         !
280      ENDIF
281
282      ! time filter with global conservation correction and array swap
283      sshbb(:,:) = sshb(:,:)
284      sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) )    &
285           &       - zfact  * 
286      sshn (:,:) = zssha(:,:)
287      empb (:,:) = emp(:,:)
288
289
290      !
291      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
292      !
293   END SUBROUTINE ssh_nxt
294
295   !!======================================================================
296END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.