source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 10 years ago

update licence of all NEMO files…

  • Property svn:keywords set to Id
File size: 18.0 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   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ssh_wzv        : after ssh & now vertical velocity
14   !!   ssh_nxt        : filter ans swap the ssh arrays
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE domvvl          ! Variable volume
20   USE divcur          ! hor. divergence and curl      (div & cur routines)
21   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
22   USE iom             ! I/O library
23   USE restart         ! only for lrst_oce
24   USE in_out_manager  ! I/O manager
25   USE prtctl          ! Print control
26   USE phycst
27   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
28   USE obc_par         ! open boundary cond. parameter
29   USE obc_oce
30   USE bdy_oce
31   USE diaar5, ONLY :   lk_diaar5
32   USE iom
33   USE sbcrnf, ONLY  : h_rnf, nk_rnf  ! River runoff
34#if defined key_asminc   
35   USE asminc          ! Assimilation increment
36#endif
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   ssh_wzv    ! called by step.F90
42   PUBLIC   ssh_nxt    ! called by step.F90
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence (NEMOGCM/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  : - Using the incompressibility hypothesis, the vertical
63      !!      velocity is computed by integrating the horizontal divergence 
64      !!      from the bottom to the surface minus the scale factor evolution.
65      !!        The boundary conditions are w=0 at the bottom (no flux) and.
66      !!
67      !! ** action  :   ssha    : after sea surface height
68      !!                wn      : now vertical velocity
69      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
70      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
71      !!
72      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
73      !!----------------------------------------------------------------------
74      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
75      !!
76      INTEGER, INTENT(in) ::   kt   ! time step
77      !!
78      INTEGER  ::   ji, jj, jk      ! dummy loop indices
79      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
80      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
81      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
82      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
83      !!----------------------------------------------------------------------
84
85      IF( kt == nit000 ) THEN
86         !
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                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
104                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
105                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
106                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
107               END DO
108            END DO
109            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
110            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
111            DO jj = 1, jpjm1
112               DO ji = 1, jpim1      ! NO Vector Opt.
113                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
114                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
115                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
116                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
117               END DO
118            END DO
119            CALL lbc_lnk( sshf_n, 'F', 1. )
120         ENDIF
121         !
122      ENDIF
123
124      !                                           !------------------------------------------!
125      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
126         !                                        !------------------------------------------!
127         DO jk = 1, jpkm1
128            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
129            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
130            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
131            !
132            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
133            fse3u (:,:,jk) = fse3u_n (:,:,jk)
134            fse3v (:,:,jk) = fse3v_n (:,:,jk)
135            fse3f (:,:,jk) = fse3f_n (:,:,jk)
136            fse3w (:,:,jk) = fse3w_n (:,:,jk)
137            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
138            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
139         END DO
140         !
141         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
142         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
143         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
144         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
145         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
146         !
147      ENDIF
148      !
149
150                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
151      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
152
153      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
154      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
155
156      !                                           !------------------------------!
157      !                                           !   After Sea Surface Height   !
158      !                                           !------------------------------!
159      zhdiv(:,:) = 0.e0
160      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
161        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
162      END DO
163      !                                                ! Sea surface elevation time stepping
164      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
165      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
166      z1_rau0 = 0.5 / rau0
167      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) &
168      &                      * tmask(:,:,1)
169
170#if defined key_obc
171      IF( Agrif_Root() ) THEN
172         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
173         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
174      ENDIF
175#endif
176#if defined key_bdy
177      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
178      CALL lbc_lnk( ssha, 'T', 1. ) 
179#endif
180
181      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
182      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
183         DO jj = 1, jpjm1
184            DO ji = 1, jpim1      ! NO Vector Opt.
185               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
186                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
187                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
188               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
189                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
190                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
191            END DO
192         END DO
193         ! Boundaries conditions
194         CALL lbc_lnk( sshu_a, 'U', 1. )
195         CALL lbc_lnk( sshv_a, 'V', 1. )
196      ENDIF
197! Include the IAU weighted SSH increment
198#if defined key_asminc
199      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
200         CALL ssh_asm_inc( kt )
201         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
202      ENDIF
203#endif
204
205      !                                           !------------------------------!
206      !                                           !     Now Vertical Velocity    !
207      !                                           !------------------------------!
208      z1_2dt = 1.e0 / z2dt
209      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
210         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
211         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
212            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
213            &                         * tmask(:,:,jk) * z1_2dt
214#if defined key_bdy
215         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
216#endif
217      END DO
218
219      !                                           !------------------------------!
220      !                                           !           outputs            !
221      !                                           !------------------------------!
222      CALL iom_put( "woce", wn                    )   ! vertical velocity
223      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
224      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
225      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
226         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
227         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
228         DO jk = 1, jpk
229            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
230         END DO
231         CALL iom_put( "w_masstr" , z3d                     ) 
232         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
233      ENDIF
234      !
235      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
236      !
237   END SUBROUTINE ssh_wzv
238
239
240   SUBROUTINE ssh_nxt( kt )
241      !!----------------------------------------------------------------------
242      !!                    ***  ROUTINE ssh_nxt  ***
243      !!
244      !! ** Purpose :   achieve the sea surface  height time stepping by
245      !!              applying Asselin time filter and swapping the arrays
246      !!              ssha  already computed in ssh_wzv 
247      !!
248      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
249      !!              from the filter, see Leclair and Madec 2010) and swap :
250      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
251      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
252      !!                sshn = ssha
253      !!
254      !! ** action  : - sshb, sshn   : before & now sea surface height
255      !!                               ready for the next time step
256      !!
257      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
258      !!----------------------------------------------------------------------
259      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
260      !!
261      INTEGER  ::   ji, jj   ! dummy loop indices
262      REAL(wp) ::   zec      ! temporary scalar
263      !!----------------------------------------------------------------------
264
265      IF( kt == nit000 ) THEN
266         IF(lwp) WRITE(numout,*)
267         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
268         IF(lwp) WRITE(numout,*) '~~~~~~~ '
269      ENDIF
270
271      !                       !--------------------------!
272      IF( lk_vvl ) THEN       !  Variable volume levels  !
273         !                    !--------------------------!
274         !
275         ! ssh at t-, u-, v, f-points
276         !===========================
277         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
278            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
279            sshu_n(:,:) = sshu_a(:,:)
280            sshv_n(:,:) = sshv_a(:,:)
281            DO jj = 1, jpjm1
282               DO ji = 1, jpim1      ! NO Vector Opt.
283                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
284                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
285                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
286                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
287               END DO
288            END DO
289            ! Boundaries conditions
290            CALL lbc_lnk( sshf_n, 'F', 1. )
291         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
292            zec = atfp * rdt / rau0
293            DO jj = 1, jpj
294               DO ji = 1, jpi                                ! before <-- now filtered
295                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
296                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
297                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
298                  sshu_n(ji,jj) = sshu_a(ji,jj)
299                  sshv_n(ji,jj) = sshv_a(ji,jj)
300               END DO
301            END DO
302            DO jj = 1, jpjm1
303               DO ji = 1, jpim1      ! NO Vector Opt.
304                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
305                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
306                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
307                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
308               END DO
309            END DO
310            ! Boundaries conditions
311            CALL lbc_lnk( sshf_n, 'F', 1. )
312            DO jj = 1, jpjm1
313               DO ji = 1, jpim1      ! NO Vector Opt.
314                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
315                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
316                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
317                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
318                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
319                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
320               END DO
321            END DO
322            ! Boundaries conditions
323            CALL lbc_lnk( sshu_b, 'U', 1. )
324            CALL lbc_lnk( sshv_b, 'V', 1. )
325         ENDIF
326         !                    !--------------------------!
327      ELSE                    !        fixed levels      !
328         !                    !--------------------------!
329         !
330         ! ssh at t-point only
331         !====================
332         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
333            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
334            !
335         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
336            DO jj = 1, jpj
337               DO ji = 1, jpi                                ! before <-- now filtered
338                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
339                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
340               END DO
341            END DO
342         ENDIF
343         !
344      ENDIF
345      !
346      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
347      !
348   END SUBROUTINE ssh_nxt
349
350   !!======================================================================
351END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.