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

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

nemo_v1_bugfix_079: CT : - add the creation of restart file (saving sshn & sshb) when using the explicit surface pressure option

  • light modifications in the way to read/write restart files in dynspg_* and associated _jki_ modules (key_mpp_omp removed)
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 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
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC dyn_spg_flt  ! routine called by step.F90
53   PUBLIC flt_rst      ! routine called by j-k-i subroutine
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !!   OPA 9.0 , LOCEAN-IPSL (2005)
60   !! $Header$
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 + rnu btda )
77      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + rnu btda )
78      !!      where sshn is the free surface elevation and btda is the after
79      !!      of the free surface elevation
80      !!       -1- compute the after sea surface elevation from the kinematic
81      !!      surface boundary condition:
82      !!              zssha = sshb + 2 rdt ( wn - emp )
83      !!           Time filter applied on now sea surface elevation to avoid
84      !!      the divergence of two consecutive time-steps and swap of free
85      !!      surface arrays to start the next time step:
86      !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
87      !!              sshn = zssha
88      !!       -2- evaluate the surface presure trend (including the addi-
89      !!      tional force) in three steps:
90      !!        a- compute the right hand side of the elliptic equation:
91      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
92      !!         where (spgu,spgv) are given by:
93      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
94      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
95      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
96      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
97      !!         and define the first guess from previous computation :
98      !!            zbtd = btda
99      !!            btda = 2 zbtd - btdb
100      !!            btdb = zbtd
101      !!        b- compute the relative accuracy to be reached by the
102      !!         iterative solver
103      !!        c- apply the solver by a call to sol... routine
104      !!       -3- compute and add the free surface pressure gradient inclu-
105      !!      ding the additional force used to stabilize the equation.
106      !!
107      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
108      !!
109      !! References : Roullet and Madec 1999, JGR.
110      !!---------------------------------------------------------------------
111      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
112      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag (<0 if not converge)
113      !!                                   
114      INTEGER  ::   ji, jj, jk                          ! dummy loop indices
115      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
116         &          znurau, zssha, zgcb, zbtd,       &  !   "          "
117         &          ztdgu, ztdgv                        !   "          "
118      !!----------------------------------------------------------------------
119      !
120      IF( kt == nit000 ) THEN
121         IF(lwp) WRITE(numout,*)
122         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
123         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
124       
125         ! set to zero free surface specific arrays
126         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
127         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
128         CALL solver_init( nit000 )           ! Elliptic solver initialisation
129
130         ! read filtered free surface arrays in restart file
131         CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
132         !                                    ! gcx, gcxb, sshb, sshn
133      ENDIF
134
135      ! Local constant initialization
136      z2dt = 2. * rdt                                       ! time step: leap-frog
137      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
138      IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat(kt)
139      z2dtg  = grav * z2dt
140      zraur  = 1. / rauw
141      znugdt =  rnu * grav * z2dt
142      znurau =  znugdt * zraur
143
144      ! Surface pressure gradient (now)
145      DO jj = 2, jpjm1
146         DO ji = fs_2, fs_jpim1   ! vector opt.
147            spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
148            spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
149         END DO
150      END DO 
151
152      ! Add the surface pressure trend to the general trend
153      DO jk = 1, jpkm1
154         DO jj = 2, jpjm1
155            DO ji = fs_2, fs_jpim1   ! vector opt.
156               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
157               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
158            END DO
159         END DO
160      END DO
161
162      ! Evaluate the masked next velocity (effect of the additional force not included)
163      DO jk = 1, jpkm1
164         DO jj = 2, jpjm1
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
167               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
168            END DO
169         END DO
170      END DO
171
172#if defined key_obc
173      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
174      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
175#endif
176#if defined key_agrif
177      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
178#endif
179#if defined key_orca_r2
180      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
181#endif
182
183      ! compute the next vertically averaged velocity (effect of the additional force not included)
184      ! ---------------------------------------------
185      DO jj = 2, jpjm1
186         DO ji = fs_2, fs_jpim1   ! vector opt.
187            spgu(ji,jj) = 0.e0
188            spgv(ji,jj) = 0.e0
189         END DO
190      END DO
191
192      ! vertical sum
193!CDIR NOLOOPCHG
194      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
195         DO jk = 1, jpkm1
196            DO ji = 1, jpij
197               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
198               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
199            END DO
200         END DO
201      ELSE                        ! No  vector opt.
202         DO jk = 1, jpkm1
203            DO jj = 2, jpjm1
204               DO ji = 2, jpim1
205                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
206                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
207               END DO
208            END DO
209         END DO
210      ENDIF
211
212      ! transport: multiplied by the horizontal scale factor
213      DO jj = 2, jpjm1
214         DO ji = fs_2, fs_jpim1   ! vector opt.
215            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
216            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
217         END DO
218      END DO
219
220      ! Boundary conditions on (spgu,spgv)
221      CALL lbc_lnk( spgu, 'U', -1. )
222      CALL lbc_lnk( spgv, 'V', -1. )
223
224      ! Right hand side of the elliptic equation and first guess
225      ! -----------------------------------------------------------
226      DO jj = 2, jpjm1
227         DO ji = fs_2, fs_jpim1   ! vector opt.
228            ! Divergence of the after vertically averaged velocity
229            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
230                  + spgv(ji,jj) - spgv(ji,jj-1)
231            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
232            ! First guess of the after barotropic transport divergence
233            zbtd = gcx(ji,jj)
234            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
235            gcxb(ji,jj) =      zbtd
236         END DO
237      END DO
238      ! applied the lateral boundary conditions
239      IF( nsolv == 4 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
240
241#if defined key_agrif
242      IF( .NOT. AGRIF_ROOT() ) THEN
243         ! add contribution of gradient of after barotropic transport divergence
244         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
245            &    gcb(3     ,:) - znugdt * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
246         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
247            &    gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
248         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
249            &    gcb(:,3     ) - znugdt * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
250         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
251            &    gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
252      ENDIF
253#endif
254
255
256      ! Relative precision (computation on one processor)
257      ! ------------------
258      rnorme =0.
259      rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
260      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
261
262      epsr = eps * eps * rnorme
263      ncut = 0
264      ! if rnorme is 0, the solution is 0, the solver is not called
265      IF( rnorme == 0.e0 ) THEN
266         gcx(:,:) = 0.e0
267         res   = 0.e0
268         niter = 0
269         ncut  = 999
270      ENDIF
271
272      ! Evaluate the next transport divergence
273      ! --------------------------------------
274      !    Iterarive solver for the elliptic equation (except IF sol.=0)
275      !    (output in gcx with boundary conditions applied)
276      kindic = 0
277      IF( ncut == 0 ) THEN
278         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
279            CALL sol_pcg( kindic )
280         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
281            CALL sol_sor( kindic )
282         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
283            CALL sol_fet( kindic )
284         ELSEIF( nsolv == 4 ) THEN     ! successive-over-relaxation with extra outer halo
285            CALL sol_sor_e( kindic )
286         ELSE                          ! e r r o r in nsolv namelist parameter
287            WRITE(ctmp1,*) ' ~~~~~~~~~~~                not = ', nsolv
288            CALL ctl_stop( ' dyn_spg_flt : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
289         ENDIF
290      ENDIF
291
292      ! Transport divergence gradient multiplied by z2dt
293      ! --------------------------------------------====
294      DO jj = 2, jpjm1
295         DO ji = fs_2, fs_jpim1   ! vector opt.
296            ! trend of Transport divergence gradient
297            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
298            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
299            ! multiplied by z2dt
300#if defined key_obc
301            ! caution : grad D = 0 along open boundaries
302            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
303            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
304#else
305            spgu(ji,jj) = z2dt * ztdgu
306            spgv(ji,jj) = z2dt * ztdgv
307#endif
308         END DO
309      END DO
310
311#if defined key_agrif     
312      IF( .NOT. Agrif_Root() ) THEN
313         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
314         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = znugdt * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
315         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
316         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = znugdt * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
317         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
318      ENDIF
319#endif     
320      ! 7.  Add the trends multiplied by z2dt to the after velocity
321      ! -----------------------------------------------------------
322      !     ( c a u t i o n : (ua,va) here are the after velocity not the
323      !                       trend, the leap-frog time stepping will not
324      !                       be done in dynnxt.F90 routine)
325      DO jk = 1, jpkm1
326         DO jj = 2, jpjm1
327            DO ji = fs_2, fs_jpim1   ! vector opt.
328               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
329               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
330            END DO
331         END DO
332      END DO
333
334      ! Sea surface elevation time stepping
335      ! -----------------------------------
336      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter
337         DO jj = 1, jpj
338            DO ji = 1, jpi
339               ! after free surface elevation
340               zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
341               ! swap of arrays
342               sshb(ji,jj) = sshn(ji,jj)
343               sshn(ji,jj) = zssha
344            END DO
345         END DO
346      ELSE                                           ! Leap-frog time stepping and time filter
347         DO jj = 1, jpj
348            DO ji = 1, jpi
349               ! after free surface elevation
350               zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
351               ! time filter and array swap
352               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
353               sshn(ji,jj) = zssha
354            END DO
355         END DO
356      ENDIF
357
358      ! write filtered free surface arrays in restart file
359      ! --------------------------------------------------
360      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
361
362      ! print sum trends (used for debugging)
363      IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
364      !
365   END SUBROUTINE dyn_spg_flt
366
367
368   SUBROUTINE flt_rst( kt, cdrw )
369     !!---------------------------------------------------------------------
370     !!                   ***  ROUTINE ts_rst  ***
371     !!
372     !! ** Purpose : Read or write filtered free surface arrays in restart file
373     !!----------------------------------------------------------------------
374     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
375     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
376     !!----------------------------------------------------------------------
377
378     IF( TRIM(cdrw) == 'READ' ) THEN
379        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
380! Caution : extra-hallow
381! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
382           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
383           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
384           CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:)         )
385           CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:)         )
386           IF( neuler == 0 ) THEN
387              sshb(:,:) = sshn(:,:)
388              gcxb(:,:) = gcx (:,:)
389           ENDIF
390        ELSE
391           gcx (:,:) = 0.e0
392           gcxb(:,:) = 0.e0
393           IF( nn_rstssh == 1 ) THEN 
394              sshb(:,:) = 0.e0
395              sshn(:,:) = 0.e0
396           ENDIF
397        ENDIF
398     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
399! Caution : extra-hallow
400! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
401        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
402        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
403        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
404        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
405     ENDIF
406     !
407   END SUBROUTINE flt_rst
408
409#else
410   !!----------------------------------------------------------------------
411   !!   Default case :   Empty module   No standart free surface cst volume
412   !!----------------------------------------------------------------------
413CONTAINS
414   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
415      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
416   END SUBROUTINE dyn_spg_flt
417#endif
418   
419   !!======================================================================
420END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.