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

source: branches/DEV_R1821_Rivers/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2101

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

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