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

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 11.9 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 sbc_oce         ! surface boundary condition: ocean
19   USE obc_oce         ! Lateral open boundary condition
20   USE phycst          ! physical constants
21   USE obc_par         ! open boundary condition parameters
22   USE obcdta          ! open boundary condition data     (obc_dta_bt routine)
23   USE in_out_manager  ! I/O manager
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 by istate.F90
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !!  OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Id$
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      IF( .NOT. lk_vvl ) THEN
113
114         ! 0. Initialization
115         ! -----------------
116         ! read or estimate sea surface height and vertically integrated velocities
117         IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
118         z2dt = 2. * rdt                                       ! time step: leap-frog
119         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
120         zraur = 1.e0 / rauw
121
122         ! 1. Surface pressure gradient (now)
123         ! ----------------------------
124         DO jj = 2, jpjm1
125            DO ji = fs_2, fs_jpim1   ! vector opt.
126               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
127               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
128            END DO
129         END DO
130
131         ! 2. Add the surface pressure trend to the general trend
132         ! ------------------------------------------------------
133         DO jk = 1, jpkm1
134            DO jj = 2, jpjm1
135               DO ji = fs_2, fs_jpim1   ! vector opt.
136                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
137                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
138               END DO
139            END DO
140         END DO
141     
142         ! 3. Vertical integration of horizontal divergence of velocities
143         ! --------------------------------
144         zhdiv(:,:) = 0.e0
145         DO jk = jpkm1, 1, -1
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  zhdiv(ji,jj) = zhdiv(ji,jj) + (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)      &
149                     &                           - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)      &
150                     &                           + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)      &
151                     &                           - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )   &
152                     &                        / ( e1t(ji,jj) * e2t(ji,jj) )
153               END DO
154            END DO
155         END DO
156
157#if defined key_obc
158         ! open boundaries (div must be zero behind the open boundary)
159         !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column
160         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
161         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
162         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
163         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
164#endif
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      ELSE !! Variable volume, ssh time-stepping already done
193
194         ! read or estimate sea surface height and vertically integrated velocities
195         IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
196
197         ! Sea surface elevation swap
198         ! -----------------------------
199         !
200         sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping
201         !
202         IF( neuler == 0 .AND. kt == nit000 ) THEN
203            ! No time filter
204            sshb(:,:) = sshn(:,:)
205            sshn(:,:) = ssha(:,:)
206         ELSE
207            ! Time filter
208            sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)
209            sshn(:,:) = ssha(:,:)
210         ENDIF
211
212      ENDIF
213
214      ! Boundary conditions on sshn
215      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
216
217      ! write filtered free surface arrays in restart file
218      ! --------------------------------------------------
219      IF( lrst_oce )   CALL exp_rst( kt, 'WRITE' )
220 
221      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
222         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
223      ENDIF
224     
225   END SUBROUTINE dyn_spg_exp
226
227   SUBROUTINE exp_rst( kt, cdrw )
228      !!---------------------------------------------------------------------
229      !!                   ***  ROUTINE exp_rst  ***
230      !!
231      !! ** Purpose : Read or write explicit arrays in restart file
232      !!----------------------------------------------------------------------
233      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
234      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
235      !
236      !!----------------------------------------------------------------------
237      !
238      IF( TRIM(cdrw) == 'READ' ) THEN
239         IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
240            CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   )
241            CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   )
242            IF( neuler == 0 ) sshb(:,:) = sshn(:,:)
243         ELSE
244            IF( nn_rstssh == 1 ) THEN 
245               sshb(:,:) = 0.e0
246               sshn(:,:) = 0.e0
247            ENDIF
248         ENDIF
249      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
250         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
251         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
252      ENDIF
253      !
254   END SUBROUTINE exp_rst
255#else
256   !!----------------------------------------------------------------------
257   !!   Default case :   Empty module   No standart explicit free surface
258   !!----------------------------------------------------------------------
259CONTAINS
260   SUBROUTINE dyn_spg_exp( kt )       ! Empty routine
261      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
262   END SUBROUTINE dyn_spg_exp
263   SUBROUTINE exp_rst( kt, cdrw )     ! Empty routine
264      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
265      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
266      WRITE(*,*) 'exp_rst: You should not have seen this print! error?', kt, cdrw
267   END SUBROUTINE exp_rst
268#endif
269   
270   !!======================================================================
271END MODULE dynspg_exp
Note: See TracBrowser for help on using the repository browser.