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 #465 – Attachment – NEMO

Ticket #465: sshwzv.F90

File sshwzv.F90, 14.0 KB (added by rachel.furner, 14 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   !!----------------------------------------------------------------------
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 iom
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ssh_wzv    ! called by step.F90
33   PUBLIC   ssh_nxt    ! called by step.F90
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
41   !! $Id: sshwzv.F90 1565 2009-07-31 14:01:08Z rblod $
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE ssh_wzv( kt ) 
48      !!----------------------------------------------------------------------
49      !!                ***  ROUTINE ssh_wzv  ***
50      !!                   
51      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
52      !!              and update the now vertical coordinate (lk_vvl=T).
53      !!
54      !! ** Method  : -
55      !!
56      !!              - Using the incompressibility hypothesis, the vertical
57      !!      velocity is computed by integrating the horizontal divergence 
58      !!      from the bottom to the surface minus the scale factor evolution.
59      !!        The boundary conditions are w=0 at the bottom (no flux) and.
60      !!
61      !! ** action  :   ssha    : after sea surface height
62      !!                wn      : now vertical velocity
63      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
64      !!                          at u-, v-, f-point s
65      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
66      !!----------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   kt   ! time step
68      !!
69      INTEGER  ::   ji, jj, jk      ! dummy loop indices
70      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
71      REAL(wp) ::   z2dt, zraur     ! temporary scalars
72      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
73      !!----------------------------------------------------------------------
74
75      IF( kt == nit000 ) THEN
76         IF(lwp) WRITE(numout,*)
77         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
78         IF(lwp) WRITE(numout,*) '~~~~~~~ '
79         !
80         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
81         !
82         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
83            DO jj = 1, jpjm1
84               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
85                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
86                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
87                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
88                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
89                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
90                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
91                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
92                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
93                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
94                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
95                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
96                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
97                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
98                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
99                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
100               END DO
101            END DO
102            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
103            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
104            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
105         ENDIF
106         !
107      ENDIF
108
109      !                                           !------------------------------!
110      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
111      !                                           !------------------------------!
112         DO jk = 1, jpkm1
113            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
114            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
115            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
116            !
117            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
118            fse3u (:,:,jk) = fse3u_n (:,:,jk)
119            fse3v (:,:,jk) = fse3v_n (:,:,jk)
120            fse3f (:,:,jk) = fse3f_n (:,:,jk)
121            fse3w (:,:,jk) = fse3w_n (:,:,jk)
122            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
123            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
124         END DO
125         !                                             ! now ocean depth (at u- and v-points)
126         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
127         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
128         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
129         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
130         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
131         !
132      ENDIF
133
134                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
135      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
136
137      ! set time step size (Euler/Leapfrog)
138      z2dt = 2. * rdt 
139      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
140
141      zraur = 1. / rauw
142
143      !                                           !------------------------------!
144      !                                           !   After Sea Surface Height   !
145      !                                           !------------------------------!
146      zhdiv(:,:) = 0.e0
147      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
148        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
149      END DO
150
151      !                                                ! Sea surface elevation time stepping
152      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
153
154#if defined key_obc
155# if defined key_agrif
156      IF ( Agrif_Root() ) THEN 
157# endif
158         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
159         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
160# if defined key_agrif
161      ENDIF
162# endif
163#endif
164
165      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
166      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
167         DO jj = 1, jpjm1
168            DO ji = 1, fs_jpim1   ! Vector Opt.
169               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
170                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
171                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
172               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
173                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
174                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
175               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
176                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
177                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
178            END DO
179         END DO
180         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
181         CALL lbc_lnk( sshv_a, 'V', 1. )
182         CALL lbc_lnk( sshf_a, 'F', 1. )
183      ENDIF
184
185      !                                           !------------------------------!
186      !                                           !     Now Vertical Velocity    !
187      !                                           !------------------------------!
188      !                                                ! integrate from the bottom the hor. divergence
189      DO jk = jpkm1, 1, -1
190         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
191              &                    - (  fse3t_a(:,:,jk)                   &
192              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
193      END DO
194      !
195      CALL iom_put( "woce", wn   )                     ! vertical velocity
196      CALL iom_put( "ssh" , sshn )                     ! sea surface height
197      !
198   END SUBROUTINE ssh_wzv
199
200
201   SUBROUTINE ssh_nxt( kt )
202      !!----------------------------------------------------------------------
203      !!                    ***  ROUTINE ssh_nxt  ***
204      !!
205      !! ** Purpose :   achieve the sea surface  height time stepping by
206      !!              applying Asselin time filter and swapping the arrays
207      !!              ssha  already computed in ssh_wzv 
208      !!
209      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
210      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
211      !!             sshn = ssha
212      !!
213      !! ** action  : - sshb, sshn   : before & now sea surface height
214      !!                               ready for the next time step
215      !!----------------------------------------------------------------------
216      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
217      !!
218      INTEGER  ::   ji, jj               ! dummy loop indices
219      !!----------------------------------------------------------------------
220
221      IF( kt == nit000 ) THEN
222         IF(lwp) WRITE(numout,*)
223         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
224         IF(lwp) WRITE(numout,*) '~~~~~~~ '
225      ENDIF
226
227      ! Time filter and swap of the ssh
228      ! -------------------------------
229
230      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
231         !                   ! ---------------------- !
232         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
233            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
234            sshu_n(:,:) = sshu_a(:,:)
235            sshv_n(:,:) = sshv_a(:,:)
236            sshf_n(:,:) = sshf_a(:,:)
237         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
238            DO jj = 1, jpj
239               DO ji = 1, jpi                                ! before <-- now filtered
240                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
241                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
242                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
243                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
244                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
245                  sshu_n(ji,jj) = sshu_a(ji,jj)
246                  sshv_n(ji,jj) = sshv_a(ji,jj)
247                  sshf_n(ji,jj) = sshf_a(ji,jj)
248               END DO
249            END DO
250         ENDIF
251         !
252      ELSE                   ! fixed levels :   ssh at t-point only
253         !                   ! ------------ !
254         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
255            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
256            !
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                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
262               END DO
263            END DO
264         ENDIF
265         !
266      ENDIF
267      !
268      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
269      !
270   END SUBROUTINE ssh_nxt
271
272   !!======================================================================
273END MODULE sshwzv