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 trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 1662

Last change on this file since 1662 was 1601, checked in by ctlod, 15 years ago

Doctor naming of OPA namelist variables , see ticket: #526

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 18.6 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 solmat          ! matrix construction for elliptic solvers
32   USE solver          ! solver initialization
33   USE solpcg          ! preconditionned conjugate gradient solver
34   USE solsor          ! Successive Over-relaxation solver
35   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
36   USE obcvol          ! ocean open boundary condition (obc_vol routines)
37   USE bdy_oce         ! Unstructured open boundaries condition
38   USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine)
39   USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine)
40   USE cla_dynspg      ! cross land advection
41   USE in_out_manager  ! I/O manager
42   USE lib_mpp         ! distributed memory computing library
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE agrif_opa_interp
46   USE iom
47   USE restart         ! only for lrst_oce
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   dyn_spg_flt  ! routine called by step.F90
53   PUBLIC   flt_rst      ! routine called by istate.F90
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
60   !! $Id$
61   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63
64CONTAINS
65
66   SUBROUTINE dyn_spg_flt( kt, kindic )
67      !!----------------------------------------------------------------------
68      !!                  ***  routine dyn_spg_flt  ***
69      !!
70      !! ** Purpose :   Compute the now trend due to the surface pressure
71      !!      gradient in case of filtered free surface formulation  and add
72      !!      it to the general trend of momentum equation.
73      !!
74      !! ** Method  :   Filtered free surface formulation. The surface
75      !!      pressure gradient is given by:
76      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
77      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
78      !!      where sshn is the free surface elevation and btda is the after
79      !!      time derivative of the free surface elevation
80      !!       -1- evaluate the surface presure trend (including the addi-
81      !!      tional force) in three steps:
82      !!        a- compute the right hand side of the elliptic equation:
83      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
84      !!         where (spgu,spgv) are given by:
85      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
86      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
87      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
88      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
89      !!         and define the first guess from previous computation :
90      !!            zbtd = btda
91      !!            btda = 2 zbtd - btdb
92      !!            btdb = zbtd
93      !!        b- compute the relative accuracy to be reached by the
94      !!         iterative solver
95      !!        c- apply the solver by a call to sol... routine
96      !!       -2- compute and add the free surface pressure gradient inclu-
97      !!      ding the additional force used to stabilize the equation.
98      !!
99      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
100      !!
101      !! References : Roullet and Madec 1999, JGR.
102      !!---------------------------------------------------------------------
103      USE oce, ONLY :   zub   => ta   ! ta used as workspace
104      USE oce, ONLY :   zvb   => sa   ! ta used as workspace
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, zraur   ! temporary scalars
111      REAL(wp) ::   zgcb, zbtd   !   -          -
112      REAL(wp) ::   ztdgu, ztdgv         !   -          -
113      !!----------------------------------------------------------------------
114      !
115      IF( kt == nit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
119       
120         ! set to zero free surface specific arrays
121         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
122         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
123         CALL solver_init( nit000 )           ! Elliptic solver initialisation
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      zraur  = 1. / rauw
137
138      !! Explicit physics with thickness weighted updates
139      IF( lk_vvl ) THEN          ! variable volume
140
141         ! Evaluate the masked next velocity (effect of the additional force not included)
142         ! -------------------   (thickness weighted velocity, surface pressure gradient already included in dyn_hpg)
143         DO jk = 1, jpkm1
144            DO jj = 2, jpjm1
145               DO ji = fs_2, fs_jpim1   ! vector opt.
146                  ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
147                     &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
148                     &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
149                  va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
150                     &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
151                     &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
152               END DO
153            END DO
154         END DO
155
156      ELSE                       ! fixed volume
157
158         ! Surface pressure gradient (now)
159         DO jj = 2, jpjm1
160            DO ji = fs_2, fs_jpim1   ! vector opt.
161               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
162               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
163            END DO
164         END DO 
165         !
166         ! add the surface pressure trend to the general trend and
167         ! evaluate the masked next velocity (effect of the additional force not included)
168         DO jk = 1, jpkm1
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.
171                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
172                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
173               END DO
174            END DO
175         END DO
176         !
177      ENDIF
178
179#if defined key_obc
180      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
181      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
182#endif
183#if defined key_bdy
184      ! Update velocities on unstructured boundary using the Flow Relaxation Scheme
185      CALL bdy_dyn_frs( kt )
186
187      IF (ln_bdy_vol) THEN
188      ! Correction of the barotropic component velocity to control the volume of the system
189        CALL bdy_vol( kt )
190      ENDIF
191#endif
192#if defined key_agrif
193      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
194#endif
195#if defined key_orca_r2
196      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
197#endif
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.e0
204            spgv(ji,jj) = 0.e0
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(ji,1,jk) * ua(ji,1,jk)
214               spgv(ji,1) = spgv(ji,1) + fse3v(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(ji,jj,jk) * ua(ji,jj,jk)
222                  spgv(ji,jj) = spgv(ji,jj) + fse3v(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. )   
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 = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
276      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
277
278      epsr = eps * eps * rnorme
279      ncut = 0
280      ! if rnorme is 0, the solution is 0, the solver is not called
281      IF( rnorme == 0.e0 ) THEN
282         gcx(:,:) = 0.e0
283         res   = 0.e0
284         niter = 0
285         ncut  = 999
286      ENDIF
287
288      ! Evaluate the next transport divergence
289      ! --------------------------------------
290      !    Iterarive solver for the elliptic equation (except IF sol.=0)
291      !    (output in gcx with boundary conditions applied)
292      kindic = 0
293      IF( ncut == 0 ) THEN
294         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
295         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
296         ENDIF
297      ENDIF
298
299      ! Transport divergence gradient multiplied by z2dt
300      ! --------------------------------------------====
301      DO jj = 2, jpjm1
302         DO ji = fs_2, fs_jpim1   ! vector opt.
303            ! trend of Transport divergence gradient
304            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
305            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
306            ! multiplied by z2dt
307#if defined key_obc
308            ! caution : grad D = 0 along open boundaries
309            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
310            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
311#elif defined key_bdy
312            ! caution : grad D = 0 along open boundaries
313            ! Remark: The filtering force could be reduced here in the FRS zone
314            !         by multiplying spgu/spgv by (1-alpha) ?? 
315            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
316            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)           
317#else
318            spgu(ji,jj) = z2dt * ztdgu
319            spgv(ji,jj) = z2dt * ztdgv
320#endif
321         END DO
322      END DO
323
324#if defined key_agrif     
325      IF( .NOT. Agrif_Root() ) THEN
326         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
327         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
328         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
329         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
330         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
331      ENDIF
332#endif     
333      ! Add the trends multiplied by z2dt to the after velocity
334      ! -------------------------------------------------------
335      !     ( c a u t i o n : (ua,va) here are the after velocity not the
336      !                       trend, the leap-frog time stepping will not
337      !                       be done in dynnxt.F90 routine)
338      DO jk = 1, jpkm1
339         DO jj = 2, jpjm1
340            DO ji = fs_2, fs_jpim1   ! vector opt.
341               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
342               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
343            END DO
344         END DO
345      END DO
346
347      ! write filtered free surface arrays in restart file
348      ! --------------------------------------------------
349      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
350      !
351   END SUBROUTINE dyn_spg_flt
352
353
354   SUBROUTINE flt_rst( kt, cdrw )
355     !!---------------------------------------------------------------------
356     !!                   ***  ROUTINE ts_rst  ***
357     !!
358     !! ** Purpose : Read or write filtered free surface arrays in restart file
359     !!----------------------------------------------------------------------
360     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
361     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
362     !!----------------------------------------------------------------------
363
364     IF( TRIM(cdrw) == 'READ' ) THEN
365        IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
366! Caution : extra-hallow
367! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
368           CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
369           CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
370           IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
371        ELSE
372           gcx (:,:) = 0.e0
373           gcxb(:,:) = 0.e0
374        ENDIF
375     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
376! Caution : extra-hallow
377! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
378        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
379        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
380     ENDIF
381     !
382   END SUBROUTINE flt_rst
383
384#else
385   !!----------------------------------------------------------------------
386   !!   Default case :   Empty module   No standart free surface cst volume
387   !!----------------------------------------------------------------------
388CONTAINS
389   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
390      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
391   END SUBROUTINE dyn_spg_flt
392   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
393      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
394      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
395      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
396   END SUBROUTINE flt_rst
397#endif
398   
399   !!======================================================================
400END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.