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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta
  • Property svn:keywords set to Id
File size: 18.1 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 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  : h_rnf, nk_rnf  ! River runoff
33#if defined key_agrif
34   USE agrif_opa_update
35   USE agrif_opa_interp
36#endif
37#if defined key_asminc   
38   USE asminc          ! Assimilation increment
39#endif
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   ssh_wzv    ! called by step.F90
45   PUBLIC   ssh_nxt    ! called by step.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE ssh_wzv( kt ) 
59      !!----------------------------------------------------------------------
60      !!                ***  ROUTINE ssh_wzv  ***
61      !!                   
62      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
63      !!              and update the now vertical coordinate (lk_vvl=T).
64      !!
65      !! ** Method  : - Using the incompressibility hypothesis, the vertical
66      !!      velocity is computed by integrating the horizontal divergence 
67      !!      from the bottom to the surface minus the scale factor evolution.
68      !!        The boundary conditions are w=0 at the bottom (no flux) and.
69      !!
70      !! ** action  :   ssha    : after sea surface height
71      !!                wn      : now vertical velocity
72      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
73      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
74      !!
75      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
76      !!----------------------------------------------------------------------
77      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
78      !!
79      INTEGER, INTENT(in) ::   kt   ! time step
80      !!
81      INTEGER  ::   ji, jj, jk      ! dummy loop indices
82      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
83      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
84      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
85      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
86      !!----------------------------------------------------------------------
87
88      IF( kt == nit000 ) THEN
89         !
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
92         IF(lwp) WRITE(numout,*) '~~~~~~~ '
93         !
94         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all)
95         !
96         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
97            DO jj = 1, jpjm1
98               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
99                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
100                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
101                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
102                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
103                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
104                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
105                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
106                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
107                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
108                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
109                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
110               END DO
111            END DO
112            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
113            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
114            DO jj = 1, jpjm1
115               DO ji = 1, jpim1      ! NO Vector Opt.
116                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
117                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
118                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
119                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
120               END DO
121            END DO
122            CALL lbc_lnk( sshf_n, 'F', 1. )
123         ENDIF
124         !
125      ENDIF
126
127      !                                           !------------------------------------------!
128      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
129         !                                        !------------------------------------------!
130         DO jk = 1, jpkm1
131            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
132            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
133            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
134            !
135            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
136            fse3u (:,:,jk) = fse3u_n (:,:,jk)
137            fse3v (:,:,jk) = fse3v_n (:,:,jk)
138            fse3f (:,:,jk) = fse3f_n (:,:,jk)
139            fse3w (:,:,jk) = fse3w_n (:,:,jk)
140            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
141            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
142         END DO
143         !
144         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
145         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
146         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
147         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
148         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
149         !
150      ENDIF
151      !
152      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
153      !
154      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
155      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
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      !                                                ! Sea surface elevation time stepping
165      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
166      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
167      z1_rau0 = 0.5 / rau0
168      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) &
169      &                      * tmask(:,:,1)
170
171#if defined key_agrif
172      CALL agrif_ssh(kt)
173#endif
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#if defined key_bdy
181      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
182      CALL lbc_lnk( ssha, 'T', 1. ) 
183#endif
184
185      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
186      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
187         DO jj = 1, jpjm1
188            DO ji = 1, jpim1      ! NO Vector Opt.
189               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
190                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
191                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
192               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
193                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
194                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
195            END DO
196         END DO
197         ! Boundaries conditions
198         CALL lbc_lnk( sshu_a, 'U', 1. )
199         CALL lbc_lnk( sshv_a, 'V', 1. )
200      ENDIF
201! Include the IAU weighted SSH increment
202#if defined key_asminc
203      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
204         CALL ssh_asm_inc( kt )
205         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
206      ENDIF
207#endif
208
209      !                                           !------------------------------!
210      !                                           !     Now Vertical Velocity    !
211      !                                           !------------------------------!
212      z1_2dt = 1.e0 / z2dt
213      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
214         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
215         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
216            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
217            &                         * tmask(:,:,jk) * z1_2dt
218#if defined key_bdy
219         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
220#endif
221      END DO
222
223      !                                           !------------------------------!
224      !                                           !           outputs            !
225      !                                           !------------------------------!
226      CALL iom_put( "woce", wn                    )   ! vertical velocity
227      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
228      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
229      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
230         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
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                     ) 
236         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
237      ENDIF
238      !
239      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
240      !
241   END SUBROUTINE ssh_wzv
242
243
244   SUBROUTINE ssh_nxt( kt )
245      !!----------------------------------------------------------------------
246      !!                    ***  ROUTINE ssh_nxt  ***
247      !!
248      !! ** Purpose :   achieve the sea surface  height time stepping by
249      !!              applying Asselin time filter and swapping the arrays
250      !!              ssha  already computed in ssh_wzv 
251      !!
252      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
253      !!              from the filter, see Leclair and Madec 2010) and swap :
254      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
255      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
256      !!                sshn = ssha
257      !!
258      !! ** action  : - sshb, sshn   : before & now sea surface height
259      !!                               ready for the next time step
260      !!
261      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
262      !!----------------------------------------------------------------------
263      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
264      !!
265      INTEGER  ::   ji, jj   ! dummy loop indices
266      REAL(wp) ::   zec      ! temporary scalar
267      !!----------------------------------------------------------------------
268
269      IF( kt == nit000 ) THEN
270         IF(lwp) WRITE(numout,*)
271         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
272         IF(lwp) WRITE(numout,*) '~~~~~~~ '
273      ENDIF
274
275      !                       !--------------------------!
276      IF( lk_vvl ) THEN       !  Variable volume levels  !
277         !                    !--------------------------!
278         !
279         ! ssh at t-, u-, v, f-points
280         !===========================
281         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
282            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
283            sshu_n(:,:) = sshu_a(:,:)
284            sshv_n(:,:) = sshv_a(:,:)
285            DO jj = 1, jpjm1
286               DO ji = 1, jpim1      ! NO Vector Opt.
287                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
288                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
289                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
290                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
291               END DO
292            END DO
293            ! Boundaries conditions
294            CALL lbc_lnk( sshf_n, 'F', 1. )
295         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
296            zec = atfp * rdt / rau0
297            DO jj = 1, jpj
298               DO ji = 1, jpi                                ! before <-- now filtered
299                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
300                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
301                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
302                  sshu_n(ji,jj) = sshu_a(ji,jj)
303                  sshv_n(ji,jj) = sshv_a(ji,jj)
304               END DO
305            END DO
306            DO jj = 1, jpjm1
307               DO ji = 1, jpim1      ! NO Vector Opt.
308                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
309                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
310                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
311                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
312               END DO
313            END DO
314            ! Boundaries conditions
315            CALL lbc_lnk( sshf_n, 'F', 1. )
316            DO jj = 1, jpjm1
317               DO ji = 1, jpim1      ! NO Vector Opt.
318                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
319                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
320                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
321                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
322                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
323                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
324               END DO
325            END DO
326            ! Boundaries conditions
327            CALL lbc_lnk( sshu_b, 'U', 1. )
328            CALL lbc_lnk( sshv_b, 'V', 1. )
329         ENDIF
330         !                    !--------------------------!
331      ELSE                    !        fixed levels      !
332         !                    !--------------------------!
333         !
334         ! ssh at t-point only
335         !====================
336         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
337            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
338            !
339         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
340            DO jj = 1, jpj
341               DO ji = 1, jpi                                ! before <-- now filtered
342                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
343                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
344               END DO
345            END DO
346         ENDIF
347         !
348      ENDIF
349      !
350      ! Update velocity at AGRIF zoom boundaries
351#if defined key_agrif
352      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
353#endif
354      !
355      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
356      !
357   END SUBROUTINE ssh_nxt
358
359   !!======================================================================
360END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.