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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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