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

source: trunk/NEMO/OPA_SRC/OBC/obcvol.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1MODULE obcvol
2   !!=================================================================================
3   !!                       ***  MODULE  obcvol  ***
4   !! Ocean dynamic :  Volume constraint when OBC and Free surface are activated
5   !!=================================================================================
6#if defined key_obc && defined key_dynspg_fsc
7   !!---------------------------------------------------------------------------------
8   !!   'key_obc'               and                           open boundary conditions
9   !!   'key_dynspg_fsc'                                  constant volume free surface
10   !!---------------------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE phycst          ! physical constants
15   USE obc_oce         ! ocean open boundary conditions
16   USE lib_mpp         ! for mppsum
17   USE in_out_manager  ! I/O manager
18   USE ocesbc          ! ocean surface boundary conditions
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Accessibility
24   PUBLIC obc_vol        ! routine called by dynspg_fsc.h90
25
26   !! * Substitutions
27#  include "domzgr_substitute.h90"
28#  include "obc_vectopt_loop_substitute.h90"
29   !!---------------------------------------------------------------------------------
30   !!   OPA 9.0 , LODYC-IPSL  (2003)
31   !!---------------------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE obc_vol ( kt )
36      !!------------------------------------------------------------------------------
37      !!                      ***  ROUTINE obcvol  ***
38      !!
39      !! ** Purpose :
40      !!      This routine is called in dynspg_fsc to control
41      !!      the volume of the system. A correction velocity is calculated
42      !!      to correct the total transport through the OBC.
43      !!      The total depth used is constant (H0) to be consistent with the
44      !!      linear free surface coded in OPA 8.2
45      !!
46      !! ** Method : 
47      !!      The correction velocity (zubtpecor here) is defined calculating
48      !!      the total transport through all open boundaries (trans_obc) minus
49      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral
50      !!      surface (obcsurftot) of these OBC.
51      !!
52      !!      zubtpecor = [trans_obc - zCflxemp ]*(1./obcsurftot)
53      !!
54      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
55      !!                       over all the domain in m3/s at each time step.
56      !!
57      !!      zCflxemp < 0 when precipitation dominate
58      !!      zCflxemp > 0 when evaporation dominate
59      !!
60      !!      There are 2 options (user's desiderata):
61      !!
62      !!         1/ The volume changes according to E-P, this is the default
63      !!            option. In this case the cumulate E-P flux are setting to
64      !!            zero (zCflxemp=0) to calculate the correction velocity. So
65      !!            it will only balance the flux through open boundaries.
66      !!            (set volemp to 0 in tne namelist for this option)
67      !!
68      !!         2/ The volume is constant even with E-P flux. In this case
69      !!            the correction velocity must balance both the flux
70      !!            through open boundaries and the ones through the free
71      !!            surface.
72      !!            (set volemp to 1 in tne namelist for this option)
73      !!
74      !! History :
75      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Original
76      !!----------------------------------------------------------------------------
77      !! * Arguments
78      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
79
80      !! * Local declarations
81      INTEGER ::   ji, jj, jk
82      REAL(wp) ::   zubtpecor
83      REAL(wp) ::   zCflxemp
84      REAL(wp) ::   ztransw, ztranse, ztransn, ztranss, ztranst
85      !!-----------------------------------------------------------------------------
86
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)'        '
89         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along OBC'
90         IF(lwp) WRITE(numout,*)'~~~~~~~'
91         IF(lwp) WRITE(numout,*)'        '
92      END IF 
93
94      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
95      ! ---------------------------------------------------------------------------
96 
97      zCflxemp = 0.e0
98
99      DO jj = 1, jpj
100         DO ji = 1, jpi
101            zCflxemp = zCflxemp + ( (emp(ji,jj)*obctmsk(ji,jj) )/rauw)*e1v(ji,jj)*e2u(ji,jj)
102         END DO
103      END DO
104
105# if defined key_mpp
106      CALL mpp_sum( zCflxemp )
107# endif
108
109      ! 2. Barotropic velocity for each open boundary
110      ! ---------------------------------------------
111
112      zubtpecor = 0.e0
113
114      ! ... West open boundary
115      IF( lpwestobc ) THEN
116
117         ! ... Total transport through the West OBC
118         DO ji = fs_niw0, fs_niw1 ! Vector opt.
119            DO jk = 1, jpkm1
120               DO jj = 1, jpj
121                  zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) &
122                                        * uwmsk(jj,jk)
123               END DO
124            END DO
125         END DO
126
127      END IF 
128
129      ! ... East open boundary
130      IF( lpeastobc ) THEN
131
132         ! ... Total transport through the East OBC
133         DO ji = fs_nie0, fs_nie1 ! Vector opt.
134            DO jk = 1, jpkm1
135               DO jj = 1, jpj
136                  zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) &
137                                        * uemsk(jj,jk)
138               END DO
139            END DO
140         END DO
141
142      END IF 
143
144      ! ... North open boundary
145      IF( lpnorthobc ) THEN
146
147         ! ... Total transport through the North OBC
148         DO jj = fs_njn0, fs_njn1 ! Vector opt.
149            DO jk = 1, jpkm1
150               DO ji = 1, jpi
151                  zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) &
152                                        * vnmsk(ji,jk)
153               END DO
154            END DO
155         END DO
156
157      END IF 
158
159      ! ... South open boundary
160      IF( lpsouthobc ) THEN
161
162         ! ... Total transport through the South OBC
163         DO jj = fs_njs0, fs_njs1 ! Vector opt.
164            DO jk = 1, jpkm1
165               DO ji = 1, jpi
166                  zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) &
167                                        * vsmsk(ji,jk)
168               END DO
169            END DO
170         END DO
171
172      END IF 
173
174# if defined key_mpp
175      CALL mpp_sum( zubtpecor )
176# endif
177
178      ! 3. The normal velocity correction
179      ! ---------------------------------
180
181      zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot)
182
183      IF( lwp .and. mod( kt, nwrite ) == 0) THEN
184         IF(lwp) WRITE(numout,*)'        '
185         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
186         IF(lwp) WRITE(numout,*)'~~~~~~~ '
187         IF(lwp) WRITE(numout,*)'        '
188         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
189         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',obcsurftot,'(m2)'
190         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
191         IF(lwp) WRITE(numout,*)'        '
192      END IF 
193
194      ! 4. Correction of the total velocity on each open
195      !    boundary torespect the mass flux conservation
196      ! -------------------------------------------------
197
198      ztransw = 0.e0
199      ztranse = 0.e0
200      ztransn = 0.e0
201      ztranss = 0.e0
202      ztranst = 0.e0
203
204      IF( lpwestobc ) THEN
205
206         ! ... correction of the west velocity
207         DO ji = fs_niw0, fs_niw1 ! Vector opt.
208            DO jk = 1, jpkm1
209               DO jj = 1, jpj
210                  ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk)
211                  ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk)
212               END DO
213            END DO
214         END DO
215
216# if defined key_mpp
217         CALL mpp_sum( ztransw )
218# endif
219
220         IF( lwp .and. mod( kt, nwrite ) == 0) THEN
221            IF(lwp) WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)'
222         END IF
223
224      END IF
225
226      IF( lpeastobc ) THEN
227
228         ! ... correction of the east velocity
229         DO ji = fs_nie0, fs_nie1 ! Vector opt.
230            DO jk = 1, jpkm1
231               DO jj = 1, jpj
232                  ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk)
233                  ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk)
234               END DO
235            END DO
236         END DO
237
238# if defined key_mpp
239         CALL mpp_sum( ztranse )
240# endif
241
242         IF( lwp .and. mod( kt, nwrite ) == 0) THEN
243            IF(lwp) WRITE(numout,*)'          East OB transport ztranse :', ztranse,'(m3/s)'
244         END IF
245
246      END IF
247
248      IF( lpnorthobc ) THEN
249
250         ! ... correction of the north velocity
251         DO jj = fs_njn0, fs_njn1 ! Vector opt.
252            DO jk = 1, jpkm1
253               DO ji =  1, jpi
254                  va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk)
255                  ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk)
256               END DO
257            END DO
258         END DO
259
260# if defined key_mpp
261         CALL mpp_sum( ztransn )
262# endif
263
264         IF( lwp .and. mod( kt, nwrite ) == 0) THEN
265            IF(lwp) WRITE(numout,*)'          North OB transport ztransn :', ztransn,'(m3/s)'
266         END IF
267
268      END IF
269
270      IF( lpsouthobc ) THEN
271
272         ! ... correction of the south velocity
273         DO jj = fs_njs0, fs_njs1 ! Vector opt.
274            DO jk = 1, jpkm1
275               DO ji =  1, jpi
276                  va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk)
277                  ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk)
278               END DO
279            END DO
280         END DO
281 
282# if defined key_mpp
283         CALL mpp_sum( ztranss )
284# endif
285
286         IF( lwp .and. mod( kt, nwrite ) == 0) THEN
287            IF(lwp) WRITE(numout,*)'          South OB transport ztranss :', ztranss,'(m3/s)'
288         END IF
289
290      END IF 
291
292      ! 5. Check the cumulate transport through OBC
293      !    once barotropic velocities corrected
294      ! -------------------------------------------
295
296      ztranst = ztransw - ztranse + ztranss - ztransn
297
298      IF( lwp .and. mod( kt, nwrite ) == 0) THEN
299         IF(lwp) WRITE(numout,*)'        '
300         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
301         IF(lwp) WRITE(numout,*)'        '
302      END IF
303
304   END SUBROUTINE obc_vol
305
306#else
307   !!---------------------------------------------------------------------------------
308   !!  Default option :                                                   Empty module
309   !!---------------------------------------------------------------------------------
310CONTAINS
311   SUBROUTINE obc_vol        ! Empty routine
312   END SUBROUTINE obc_vol
313#endif
314
315   !!=================================================================================
316END MODULE obcvol
Note: See TracBrowser for help on using the repository browser.