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_exp.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90 @ 579

Last change on this file since 579 was 575, checked in by opalod, 18 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: 10.6 KB
Line 
1MODULE dynspg_exp
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_exp  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if defined key_dynspg_exp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_exp'           free sfce cst vol. without filter nor ts
9   !!----------------------------------------------------------------------
10   !!   dyn_spg_exp  : update the momentum trend with the surface
11   !!                      pressure gradient in the free surface constant 
12   !!                      volume case with vector optimization
13   !!   exp_rst      : read/write the explicit restart fields in the ocean restart file
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE phycst          ! physical constants
20   USE ocesbc          ! ocean surface boundary condition
21   USE obc_oce         ! Lateral open boundary condition
22   USE obc_par         ! open boundary condition parameters
23   USE obcdta          ! open boundary condition data     (obc_dta_bt routine)
24   USE lib_mpp         ! distributed memory computing library
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE prtctl          ! Print control
27   USE iom             ! I/O library
28   USE restart         ! only for lrst_oce
29
30   IMPLICIT NONE
31   PRIVATE
32
33   !! * Accessibility
34   PUBLIC dyn_spg_exp  ! routine called by step.F90
35   PUBLIC exp_rst      ! routine called j-k-i subroutine
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !!  OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Header$
43   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
44   !!----------------------------------------------------------------------
45
46
47CONTAINS
48
49   SUBROUTINE dyn_spg_exp( kt )
50      !!----------------------------------------------------------------------
51      !!                  ***  routine dyn_spg_exp  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the surface pressure
54      !!      gradient in case of explicit free surface formulation and
55      !!      add it to the general trend of momentum equation. Compute
56      !!      the free surface.
57      !!
58      !! ** Method  :   Explicit free surface formulation. The surface pressure
59      !!      gradient is given by:
60      !!         spgu = 1/rau0 d/dx(ps) =  g/e1u di( sshn )
61      !!         spgv = 1/rau0 d/dy(ps) =  g/e2v dj( sshn )
62      !!      -1- Compute the now surface pressure gradient
63      !!      -2- Add it to the general trend
64      !!      -3- Compute the horizontal divergence of velocities
65      !!      - the now divergence is given by :
66      !!         zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
67      !!      - integrate the horizontal divergence from the bottom
68      !!         to the surface
69      !!      - apply lateral boundary conditions on zhdivn
70      !!      -4- Estimate the after sea surface elevation from the kinematic
71      !!         surface boundary condition:
72      !!         zssha = sshb - 2 rdt ( zhdiv + emp )
73      !!      - Time filter applied on now sea surface elevation to avoid
74      !!         the divergence of two consecutive time-steps and swap of free
75      !!         surface arrays to start the next time step:
76      !!         sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
77      !!         sshn = zssha
78      !!      - apply lateral boundary conditions on sshn
79      !!
80      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
81      !!
82      !! References :
83      !!
84      !! History :
85      !!   9.0  !  05-11  (V. Garnier, G. Madec, L. Bessieres) Original code
86      !!
87      !!---------------------------------------------------------------------
88      !! * Arguments
89      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
90
91      !! * Local declarations
92      INTEGER  ::   ji, jj, jk               ! dummy loop indices
93      REAL(wp) ::   z2dt, zraur, zssha       ! temporary scalars
94      REAL(wp), DIMENSION(jpi,jpj)    ::  &  ! temporary arrays
95         &         zhdiv
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'dyn_spg_exp : surface pressure gradient trend'
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (explicit free surface)'
102
103         ! set to zero free surface specific arrays
104         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
105         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
106
107         CALL exp_rst( nit000, 'READ' )       ! read or initialize the following fields:
108         !                                    ! sshb, sshn
109
110      ENDIF
111
112      ! 0. Initialization
113      ! -----------------
114      ! read or estimate sea surface height and vertically integrated velocities
115      IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
116      z2dt = 2. * rdt                                       ! time step: leap-frog
117      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
118      zraur = 1.e0 / rauw
119
120
121      ! 1. Surface pressure gradient (now)
122      ! ----------------------------
123      DO jj = 2, jpjm1
124         DO ji = fs_2, fs_jpim1   ! vector opt.
125            spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
126            spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
127         END DO
128      END DO
129
130      ! 2. Add the surface pressure trend to the general trend
131      ! ------------------------------------------------------
132      DO jk = 1, jpkm1
133         DO jj = 2, jpjm1
134            DO ji = fs_2, fs_jpim1   ! vector opt.
135               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
136               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
137            END DO
138         END DO
139      END DO
140     
141      ! 3. Vertical integration of horizontal divergence of velocities
142      ! --------------------------------
143      zhdiv(:,:) = 0.e0
144      DO jk = jpkm1, 1, -1
145         DO jj = 2, jpjm1
146            DO ji = fs_2, fs_jpim1   ! vector opt.
147               zhdiv(ji,jj) = zhdiv(ji,jj) + (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)      &
148                  &                           - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)      &
149                  &                           + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)      &
150                  &                           - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )   &
151                  &                        / ( e1t(ji,jj) * e2t(ji,jj) )
152            END DO
153         END DO
154      END DO
155
156#if defined key_obc
157      ! open boundaries (div must be zero behind the open boundary)
158      !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column
159      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
160      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
161      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
162      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
163#endif
164
165
166      ! 4. Sea surface elevation time stepping
167      ! --------------------------------------
168      ! Euler (forward) time stepping, no time filter
169      IF( neuler == 0 .AND. kt == nit000 ) THEN
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172               ! after free surface elevation
173               zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)
174               ! swap of arrays
175               sshb(ji,jj) = sshn(ji,jj)
176               sshn(ji,jj) = zssha
177            END DO
178         END DO
179      ELSE
180         ! Leap-frog time stepping and time filter
181         DO jj = 1, jpj
182            DO ji = 1, jpi
183               ! after free surface elevation
184               zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)
185               ! time filter and array swap
186               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
187               sshn(ji,jj) = zssha
188            END DO
189         END DO
190      ENDIF
191
192      ! Boundary conditions on sshn
193      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
194
195      ! write filtered free surface arrays in restart file
196      ! --------------------------------------------------
197      IF( lrst_oce )   CALL exp_rst( kt, 'WRITE' )
198 
199      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
200         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
201      ENDIF
202     
203   END SUBROUTINE dyn_spg_exp
204
205   SUBROUTINE exp_rst( kt, cdrw )
206      !!---------------------------------------------------------------------
207      !!                   ***  ROUTINE exp_rst  ***
208      !!
209      !! ** Purpose : Read or write explicit arrays in restart file
210      !!----------------------------------------------------------------------
211      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
212      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
213      !
214      !!----------------------------------------------------------------------
215      !
216      IF( TRIM(cdrw) == 'READ' ) THEN
217         IF( iom_varid( numror, 'sshn' ) > 0 ) THEN
218            CALL iom_get( numror, jpdom_local, 'sshb'  , sshb(:,:)   )
219            CALL iom_get( numror, jpdom_local, 'sshn'  , sshn(:,:)   )
220            IF( neuler == 0 ) sshb(:,:) = sshn(:,:)
221         ELSE
222            IF( nn_rstssh == 1 ) THEN 
223               sshb(:,:) = 0.e0
224               sshn(:,:) = 0.e0
225            ENDIF
226         ENDIF
227      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
228         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
229         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
230      ENDIF
231      !
232   END SUBROUTINE exp_rst
233#else
234   !!----------------------------------------------------------------------
235   !!   Default case :   Empty module   No standart explicit free surface
236   !!----------------------------------------------------------------------
237CONTAINS
238   SUBROUTINE dyn_spg_exp( kt )       ! Empty routine
239      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
240   END SUBROUTINE dyn_spg_exp
241#endif
242   
243   !!======================================================================
244END MODULE dynspg_exp
Note: See TracBrowser for help on using the repository browser.