source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyvol.F90 @ 10455

Last change on this file since 10455 was 10455, checked in by mathiot, 22 months ago

suggestion for ticket #2200 (bdyvol)

  • Property svn:keywords set to Id
File size: 10.7 KB
Line 
1MODULE bdyvol
2   !!======================================================================
3   !!                       ***  MODULE  bdyvol  ***
4   !! Ocean dynamic :  Volume constraint when unstructured boundary
5   !!                  and filtered free surface are used
6   !!======================================================================
7   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
8   !!             -   !  2006-01  (J. Chanut) Bug correction
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
11   !!            4.0  !  2019-01  (P. Mathiot) adapted to time splitting     
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE bdy_oce        ! ocean open boundary conditions
15   USE sbc_oce        ! ocean surface boundary conditions
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE sbcisf         ! ice shelf
19   !
20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! for mppsum
22   USE lib_fortran    ! Fortran routines library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   bdy_vol2d    ! called by dynspg_ts
28
29   !!----------------------------------------------------------------------
30   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
31   !! $Id$
32   !! Software governed by the CeCILL license (see ./LICENSE)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE bdy_vol2d( kt, kc, pua2d, pva2d, phu, phv )
37      !!----------------------------------------------------------------------
38      !!                      ***  ROUTINE bdyvol  ***
39      !!
40      !! ** Purpose :   This routine controls the volume of the system.
41      !!      A correction velocity is calculated to correct the total transport
42      !!      through the unstructured OBC.
43      !!
44      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
45      !!      the total transport through all open boundaries (trans_bdy) minus
46      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral
47      !!      surface (bdysurftot) of the unstructured boundary.
48      !!         zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot)
49      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
50      !!                       over all the domain in m3/s at each time step.
51      !!      z_cflxemp < 0 when precipitation dominate
52      !!      z_cflxemp > 0 when evaporation dominate
53      !!
54      !!      There are 2 options (user's desiderata):
55      !!         1/ The volume changes according to E-P, this is the default
56      !!            option. In this case the cumulate E-P flux are setting to
57      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
58      !!            it will only balance the flux through open boundaries.
59      !!            (set nn_volctl to 0 in tne namelist for this option)
60      !!         2/ The volume is constant even with E-P flux. In this case
61      !!            the correction velocity must balance both the flux
62      !!            through open boundaries and the ones through the free
63      !!            surface.
64      !!            (set nn_volctl to 1 in tne namelist for this option)
65      !!----------------------------------------------------------------------
66      INTEGER, INTENT(in) ::   kt, kc   ! ocean time-step index, cycle time step
67      !
68      INTEGER  ::   ji, jj, jk, jb, jgrd
69      INTEGER  ::   ib_bdy, ii, ij
70      REAL(wp) ::   zubtpecor, ztranst, zeh
71      REAL(wp), SAVE :: z_cflxemp   
72      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv
74      TYPE(OBC_INDEX), POINTER :: idx
75      !!-----------------------------------------------------------------------------
76      !
77      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
78      ! -----------------------------------------------------------------------
79      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
80
81      ! Compute bdy surface each cycle if non linear free surface
82      ! ---------------------------------------------------------
83      IF ( .NOT. ln_linssh ) THEN
84         ! compute area each time step
85         bdysurftot = bdy_segs_surf( phu, phv )
86      ELSE
87         ! compute area only the first time
88         IF ( ( kt == nit000 ) .AND. ( kc == 1 ) ) bdysurftot = bdy_segs_surf( phu, phv )
89      END IF
90
91      ! Transport through the unstructured open boundary
92      ! ------------------------------------------------
93      zubtpecor = 0._wp
94      DO ib_bdy = 1, nb_bdy
95         idx => idx_bdy(ib_bdy)
96         !
97         jgrd = 2                               ! cumulate u component contribution first
98         DO jb = 1, idx%nblenrim(jgrd)
99            ii = idx%nbi(jb,jgrd)
100            ij = idx%nbj(jb,jgrd)
101            zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
102         END DO
103         jgrd = 3                               ! then add v component contribution
104         DO jb = 1, idx%nblenrim(jgrd)
105            ii = idx%nbi(jb,jgrd)
106            ij = idx%nbj(jb,jgrd)
107            zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
108         END DO
109         !
110      END DO
111      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
112
113      ! The normal velocity correction
114      ! ------------------------------
115      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  ! maybe should be apply only once at the end
116      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
117      END IF
118
119      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
120      ! -------------------------------------------------------------
121      DO ib_bdy = 1, nb_bdy
122         idx => idx_bdy(ib_bdy)
123         !
124         jgrd = 2                               ! correct u component
125         DO jb = 1, idx%nblenrim(jgrd)
126               ii = idx%nbi(jb,jgrd)
127               ij = idx%nbj(jb,jgrd)
128               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij)
129         END DO
130         jgrd = 3                              ! correct v component
131         DO jb = 1, idx%nblenrim(jgrd)
132               ii = idx%nbi(jb,jgrd)
133               ij = idx%nbj(jb,jgrd)
134               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1)
135         END DO
136         !
137      END DO
138      !
139      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
140      ! ------------------------------------------------------
141      IF( MOD( kt, nwrite ) == 0 .AND. ( kc == 1 ) ) THEN
142         !
143         ! compute residual transport across boundary
144         ztranst = 0._wp
145         DO ib_bdy = 1, nb_bdy
146            idx => idx_bdy(ib_bdy)
147            !
148            jgrd = 2                               ! correct u component
149            DO jb = 1, idx%nblenrim(jgrd)
150                  ii = idx%nbi(jb,jgrd)
151                  ij = idx%nbj(jb,jgrd)
152                  ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
153            END DO
154            jgrd = 3                              ! correct v component
155            DO jb = 1, idx%nblenrim(jgrd)
156                  ii = idx%nbi(jb,jgrd)
157                  ij = idx%nbj(jb,jgrd)
158                  ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
159            END DO
160            !
161         END DO
162         IF( lk_mpp )   CALL mpp_sum('bdyvol', ztranst )   ! sum over the global domain
163
164
165         IF(lwp) WRITE(numout,*)
166         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
167         IF(lwp) WRITE(numout,*)'~~~~~~~ '
168         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
169         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)'
170         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
171         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
172      END IF 
173      !
174   END SUBROUTINE bdy_vol2d
175   !
176   REAL(wp) FUNCTION bdy_segs_surf(phu, phv)
177      !!----------------------------------------------------------------------
178      !!                 ***  ROUTINE bdy_ctl_seg  ***
179      !!
180      !! ** Purpose :   Compute total lateral surface for volume correction
181      !!
182      !!----------------------------------------------------------------------
183
184      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv ! water column thickness at U and V points
185      INTEGER            ::  igrd, ib_bdy, ib                ! loop indexes
186      INTEGER , POINTER  ::  nbi, nbj                        ! short cuts
187      REAL(wp), POINTER  ::  zflagu, zflagv                  !    -   -
188
189      ! Compute total lateral surface for volume correction:
190      ! ----------------------------------------------------
191      bdy_segs_surf = 0._wp
192      igrd = 2      ! Lateral surface at U-points
193      DO ib_bdy = 1, nb_bdy
194         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
195            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
196            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
197            zflagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
198            bdy_segs_surf = bdy_segs_surf + phu(nbi, nbj)                              &
199               &                            * e2u(nbi, nbj) * ABS( zflagu )            &
200               &                            * tmask_i(nbi, nbj) * tmask_i(nbi+1, nbj)
201         END DO
202      END DO
203
204      igrd=3 ! Add lateral surface at V-points
205      DO ib_bdy = 1, nb_bdy
206         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
207            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
208            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
209            zflagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
210            bdy_segs_surf = bdy_segs_surf + phv(nbi, nbj)                              &
211               &                            * e1v(nbi, nbj) * ABS( zflagv )            &
212               &                            * tmask_i(nbi, nbj) * tmask_i(nbi, nbj+1)
213         END DO
214      END DO
215      !
216      ! redirect the time to bdyvol as this variable is only used by bdyvol
217      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', bdy_segs_surf ) ! sum over the global domain
218      !
219   END FUNCTION bdy_segs_surf
220   !!======================================================================
221END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.