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

Last change on this file since 719 was 719, checked in by ctlod, 14 years ago

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 22.6 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: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DYN/dynspg_flt.F90,v 1.14 2007/06/05 10:38:27 opalod Exp $
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         CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua )
190         CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va )
191
192
193         ! Thickness weighting
194         ! -------------------
195         DO jk = 1, jpkm1
196            DO jj = 1, jpj
197               DO ji = 1, jpi
198                  ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
199                  va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
200
201                  zub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
202                  zvb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
203               END DO
204            END DO
205         END DO
206
207         ! Evaluate the masked next velocity (effect of the additional force not included)
208         DO jk = 1, jpkm1
209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1   ! vector opt.
211                  ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk)
212                  va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk)
213               END DO
214            END DO
215         END DO
216
217      ELSE                       ! fixed volume
218
219         ! Surface pressure gradient (now)
220         DO jj = 2, jpjm1
221            DO ji = fs_2, fs_jpim1   ! vector opt.
222               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
223               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
224            END DO
225         END DO 
226
227         ! Add the surface pressure trend to the general trend
228         DO jk = 1, jpkm1
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
232                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
233               END DO
234            END DO
235         END DO
236
237         ! Evaluate the masked next velocity (effect of the additional force not included)
238         DO jk = 1, jpkm1
239            DO jj = 2, jpjm1
240               DO ji = fs_2, fs_jpim1   ! vector opt.
241                  ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
242                  va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
243               END DO
244            END DO
245         END DO
246
247      ENDIF
248
249#if defined key_obc
250      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
251      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
252#endif
253#if defined key_agrif
254      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
255#endif
256#if defined key_orca_r2
257      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
258#endif
259
260      ! compute the next vertically averaged velocity (effect of the additional force not included)
261      ! ---------------------------------------------
262      DO jj = 2, jpjm1
263         DO ji = fs_2, fs_jpim1   ! vector opt.
264            spgu(ji,jj) = 0.e0
265            spgv(ji,jj) = 0.e0
266         END DO
267      END DO
268
269      ! vertical sum
270!CDIR NOLOOPCHG
271      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
272         DO jk = 1, jpkm1
273            DO ji = 1, jpij
274               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
275               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
276            END DO
277         END DO
278      ELSE                        ! No  vector opt.
279         DO jk = 1, jpkm1
280            DO jj = 2, jpjm1
281               DO ji = 2, jpim1
282                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
283                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
284               END DO
285            END DO
286         END DO
287      ENDIF
288
289      ! transport: multiplied by the horizontal scale factor
290      DO jj = 2, jpjm1
291         DO ji = fs_2, fs_jpim1   ! vector opt.
292            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
293            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
294         END DO
295      END DO
296
297      ! Boundary conditions on (spgu,spgv)
298      CALL lbc_lnk( spgu, 'U', -1. )
299      CALL lbc_lnk( spgv, 'V', -1. )
300
301      IF( lk_vvl ) CALL sol_mat( kt )
302
303      ! Right hand side of the elliptic equation and first guess
304      ! -----------------------------------------------------------
305      DO jj = 2, jpjm1
306         DO ji = fs_2, fs_jpim1   ! vector opt.
307            ! Divergence of the after vertically averaged velocity
308            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
309                  + spgv(ji,jj) - spgv(ji,jj-1)
310            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
311            ! First guess of the after barotropic transport divergence
312            zbtd = gcx(ji,jj)
313            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
314            gcxb(ji,jj) =      zbtd
315         END DO
316      END DO
317      ! applied the lateral boundary conditions
318      IF( nsolv == 4 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
319
320#if defined key_agrif
321      IF( .NOT. AGRIF_ROOT() ) THEN
322         ! add contribution of gradient of after barotropic transport divergence
323         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
324            &    gcb(3     ,:) - znugdt * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
325         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
326            &    gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
327         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
328            &    gcb(:,3     ) - znugdt * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
329         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
330            &    gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
331      ENDIF
332#endif
333
334
335      ! Relative precision (computation on one processor)
336      ! ------------------
337      rnorme =0.
338      rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
339      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
340
341      epsr = eps * eps * rnorme
342      ncut = 0
343      ! if rnorme is 0, the solution is 0, the solver is not called
344      IF( rnorme == 0.e0 ) THEN
345         gcx(:,:) = 0.e0
346         res   = 0.e0
347         niter = 0
348         ncut  = 999
349      ENDIF
350
351      ! Evaluate the next transport divergence
352      ! --------------------------------------
353      !    Iterarive solver for the elliptic equation (except IF sol.=0)
354      !    (output in gcx with boundary conditions applied)
355      kindic = 0
356      IF( ncut == 0 ) THEN
357         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
358            CALL sol_pcg( kindic )
359         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
360            CALL sol_sor( kindic )
361         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
362            CALL sol_fet( kindic )
363         ELSEIF( nsolv == 4 ) THEN     ! successive-over-relaxation with extra outer halo
364            CALL sol_sor_e( kindic )
365         ELSE                          ! e r r o r in nsolv namelist parameter
366            WRITE(ctmp1,*) ' ~~~~~~~~~~~                not = ', nsolv
367            CALL ctl_stop( ' dyn_spg_flt : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
368         ENDIF
369      ENDIF
370
371      ! Transport divergence gradient multiplied by z2dt
372      ! --------------------------------------------====
373      DO jj = 2, jpjm1
374         DO ji = fs_2, fs_jpim1   ! vector opt.
375            ! trend of Transport divergence gradient
376            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
377            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
378            ! multiplied by z2dt
379#if defined key_obc
380            ! caution : grad D = 0 along open boundaries
381            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
382            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
383#else
384            spgu(ji,jj) = z2dt * ztdgu
385            spgv(ji,jj) = z2dt * ztdgv
386#endif
387         END DO
388      END DO
389
390#if defined key_agrif     
391      IF( .NOT. Agrif_Root() ) THEN
392         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
393         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = znugdt * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
394         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
395         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = znugdt * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
396         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
397      ENDIF
398#endif     
399      ! 7.  Add the trends multiplied by z2dt to the after velocity
400      ! -----------------------------------------------------------
401      !     ( c a u t i o n : (ua,va) here are the after velocity not the
402      !                       trend, the leap-frog time stepping will not
403      !                       be done in dynnxt.F90 routine)
404      DO jk = 1, jpkm1
405         DO jj = 2, jpjm1
406            DO ji = fs_2, fs_jpim1   ! vector opt.
407               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
408               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
409            END DO
410         END DO
411      END DO
412
413      ! Sea surface elevation time stepping
414      ! -----------------------------------
415      IF( lk_vvl ) THEN   ! after free surface elevation
416         zssha(:,:) = ssha(:,:)
417      ELSE
418         zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1)
419      ENDIF
420
421      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter
422         ! swap of arrays
423         sshb(:,:) = sshn (:,:)
424         sshn(:,:) = zssha(:,:)
425      ELSE                                           ! Leap-frog time stepping and time filter
426         ! time filter and array swap
427         sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:)
428         sshn(:,:) = zssha(:,:)
429      ENDIF
430
431      ! write filtered free surface arrays in restart file
432      ! --------------------------------------------------
433      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
434
435      ! print sum trends (used for debugging)
436      IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
437      !
438   END SUBROUTINE dyn_spg_flt
439
440
441   SUBROUTINE flt_rst( kt, cdrw )
442     !!---------------------------------------------------------------------
443     !!                   ***  ROUTINE ts_rst  ***
444     !!
445     !! ** Purpose : Read or write filtered free surface arrays in restart file
446     !!----------------------------------------------------------------------
447     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
448     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
449     !!----------------------------------------------------------------------
450
451     IF( TRIM(cdrw) == 'READ' ) THEN
452        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
453! Caution : extra-hallow
454! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
455           CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
456           CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
457           CALL iom_get( numror, jpdom_autoglo, 'sshb', sshb(:,:)         )
458           CALL iom_get( numror, jpdom_autoglo, 'sshn', sshn(:,:)         )
459           IF( neuler == 0 ) THEN
460              sshb(:,:) = sshn(:,:)
461              gcxb(:,:) = gcx (:,:)
462           ENDIF
463        ELSE
464           gcx (:,:) = 0.e0
465           gcxb(:,:) = 0.e0
466           IF( nn_rstssh == 1 ) THEN 
467              sshb(:,:) = 0.e0
468              sshn(:,:) = 0.e0
469           ENDIF
470        ENDIF
471     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
472! Caution : extra-hallow
473! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
474        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
475        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
476        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
477        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
478     ENDIF
479     !
480   END SUBROUTINE flt_rst
481
482#else
483   !!----------------------------------------------------------------------
484   !!   Default case :   Empty module   No standart free surface cst volume
485   !!----------------------------------------------------------------------
486CONTAINS
487   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
488      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
489   END SUBROUTINE dyn_spg_flt
490   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
491      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
492      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
493      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
494   END SUBROUTINE flt_rst
495#endif
496   
497   !!======================================================================
498END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.