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.
dynspg_flt.F90 in branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

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

Last change on this file since 3764 was 3764, checked in by smasson, 11 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012

  • 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 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            ! caution : grad D = 0 along open boundaries
310            ! Remark: The filtering force could be reduced here in the FRS zone
311            !         by multiplying spgu/spgv by (1-alpha) ?? 
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            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
317            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
318#else
319            spgu(ji,jj) = z2dt * ztdgu
320            spgv(ji,jj) = z2dt * ztdgv
321#endif
322         END DO
323      END DO
324
325#if defined key_agrif     
326      IF( .NOT. Agrif_Root() ) THEN
327         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
328         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
329         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
330         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
331         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
332      ENDIF
333#endif     
334      ! Add the trends multiplied by z2dt to the after velocity
335      ! -------------------------------------------------------
336      !     ( c a u t i o n : (ua,va) here are the after velocity not the
337      !                       trend, the leap-frog time stepping will not
338      !                       be done in dynnxt.F90 routine)
339      DO jk = 1, jpkm1
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1   ! vector opt.
342               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
343               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
344            END DO
345         END DO
346      END DO
347
348      ! write filtered free surface arrays in restart file
349      ! --------------------------------------------------
350      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
351      !
352      CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb )
353      !
354      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
355      !
356   END SUBROUTINE dyn_spg_flt
357
358
359   SUBROUTINE flt_rst( kt, cdrw )
360      !!---------------------------------------------------------------------
361      !!                   ***  ROUTINE ts_rst  ***
362      !!
363      !! ** Purpose : Read or write filtered free surface arrays in restart file
364      !!----------------------------------------------------------------------
365      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
366      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
367      !!----------------------------------------------------------------------
368      !
369      IF( TRIM(cdrw) == 'READ' ) THEN
370         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
371! Caution : extra-hallow
372! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
373            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
374            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
375            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
376         ELSE
377            gcx (:,:) = 0.e0
378            gcxb(:,:) = 0.e0
379         ENDIF
380      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
381! Caution : extra-hallow
382! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
383         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
384         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
385      ENDIF
386      !
387   END SUBROUTINE flt_rst
388
389#else
390   !!----------------------------------------------------------------------
391   !!   Default case :   Empty module   No standart free surface cst volume
392   !!----------------------------------------------------------------------
393CONTAINS
394   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
395      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
396   END SUBROUTINE dyn_spg_flt
397   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
398      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
399      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
400      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
401   END SUBROUTINE flt_rst
402#endif
403   
404   !!======================================================================
405END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.