source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 4153

Last change on this file since 4153 was 4153, checked in by cetlod, 7 years ago

dev_LOCEAN_2013: merge in trunk changes between r3940 and r4028, see ticket #1169

  • 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      !!----------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
114      !
115      !
116      IF( kt == nit000 ) THEN
117         IF(lwp) WRITE(numout,*)
118         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
119         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
120       
121         ! set to zero free surface specific arrays
122         spgu(:,:) = 0._wp                     ! surface pressure gradient (i-direction)
123         spgv(:,:) = 0._wp                     ! surface pressure gradient (j-direction)
124
125         ! read filtered free surface arrays in restart file
126         ! when using agrif, sshn, gcx have to be read in istate
127         IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields:
128         !                                                        ! gcx, gcxb
129      ENDIF
130
131      ! Local constant initialization
132      z2dt = 2. * rdt                                             ! time step: leap-frog
133      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
134      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
135      z2dtg  = grav * z2dt
136
137      ! Evaluate the masked next velocity (effect of the additional force not included)
138      ! --------------------------------- 
139      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
140         !
141         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
142            DO jk = 1, jpkm1
143               DO jj = 2, jpjm1
144                  DO ji = fs_2, fs_jpim1   ! vector opt.
145                     ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)  ) * umask(ji,jj,jk)
146                     va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * va(ji,jj,jk)  ) * vmask(ji,jj,jk)
147                  END DO
148               END DO
149            END DO
150            !
151         ELSE                          ! flux form : applied on thickness weighted velocity
152            DO jk = 1, jpkm1
153               DO jj = 2, jpjm1
154                  DO ji = fs_2, fs_jpim1   ! vector opt.
155                     ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
156                        &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
157                        &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
158                     va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
159                        &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
160                        &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
161                 END DO
162               END DO
163            END DO
164            !
165         ENDIF
166         !
167      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
168         !
169         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
170            DO ji = fs_2, fs_jpim1   ! vector opt.
171               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
172               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
173            END DO
174         END DO
175         DO jk = 1, jpkm1              ! unweighted time stepping
176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.
178                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
179                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
180               END DO
181            END DO
182         END DO
183         !
184      ENDIF
185
186#if defined key_obc
187      IF( lk_obc ) CALL obc_dyn( kt )   ! Update velocities on each open boundary with the radiation algorithm
188      IF( lk_obc ) CALL obc_vol( kt )   ! Correction of the barotropic componant velocity to control the volume of the system
189#endif
190#if defined key_bdy
191      IF( lk_bdy ) CALL bdy_dyn( kt )   ! Update velocities on each open boundary
192      IF( lk_bdy ) CALL bdy_vol( kt )   ! Correction of the barotropic component velocity to control the volume of the system
193#endif
194#if defined key_agrif
195      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
196#endif
197      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
198
199      ! compute the next vertically averaged velocity (effect of the additional force not included)
200      ! ---------------------------------------------
201      DO jj = 2, jpjm1
202         DO ji = fs_2, fs_jpim1   ! vector opt.
203            spgu(ji,jj) = 0._wp
204            spgv(ji,jj) = 0._wp
205         END DO
206      END DO
207
208      ! vertical sum
209!CDIR NOLOOPCHG
210      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
211         DO jk = 1, jpkm1
212            DO ji = 1, jpij
213               spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk)
214               spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk)
215            END DO
216         END DO
217      ELSE                        ! No  vector opt.
218         DO jk = 1, jpkm1
219            DO jj = 2, jpjm1
220               DO ji = 2, jpim1
221                  spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
222                  spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
223               END DO
224            END DO
225         END DO
226      ENDIF
227
228      ! transport: multiplied by the horizontal scale factor
229      DO jj = 2, jpjm1
230         DO ji = fs_2, fs_jpim1   ! vector opt.
231            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
232            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
233         END DO
234      END DO
235      CALL lbc_lnk( spgu, 'U', -1. )       ! lateral boundary conditions
236      CALL lbc_lnk( spgv, 'V', -1. )
237
238      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only)
239
240      ! Right hand side of the elliptic equation and first guess
241      ! --------------------------------------------------------
242      DO jj = 2, jpjm1
243         DO ji = fs_2, fs_jpim1   ! vector opt.
244            ! Divergence of the after vertically averaged velocity
245            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
246                  + spgv(ji,jj) - spgv(ji,jj-1)
247            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
248            ! First guess of the after barotropic transport divergence
249            zbtd = gcx(ji,jj)
250            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
251            gcxb(ji,jj) =      zbtd
252         END DO
253      END DO
254      ! applied the lateral boundary conditions
255      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
256
257#if defined key_agrif
258      IF( .NOT. AGRIF_ROOT() ) THEN
259         ! add contribution of gradient of after barotropic transport divergence
260         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
261            &    gcb(3     ,:) - z2dtg * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
262         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
263            &    gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
264         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
265            &    gcb(:,3     ) - z2dtg * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
266         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
267            &    gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
268      ENDIF
269#endif
270
271
272      ! Relative precision (computation on one processor)
273      ! ------------------
274      rnorme =0.e0
275      rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
276
277      epsr = eps * eps * rnorme
278      ncut = 0
279      ! if rnorme is 0, the solution is 0, the solver is not called
280      IF( rnorme == 0._wp ) THEN
281         gcx(:,:) = 0._wp
282         res   = 0._wp
283         niter = 0
284         ncut  = 999
285      ENDIF
286
287      ! Evaluate the next transport divergence
288      ! --------------------------------------
289      !    Iterarive solver for the elliptic equation (except IF sol.=0)
290      !    (output in gcx with boundary conditions applied)
291      kindic = 0
292      IF( ncut == 0 ) THEN
293         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
294         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
295         ENDIF
296      ENDIF
297
298      ! Transport divergence gradient multiplied by z2dt
299      ! --------------------------------------------====
300      DO jj = 2, jpjm1
301         DO ji = fs_2, fs_jpim1   ! vector opt.
302            ! trend of Transport divergence gradient
303            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
304            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
305            ! multiplied by z2dt
306#if defined key_obc
307            IF(lk_obc) THEN
308            ! caution : grad D = 0 along open boundaries
309            ! Remark: The filtering force could be reduced here in the FRS zone
310            !         by multiplying spgu/spgv by (1-alpha) ?? 
311               spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
312               spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
313            ELSE
314               spgu(ji,jj) = z2dt * ztdgu
315               spgv(ji,jj) = z2dt * ztdgv
316            ENDIF
317#elif defined key_bdy
318            IF(lk_bdy) THEN
319            ! caution : grad D = 0 along open boundaries
320               spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
321               spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
322            ELSE
323               spgu(ji,jj) = z2dt * ztdgu
324               spgv(ji,jj) = z2dt * ztdgv
325            ENDIF
326#else
327            spgu(ji,jj) = z2dt * ztdgu
328            spgv(ji,jj) = z2dt * ztdgv
329#endif
330         END DO
331      END DO
332
333#if defined key_agrif     
334      IF( .NOT. Agrif_Root() ) THEN
335         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
336         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
337         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
338         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
339         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
340      ENDIF
341#endif     
342      ! Add the trends multiplied by z2dt to the after velocity
343      ! -------------------------------------------------------
344      !     ( c a u t i o n : (ua,va) here are the after velocity not the
345      !                       trend, the leap-frog time stepping will not
346      !                       be done in dynnxt.F90 routine)
347      DO jk = 1, jpkm1
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
351               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
352            END DO
353         END DO
354      END DO
355
356      ! write filtered free surface arrays in restart file
357      ! --------------------------------------------------
358      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
359      !
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
362      !
363   END SUBROUTINE dyn_spg_flt
364
365
366   SUBROUTINE flt_rst( kt, cdrw )
367      !!---------------------------------------------------------------------
368      !!                   ***  ROUTINE ts_rst  ***
369      !!
370      !! ** Purpose : Read or write filtered free surface arrays in restart file
371      !!----------------------------------------------------------------------
372      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
373      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
374      !!----------------------------------------------------------------------
375      !
376      IF( TRIM(cdrw) == 'READ' ) THEN
377         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
378! Caution : extra-hallow
379! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
380            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
381            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
382            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
383         ELSE
384            gcx (:,:) = 0.e0
385            gcxb(:,:) = 0.e0
386         ENDIF
387      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
388! Caution : extra-hallow
389! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
390         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
391         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
392      ENDIF
393      !
394   END SUBROUTINE flt_rst
395
396#else
397   !!----------------------------------------------------------------------
398   !!   Default case :   Empty module   No standart free surface cst volume
399   !!----------------------------------------------------------------------
400CONTAINS
401   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
402      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
403   END SUBROUTINE dyn_spg_flt
404   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
405      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
406      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
407      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
408   END SUBROUTINE flt_rst
409#endif
410   
411   !!======================================================================
412END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.