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

source: branches/devukmo2010/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2128

Last change on this file since 2128 was 2128, checked in by rfurner, 14 years ago

merged branches OBS, ASM, Rivers, BDY & mixed_dynldf ready for vn3.3 merge
  • 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-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
8   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
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 bdy_oce
30   USE diaar5, ONLY :   lk_diaar5
31   USE iom
32   USE sbcrnf, ONLY  : rnf_dep, rnf_mod_dep  ! River runoff
33#if defined key_asminc   
34   USE asminc          ! Assimilation increment
35#endif
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   ssh_wzv    ! called by step.F90
41   PUBLIC   ssh_nxt    ! called by step.F90
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
49   !! $Id$
50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE ssh_wzv( kt ) 
56      !!----------------------------------------------------------------------
57      !!                ***  ROUTINE ssh_wzv  ***
58      !!                   
59      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
60      !!              and update the now vertical coordinate (lk_vvl=T).
61      !!
62      !! ** Method  : -
63      !!
64      !!              - Using the incompressibility hypothesis, the vertical
65      !!      velocity is computed by integrating the horizontal divergence 
66      !!      from the bottom to the surface minus the scale factor evolution.
67      !!        The boundary conditions are w=0 at the bottom (no flux) and.
68      !!
69      !! ** action  :   ssha    : after sea surface height
70      !!                wn      : now vertical velocity
71      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
72      !!                          at u-, v-, f-point s
73      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
74      !!----------------------------------------------------------------------
75      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
76      !!
77      INTEGER, INTENT(in) ::   kt   ! time step
78      !!
79      INTEGER  ::   ji, jj, jk      ! dummy loop indices
80      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
81      REAL(wp) ::   z2dt, zraur     ! temporary scalars
82      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
83      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
84      !!----------------------------------------------------------------------
85
86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
89         IF(lwp) WRITE(numout,*) '~~~~~~~ '
90         !
91         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
92         !
93         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
94            DO jj = 1, jpjm1
95               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
96                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
97                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
98                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
99                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
100                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
101                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
102                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
103                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
104                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
105                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
106                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
107                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
108                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
109                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
110                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
111               END DO
112            END DO
113            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
114            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
115            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
116         ENDIF
117         !
118      ENDIF
119
120      !                                           !------------------------------!
121      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
122      !                                           !------------------------------!
123         DO jk = 1, jpkm1
124            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
125            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
126            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
127            !
128            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
129            fse3u (:,:,jk) = fse3u_n (:,:,jk)
130            fse3v (:,:,jk) = fse3v_n (:,:,jk)
131            fse3f (:,:,jk) = fse3f_n (:,:,jk)
132            fse3w (:,:,jk) = fse3w_n (:,:,jk)
133            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
134            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
135         END DO
136         !                                             ! now ocean depth (at u- and v-points)
137         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
138         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
139         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
140         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
141         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
142         !
143         DO jj=1,jpj 
144           DO ji=1,jpi 
145             rnf_dep(ji,jj)=0 
146             DO jk=1,rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth
147               rnf_dep(ji,jj)=rnf_dep(ji,jj)+fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box
148             ENDDO 
149           ENDDO 
150         ENDDO 
151         !
152      ENDIF
153
154                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
155      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
156
157      ! set time step size (Euler/Leapfrog)
158      z2dt = 2. * rdt 
159      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
160
161      zraur = 1. / rau0
162
163      !                                           !------------------------------!
164      !                                           !   After Sea Surface Height   !
165      !                                           !------------------------------!
166      zhdiv(:,:) = 0.e0
167      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
168        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
169      END DO
170
171      !                                                ! Sea surface elevation time stepping
172      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
173
174#if defined key_obc
175      IF ( Agrif_Root() ) THEN
176         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
177         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
178      ENDIF
179#endif
180
181#if defined key_bdy
182      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
183      CALL lbc_lnk( ssha, 'T', 1. ) 
184#endif
185
186      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
187      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
188         DO jj = 1, jpjm1
189            DO ji = 1, jpim1      ! NO Vector Opt.
190               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
191                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
192                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
193               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
194                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
195                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
196               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
197                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
198                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
199            END DO
200         END DO
201         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
202         CALL lbc_lnk( sshv_a, 'V', 1. )
203         CALL lbc_lnk( sshf_a, 'F', 1. )
204      ENDIF
205
206! Include the IAU weighted SSH increment
207#if defined key_asminc
208      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
209         CALL ssh_asm_inc( kt )
210         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
211      ENDIF
212#endif
213
214      !                                           !------------------------------!
215      !                                           !     Now Vertical Velocity    !
216      !                                           !------------------------------!
217      !                                                ! integrate from the bottom the hor. divergence
218      DO jk = jpkm1, 1, -1
219         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
220              &                    - (  fse3t_a(:,:,jk)                   &
221              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
222#if defined key_bdy
223         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
224#endif
225      END DO
226      !
227      CALL iom_put( "woce", wn                    )   ! vertical velocity
228      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
229      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
230      IF( lk_diaar5 ) THEN
231         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
232         DO jk = 1, jpk
233            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
234         END DO
235         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
236         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
237      ENDIF
238      !
239   END SUBROUTINE ssh_wzv
240
241
242   SUBROUTINE ssh_nxt( kt )
243      !!----------------------------------------------------------------------
244      !!                    ***  ROUTINE ssh_nxt  ***
245      !!
246      !! ** Purpose :   achieve the sea surface  height time stepping by
247      !!              applying Asselin time filter and swapping the arrays
248      !!              ssha  already computed in ssh_wzv 
249      !!
250      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
251      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
252      !!             sshn = ssha
253      !!
254      !! ** action  : - sshb, sshn   : before & now sea surface height
255      !!                               ready for the next time step
256      !!----------------------------------------------------------------------
257      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
258      !!
259      INTEGER  ::   ji, jj               ! dummy loop indices
260      !!----------------------------------------------------------------------
261
262      IF( kt == nit000 ) THEN
263         IF(lwp) WRITE(numout,*)
264         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
265         IF(lwp) WRITE(numout,*) '~~~~~~~ '
266      ENDIF
267
268      ! Time filter and swap of the ssh
269      ! -------------------------------
270
271      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
272         !                   ! ---------------------- !
273         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
274            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
275            sshu_n(:,:) = sshu_a(:,:)
276            sshv_n(:,:) = sshv_a(:,:)
277            sshf_n(:,:) = sshf_a(:,:)
278         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
279            DO jj = 1, jpj
280               DO ji = 1, jpi                                ! before <-- now filtered
281                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
282                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
283                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
284                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
285                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
286                  sshu_n(ji,jj) = sshu_a(ji,jj)
287                  sshv_n(ji,jj) = sshv_a(ji,jj)
288                  sshf_n(ji,jj) = sshf_a(ji,jj)
289               END DO
290            END DO
291         ENDIF
292         !
293      ELSE                   ! fixed levels :   ssh at t-point only
294         !                   ! ------------ !
295         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
296            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
297            !
298         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
299            DO jj = 1, jpj
300               DO ji = 1, jpi                                ! before <-- now filtered
301                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   
302                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
303               END DO
304            END DO
305         ENDIF
306         !
307      ENDIF
308      !
309      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
310      !
311   END SUBROUTINE ssh_nxt
312
313   !!======================================================================
314END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.