source: trunk/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90 @ 23

Last change on this file since 23 was 23, checked in by rblod, 12 years ago

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