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 on 2009WP/2009Stream3/VVL – Attachment – NEMO

2009WP/2009Stream3/VVL: dynspg_flt.F90

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