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 @ 592

Last change on this file since 592 was 592, checked in by opalod, 17 years ago

nemo_v2_update_001 : CT : - add non linear free surface (variable volume) with new cpp key key_vvl

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 22.2 KB
Line 
1MODULE dynspg_flt
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_flt  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6   !! History    8.0  !  98-05  (G. Roullet)  free surface
7   !!                 !  98-10  (G. Madec, M. Imbard)  release 8.2
8   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
9   !!            " "  !  02-11  (C. Talandier, A-M Treguier) Open boundaries
10   !!            9.0  !  04-08  (C. Talandier) New trends organization
11   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
12   !!            " "  !  06-07  (S. Masson)  distributed restart using iom
13   !!----------------------------------------------------------------------
14#if defined key_dynspg_flt   ||   defined key_esopa 
15   !!----------------------------------------------------------------------
16   !!   'key_dynspg_flt'                              filtered free surface
17   !!----------------------------------------------------------------------
18   !!----------------------------------------------------------------------
19   !!   dyn_spg_flt  : update the momentum trend with the surface pressure
20   !!                  gradient in the filtered free surface case with
21   !!                  vector optimization
22   !!   flt_rst      : read/write the time-splitting restart fields in the ocean restart file
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE zdf_oce         ! ocean vertical physics
27   USE phycst          ! physical constants
28   USE ocesbc          ! ocean surface boundary condition
29   USE flxrnf          ! ocean runoffs
30   USE sol_oce         ! ocean elliptic solver
31   USE solver          ! solver initialization
32   USE solpcg          ! preconditionned conjugate gradient solver
33   USE solsor          ! Successive Over-relaxation solver
34   USE solfet          ! FETI solver
35   USE solsor_e        ! Successive Over-relaxation solver with MPP optimization
36   USE obc_oce         ! Lateral open boundary condition
37   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
38   USE obcvol          ! ocean open boundary condition (obc_vol routines)
39   USE lib_mpp         ! distributed memory computing library
40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
41   USE cla_dynspg      ! cross land advection
42   USE prtctl          ! Print control
43   USE solmat          ! matrix construction for elliptic solvers
44   USE agrif_opa_interp
45   USE in_out_manager  ! I/O manager
46   USE iom
47   USE restart         ! only for lrst_oce
48   USE domvvl          ! variable volume
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC dyn_spg_flt  ! routine called by step.F90
54   PUBLIC flt_rst      ! routine called by j-k-i subroutine
55
56   !! * Substitutions
57#  include "domzgr_substitute.h90"
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !!   OPA 9.0 , LOCEAN-IPSL (2005)
61   !! $Header$
62   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64
65CONTAINS
66
67   SUBROUTINE dyn_spg_flt( kt, kindic )
68      !!----------------------------------------------------------------------
69      !!                  ***  routine dyn_spg_flt  ***
70      !!
71      !! ** Purpose :   Compute the now trend due to the surface pressure
72      !!      gradient in case of filtered free surface formulation  and add
73      !!      it to the general trend of momentum equation.
74      !!
75      !! ** Method  :   Filtered free surface formulation. The surface
76      !!      pressure gradient is given by:
77      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + rnu btda )
78      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + rnu btda )
79      !!      where sshn is the free surface elevation and btda is the after
80      !!      of the free surface elevation
81      !!       -1- compute the after sea surface elevation from the kinematic
82      !!      surface boundary condition:
83      !!              zssha = sshb + 2 rdt ( wn - emp )
84      !!           Time filter applied on now sea surface elevation to avoid
85      !!      the divergence of two consecutive time-steps and swap of free
86      !!      surface arrays to start the next time step:
87      !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
88      !!              sshn = zssha
89      !!       -2- evaluate the surface presure trend (including the addi-
90      !!      tional force) in three steps:
91      !!        a- compute the right hand side of the elliptic equation:
92      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
93      !!         where (spgu,spgv) are given by:
94      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
95      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
96      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
97      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
98      !!         and define the first guess from previous computation :
99      !!            zbtd = btda
100      !!            btda = 2 zbtd - btdb
101      !!            btdb = zbtd
102      !!        b- compute the relative accuracy to be reached by the
103      !!         iterative solver
104      !!        c- apply the solver by a call to sol... routine
105      !!       -3- compute and add the free surface pressure gradient inclu-
106      !!      ding the additional force used to stabilize the equation.
107      !!
108      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
109      !!
110      !! References : Roullet and Madec 1999, JGR.
111      !!---------------------------------------------------------------------
112      !! * Modules used
113      USE oce    , ONLY :   zub   => ta,             &  ! ta used as workspace
114                            zvb   => sa                 ! sa   "          "
115
116      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
117      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag (<0 if not converge)
118      !!                                   
119      INTEGER  ::   ji, jj, jk                          ! dummy loop indices
120      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
121         &          znurau, zgcb, zbtd,              &  !   "          "
122         &          ztdgu, ztdgv                        !   "          "
123      !! Variable volume
124      REAL(wp), DIMENSION(jpi,jpj)     ::  &
125         zsshub, zsshua, zsshvb, zsshva, zssha          ! 2D workspace
126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &            ! 3D workspace
127         zfse3ub, zfse3ua, zfse3vb, zfse3va             ! 3D workspace
128      !!----------------------------------------------------------------------
129      !
130      IF( kt == nit000 ) THEN
131         IF(lwp) WRITE(numout,*)
132         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
133         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
134       
135         ! set to zero free surface specific arrays
136         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
137         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
138         CALL solver_init( nit000 )           ! Elliptic solver initialisation
139
140         ! read filtered free surface arrays in restart file
141         CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
142         !                                    ! gcx, gcxb, sshb, sshn
143      ENDIF
144
145      ! Local constant initialization
146      z2dt = 2. * rdt                                       ! time step: leap-frog
147      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
148      IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat(kt)
149      z2dtg  = grav * z2dt
150      zraur  = 1. / rauw
151      znugdt =  rnu * grav * z2dt
152      znurau =  znugdt * zraur
153
154      !! Explicit physics with thickness weighted updates
155      IF( lk_vvl ) THEN          ! variable volume
156
157         DO jj = 1, jpjm1
158            DO ji = 1,jpim1
159
160               ! Sea Surface Height at u-point before
161               zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  & 
162                  &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)       &
163                  &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
164
165               ! Sea Surface Height at v-point before
166               zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  &
167                  &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )       &
168                  &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
169
170               ! Sea Surface Height at u-point after
171               zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  &
172                  &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)       &
173                  &            + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
174
175               ! Sea Surface Height at v-point after
176               zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  &
177                  &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )       &
178                  &            + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
179
180            END DO
181         END DO
182
183         ! Boundaries conditions
184         CALL lbc_lnk( zsshub, 'U', 1. )    ;    CALL lbc_lnk( zsshvb, 'V', 1. )
185         CALL lbc_lnk( zsshua, 'U', 1. )    ;    CALL lbc_lnk( zsshva, 'V', 1. )
186
187         ! Scale factors at before and after time step
188         ! -------------------------------------------
189         DO jk = 1, jpkm1
190            zfse3ua(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zsshua(:,:) * muu(:,:,jk) )
191            zfse3va(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zsshva(:,:) * muv(:,:,jk) )
192            zfse3ub(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zsshub(:,:) * muu(:,:,jk) )
193            zfse3vb(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zsshvb(:,:) * muv(:,:,jk) )
194         END DO
195
196         ! Thickness weighting
197         ! -------------------
198         ua(:,:,1:jpkm1) = ua(:,:,1:jpkm1) * fse3u(:,:,1:jpkm1)
199         va(:,:,1:jpkm1) = va(:,:,1:jpkm1) * fse3v(:,:,1:jpkm1)
200
201         zub(:,:,1:jpkm1) = ub(:,:,1:jpkm1) * zfse3ub(:,:,1:jpkm1)
202         zvb(:,:,1:jpkm1) = vb(:,:,1:jpkm1) * zfse3vb(:,:,1:jpkm1)
203
204         ! Evaluate the masked next velocity (effect of the additional force not included)
205         DO jk = 1, jpkm1
206            DO jj = 2, jpjm1
207               DO ji = fs_2, fs_jpim1   ! vector opt.
208                  ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk)
209                  va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk)
210               END DO
211            END DO
212         END DO
213
214      ELSE                       ! fixed volume
215
216         ! Surface pressure gradient (now)
217         DO jj = 2, jpjm1
218            DO ji = fs_2, fs_jpim1   ! vector opt.
219               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
220               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
221            END DO
222         END DO 
223
224         ! Add the surface pressure trend to the general trend
225         DO jk = 1, jpkm1
226            DO jj = 2, jpjm1
227               DO ji = fs_2, fs_jpim1   ! vector opt.
228                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
229                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
230               END DO
231            END DO
232         END DO
233
234         ! Evaluate the masked next velocity (effect of the additional force not included)
235         DO jk = 1, jpkm1
236            DO jj = 2, jpjm1
237               DO ji = fs_2, fs_jpim1   ! vector opt.
238                  ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
239                  va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
240               END DO
241            END DO
242         END DO
243
244      ENDIF
245
246#if defined key_obc
247      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
248      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
249#endif
250#if defined key_agrif
251      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
252#endif
253#if defined key_orca_r2
254      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
255#endif
256
257      ! compute the next vertically averaged velocity (effect of the additional force not included)
258      ! ---------------------------------------------
259      DO jj = 2, jpjm1
260         DO ji = fs_2, fs_jpim1   ! vector opt.
261            spgu(ji,jj) = 0.e0
262            spgv(ji,jj) = 0.e0
263         END DO
264      END DO
265
266      ! vertical sum
267!CDIR NOLOOPCHG
268      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
269         DO jk = 1, jpkm1
270            DO ji = 1, jpij
271               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
272               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
273            END DO
274         END DO
275      ELSE                        ! No  vector opt.
276         DO jk = 1, jpkm1
277            DO jj = 2, jpjm1
278               DO ji = 2, jpim1
279                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
280                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
281               END DO
282            END DO
283         END DO
284      ENDIF
285
286      ! transport: multiplied by the horizontal scale factor
287      DO jj = 2, jpjm1
288         DO ji = fs_2, fs_jpim1   ! vector opt.
289            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
290            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
291         END DO
292      END DO
293
294      ! Boundary conditions on (spgu,spgv)
295      CALL lbc_lnk( spgu, 'U', -1. )
296      CALL lbc_lnk( spgv, 'V', -1. )
297
298      IF( lk_vvl ) CALL sol_mat( kt )
299
300      ! Right hand side of the elliptic equation and first guess
301      ! -----------------------------------------------------------
302      DO jj = 2, jpjm1
303         DO ji = fs_2, fs_jpim1   ! vector opt.
304            ! Divergence of the after vertically averaged velocity
305            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
306                  + spgv(ji,jj) - spgv(ji,jj-1)
307            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
308            ! First guess of the after barotropic transport divergence
309            zbtd = gcx(ji,jj)
310            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
311            gcxb(ji,jj) =      zbtd
312         END DO
313      END DO
314      ! applied the lateral boundary conditions
315      IF( nsolv == 4 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
316
317#if defined key_agrif
318      IF( .NOT. AGRIF_ROOT() ) THEN
319         ! add contribution of gradient of after barotropic transport divergence
320         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
321            &    gcb(3     ,:) - znugdt * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
322         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
323            &    gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
324         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
325            &    gcb(:,3     ) - znugdt * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
326         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
327            &    gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
328      ENDIF
329#endif
330
331
332      ! Relative precision (computation on one processor)
333      ! ------------------
334      rnorme =0.
335      rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
336      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
337
338      epsr = eps * eps * rnorme
339      ncut = 0
340      ! if rnorme is 0, the solution is 0, the solver is not called
341      IF( rnorme == 0.e0 ) THEN
342         gcx(:,:) = 0.e0
343         res   = 0.e0
344         niter = 0
345         ncut  = 999
346      ENDIF
347
348      ! Evaluate the next transport divergence
349      ! --------------------------------------
350      !    Iterarive solver for the elliptic equation (except IF sol.=0)
351      !    (output in gcx with boundary conditions applied)
352      kindic = 0
353      IF( ncut == 0 ) THEN
354         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
355            CALL sol_pcg( kindic )
356         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
357            CALL sol_sor( kindic )
358         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
359            CALL sol_fet( kindic )
360         ELSEIF( nsolv == 4 ) THEN     ! successive-over-relaxation with extra outer halo
361            CALL sol_sor_e( kindic )
362         ELSE                          ! e r r o r in nsolv namelist parameter
363            WRITE(ctmp1,*) ' ~~~~~~~~~~~                not = ', nsolv
364            CALL ctl_stop( ' dyn_spg_flt : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
365         ENDIF
366      ENDIF
367
368      ! Transport divergence gradient multiplied by z2dt
369      ! --------------------------------------------====
370      DO jj = 2, jpjm1
371         DO ji = fs_2, fs_jpim1   ! vector opt.
372            ! trend of Transport divergence gradient
373            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
374            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
375            ! multiplied by z2dt
376#if defined key_obc
377            ! caution : grad D = 0 along open boundaries
378            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
379            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
380#else
381            spgu(ji,jj) = z2dt * ztdgu
382            spgv(ji,jj) = z2dt * ztdgv
383#endif
384         END DO
385      END DO
386
387#if defined key_agrif     
388      IF( .NOT. Agrif_Root() ) THEN
389         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
390         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = znugdt * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
391         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
392         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = znugdt * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
393         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
394      ENDIF
395#endif     
396      ! 7.  Add the trends multiplied by z2dt to the after velocity
397      ! -----------------------------------------------------------
398      !     ( c a u t i o n : (ua,va) here are the after velocity not the
399      !                       trend, the leap-frog time stepping will not
400      !                       be done in dynnxt.F90 routine)
401      DO jk = 1, jpkm1
402         DO jj = 2, jpjm1
403            DO ji = fs_2, fs_jpim1   ! vector opt.
404               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
405               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
406            END DO
407         END DO
408      END DO
409
410      ! Sea surface elevation time stepping
411      ! -----------------------------------
412      IF( lk_vvl ) THEN   ! after free surface elevation
413         zssha(:,:) = ssha(:,:)
414      ELSE
415         zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1)
416      ENDIF
417
418      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter
419         ! swap of arrays
420         sshb(:,:) = sshn (:,:)
421         sshn(:,:) = zssha(:,:)
422      ELSE                                           ! Leap-frog time stepping and time filter
423         ! time filter and array swap
424         sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:)
425         sshn(:,:) = zssha(:,:)
426      ENDIF
427
428      ! write filtered free surface arrays in restart file
429      ! --------------------------------------------------
430      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
431
432      ! print sum trends (used for debugging)
433      IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
434      !
435   END SUBROUTINE dyn_spg_flt
436
437
438   SUBROUTINE flt_rst( kt, cdrw )
439     !!---------------------------------------------------------------------
440     !!                   ***  ROUTINE ts_rst  ***
441     !!
442     !! ** Purpose : Read or write filtered free surface arrays in restart file
443     !!----------------------------------------------------------------------
444     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
445     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
446     !!----------------------------------------------------------------------
447
448     IF( TRIM(cdrw) == 'READ' ) THEN
449        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
450! Caution : extra-hallow
451! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
452           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
453           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
454           CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:)         )
455           CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:)         )
456           IF( neuler == 0 ) THEN
457              sshb(:,:) = sshn(:,:)
458              gcxb(:,:) = gcx (:,:)
459           ENDIF
460        ELSE
461           gcx (:,:) = 0.e0
462           gcxb(:,:) = 0.e0
463           IF( nn_rstssh == 1 ) THEN 
464              sshb(:,:) = 0.e0
465              sshn(:,:) = 0.e0
466           ENDIF
467        ENDIF
468     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
469! Caution : extra-hallow
470! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
471        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
472        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
473        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
474        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
475     ENDIF
476     !
477   END SUBROUTINE flt_rst
478
479#else
480   !!----------------------------------------------------------------------
481   !!   Default case :   Empty module   No standart free surface cst volume
482   !!----------------------------------------------------------------------
483CONTAINS
484   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
485      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
486   END SUBROUTINE dyn_spg_flt
487#endif
488   
489   !!======================================================================
490END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.