source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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