source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 3765

Last change on this file since 3765 was 3765, checked in by smasson, 8 years ago

dev_MERGE_2012: add missing bugfix in r3692 see ticket #1029

  • Property svn:keywords set to Id
File size: 19.3 KB
Line 
1MODULE dynspg_flt
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_flt  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6   !! History    OPA  !  1998-05  (G. Roullet)  free surface
7   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2
8   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module
9   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
10   !!            1.0  !  2004-08  (C. Talandier) New trends organization
11   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
12   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
13   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.
14   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
15   !!----------------------------------------------------------------------
16#if defined key_dynspg_flt   ||   defined key_esopa 
17   !!----------------------------------------------------------------------
18   !!   'key_dynspg_flt'                              filtered free surface
19   !!----------------------------------------------------------------------
20   !!   dyn_spg_flt  : update the momentum trend with the surface pressure gradient in the filtered free surface case
21   !!   flt_rst      : read/write the time-splitting restart fields in the ocean restart file
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE zdf_oce         ! ocean vertical physics
26   USE sbc_oce         ! surface boundary condition: ocean
27   USE obc_oce         ! Lateral open boundary condition
28   USE bdy_oce         ! Lateral open boundary condition
29   USE sol_oce         ! ocean elliptic solver
30   USE phycst          ! physical constants
31   USE domvvl          ! variable volume
32   USE dynadv          ! advection
33   USE solmat          ! matrix construction for elliptic solvers
34   USE solpcg          ! preconditionned conjugate gradient solver
35   USE solsor          ! Successive Over-relaxation solver
36   USE obcdyn          ! ocean open boundary condition on dynamics
37   USE obcvol          ! ocean open boundary condition (obc_vol routine)
38   USE bdydyn          ! ocean open boundary condition on dynamics
39   USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
40   USE cla             ! cross land advection
41   USE in_out_manager  ! I/O manager
42   USE lib_mpp         ! distributed memory computing library
43   USE wrk_nemo        ! Memory Allocation
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl          ! Print control
46   USE iom
47   USE lib_fortran
48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
51   USE timing          ! Timing
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC   dyn_spg_flt  ! routine called by step.F90
57   PUBLIC   flt_rst      ! routine called by istate.F90
58
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
64   !! $Id$
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE dyn_spg_flt( kt, kindic )
70      !!----------------------------------------------------------------------
71      !!                  ***  routine dyn_spg_flt  ***
72      !!
73      !! ** Purpose :   Compute the now trend due to the surface pressure
74      !!      gradient in case of filtered free surface formulation  and add
75      !!      it to the general trend of momentum equation.
76      !!
77      !! ** Method  :   Filtered free surface formulation. The surface
78      !!      pressure gradient is given by:
79      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
80      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
81      !!      where sshn is the free surface elevation and btda is the after
82      !!      time derivative of the free surface elevation
83      !!       -1- evaluate the surface presure trend (including the addi-
84      !!      tional force) in three steps:
85      !!        a- compute the right hand side of the elliptic equation:
86      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
87      !!         where (spgu,spgv) are given by:
88      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
89      !!                 - grav 2 rdt hu /e1u di[sshn + (emp-rnf)]
90      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
91      !!                 - grav 2 rdt hv /e2v dj[sshn + (emp-rnf)]
92      !!         and define the first guess from previous computation :
93      !!            zbtd = btda
94      !!            btda = 2 zbtd - btdb
95      !!            btdb = zbtd
96      !!        b- compute the relative accuracy to be reached by the
97      !!         iterative solver
98      !!        c- apply the solver by a call to sol... routine
99      !!       -2- compute and add the free surface pressure gradient inclu-
100      !!      ding the additional force used to stabilize the equation.
101      !!
102      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
103      !!
104      !! References : Roullet and Madec 1999, JGR.
105      !!---------------------------------------------------------------------
106      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
107      INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge)
108      !!                                   
109      INTEGER  ::   ji, jj, jk   ! dummy loop indices
110      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars
111      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zub, zvb
112      !!----------------------------------------------------------------------
113      !
114      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
115      !
116      CALL wrk_alloc( jpi,jpj,jpk, zub, zvb )
117      !
118      IF( kt == nit000 ) THEN
119         IF(lwp) WRITE(numout,*)
120         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
121         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
122       
123         ! set to zero free surface specific arrays
124         spgu(:,:) = 0._wp                     ! surface pressure gradient (i-direction)
125         spgv(:,:) = 0._wp                     ! surface pressure gradient (j-direction)
126
127         ! read filtered free surface arrays in restart file
128         ! when using agrif, sshn, gcx have to be read in istate
129         IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields:
130         !                                                        ! gcx, gcxb
131      ENDIF
132
133      ! Local constant initialization
134      z2dt = 2. * rdt                                             ! time step: leap-frog
135      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
136      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
137      z2dtg  = grav * z2dt
138
139      ! Evaluate the masked next velocity (effect of the additional force not included)
140      ! --------------------------------- 
141      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
142         !
143         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
144            DO jk = 1, jpkm1
145               DO jj = 2, jpjm1
146                  DO ji = fs_2, fs_jpim1   ! vector opt.
147                     ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)  ) * umask(ji,jj,jk)
148                     va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * va(ji,jj,jk)  ) * vmask(ji,jj,jk)
149                  END DO
150               END DO
151            END DO
152            !
153         ELSE                          ! flux form : applied on thickness weighted velocity
154            DO jk = 1, jpkm1
155               DO jj = 2, jpjm1
156                  DO ji = fs_2, fs_jpim1   ! vector opt.
157                     ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
158                        &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
159                        &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
160                     va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
161                        &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
162                        &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
163                 END DO
164               END DO
165            END DO
166            !
167         ENDIF
168         !
169      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
170         !
171         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
174               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
175            END DO
176         END DO
177         DO jk = 1, jpkm1              ! unweighted time stepping
178            DO jj = 2, jpjm1
179               DO ji = fs_2, fs_jpim1   ! vector opt.
180                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
181                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
182               END DO
183            END DO
184         END DO
185         !
186      ENDIF
187
188#if defined key_obc
189      IF( lk_obc ) CALL obc_dyn( kt )   ! Update velocities on each open boundary with the radiation algorithm
190      IF( lk_obc ) CALL obc_vol( kt )   ! Correction of the barotropic componant velocity to control the volume of the system
191#endif
192#if defined key_bdy
193      IF( lk_bdy ) CALL bdy_dyn( kt )   ! Update velocities on each open boundary
194      IF( lk_bdy ) CALL bdy_vol( kt )   ! Correction of the barotropic component velocity to control the volume of the system
195#endif
196#if defined key_agrif
197      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
198#endif
199      IF( nn_cla == 1 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
200
201      ! compute the next vertically averaged velocity (effect of the additional force not included)
202      ! ---------------------------------------------
203      DO jj = 2, jpjm1
204         DO ji = fs_2, fs_jpim1   ! vector opt.
205            spgu(ji,jj) = 0._wp
206            spgv(ji,jj) = 0._wp
207         END DO
208      END DO
209
210      ! vertical sum
211!CDIR NOLOOPCHG
212      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
213         DO jk = 1, jpkm1
214            DO ji = 1, jpij
215               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
216               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
217            END DO
218         END DO
219      ELSE                        ! No  vector opt.
220         DO jk = 1, jpkm1
221            DO jj = 2, jpjm1
222               DO ji = 2, jpim1
223                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
224                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
225               END DO
226            END DO
227         END DO
228      ENDIF
229
230      ! transport: multiplied by the horizontal scale factor
231      DO jj = 2, jpjm1
232         DO ji = fs_2, fs_jpim1   ! vector opt.
233            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
234            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
235         END DO
236      END DO
237      CALL lbc_lnk( spgu, 'U', -1. )       ! lateral boundary conditions
238      CALL lbc_lnk( spgv, 'V', -1. )
239
240      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only)
241
242      ! Right hand side of the elliptic equation and first guess
243      ! --------------------------------------------------------
244      DO jj = 2, jpjm1
245         DO ji = fs_2, fs_jpim1   ! vector opt.
246            ! Divergence of the after vertically averaged velocity
247            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
248                  + spgv(ji,jj) - spgv(ji,jj-1)
249            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
250            ! First guess of the after barotropic transport divergence
251            zbtd = gcx(ji,jj)
252            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
253            gcxb(ji,jj) =      zbtd
254         END DO
255      END DO
256      ! applied the lateral boundary conditions
257      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
258
259#if defined key_agrif
260      IF( .NOT. AGRIF_ROOT() ) THEN
261         ! add contribution of gradient of after barotropic transport divergence
262         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
263            &    gcb(3     ,:) - z2dtg * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
264         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
265            &    gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
266         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
267            &    gcb(:,3     ) - z2dtg * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
268         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
269            &    gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
270      ENDIF
271#endif
272
273
274      ! Relative precision (computation on one processor)
275      ! ------------------
276      rnorme =0.e0
277      rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
278
279      epsr = eps * eps * rnorme
280      ncut = 0
281      ! if rnorme is 0, the solution is 0, the solver is not called
282      IF( rnorme == 0._wp ) THEN
283         gcx(:,:) = 0._wp
284         res   = 0._wp
285         niter = 0
286         ncut  = 999
287      ENDIF
288
289      ! Evaluate the next transport divergence
290      ! --------------------------------------
291      !    Iterarive solver for the elliptic equation (except IF sol.=0)
292      !    (output in gcx with boundary conditions applied)
293      kindic = 0
294      IF( ncut == 0 ) THEN
295         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
296         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
297         ENDIF
298      ENDIF
299
300      ! Transport divergence gradient multiplied by z2dt
301      ! --------------------------------------------====
302      DO jj = 2, jpjm1
303         DO ji = fs_2, fs_jpim1   ! vector opt.
304            ! trend of Transport divergence gradient
305            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
306            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
307            ! multiplied by z2dt
308#if defined key_obc
309            IF(lk_obc) THEN
310            ! caution : grad D = 0 along open boundaries
311            ! Remark: The filtering force could be reduced here in the FRS zone
312            !         by multiplying spgu/spgv by (1-alpha) ?? 
313               spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
314               spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
315            ELSE
316               spgu(ji,jj) = z2dt * ztdgu
317               spgv(ji,jj) = z2dt * ztdgv
318            ENDIF
319#elif defined key_bdy
320            IF(lk_bdy) THEN
321            ! caution : grad D = 0 along open boundaries
322               spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
323               spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
324            ELSE
325               spgu(ji,jj) = z2dt * ztdgu
326               spgv(ji,jj) = z2dt * ztdgv
327            ENDIF
328#else
329            spgu(ji,jj) = z2dt * ztdgu
330            spgv(ji,jj) = z2dt * ztdgv
331#endif
332         END DO
333      END DO
334
335#if defined key_agrif     
336      IF( .NOT. Agrif_Root() ) THEN
337         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
338         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
339         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
340         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
341         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
342      ENDIF
343#endif     
344      ! Add the trends multiplied by z2dt to the after velocity
345      ! -------------------------------------------------------
346      !     ( c a u t i o n : (ua,va) here are the after velocity not the
347      !                       trend, the leap-frog time stepping will not
348      !                       be done in dynnxt.F90 routine)
349      DO jk = 1, jpkm1
350         DO jj = 2, jpjm1
351            DO ji = fs_2, fs_jpim1   ! vector opt.
352               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
353               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
354            END DO
355         END DO
356      END DO
357
358      ! write filtered free surface arrays in restart file
359      ! --------------------------------------------------
360      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
361      !
362      CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb )
363      !
364      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
365      !
366   END SUBROUTINE dyn_spg_flt
367
368
369   SUBROUTINE flt_rst( kt, cdrw )
370      !!---------------------------------------------------------------------
371      !!                   ***  ROUTINE ts_rst  ***
372      !!
373      !! ** Purpose : Read or write filtered free surface arrays in restart file
374      !!----------------------------------------------------------------------
375      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
376      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
377      !!----------------------------------------------------------------------
378      !
379      IF( TRIM(cdrw) == 'READ' ) THEN
380         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
381! Caution : extra-hallow
382! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
383            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
384            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
385            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
386         ELSE
387            gcx (:,:) = 0.e0
388            gcxb(:,:) = 0.e0
389         ENDIF
390      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
391! Caution : extra-hallow
392! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
393         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
394         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
395      ENDIF
396      !
397   END SUBROUTINE flt_rst
398
399#else
400   !!----------------------------------------------------------------------
401   !!   Default case :   Empty module   No standart free surface cst volume
402   !!----------------------------------------------------------------------
403CONTAINS
404   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
405      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
406   END SUBROUTINE dyn_spg_flt
407   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
408      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
409      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
410      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
411   END SUBROUTINE flt_rst
412#endif
413   
414   !!======================================================================
415END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.