source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90 @ 11048

Last change on this file since 11048 was 11048, checked in by girrmann, 19 months ago

dev_r10984_HPC-13 : Step 1, boundary is now detected all over the local domain, this does not change the result. Improve bdy treatment for bdy_rnf in bdytra.F90, this changes the result when keyword runoff is specified in namelist

  • Property svn:keywords set to Id
File size: 11.6 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
71      REAL(wp), SAVE :: z_cflxemp                                  ! cumulated emp flux
72      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d  ! Barotropic velocities
73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv         ! Ocean depth at U- and V-points
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            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? check tmask_i definition...
102            zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
103         END DO
104         jgrd = 3                               ! then add v component contribution
105         DO jb = 1, idx%nblenrim(jgrd)
106            ii = idx%nbi(jb,jgrd)
107            ij = idx%nbj(jb,jgrd)
108            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
109            zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
110         END DO
111         !
112      END DO
113      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
114
115      ! The normal velocity correction
116      ! ------------------------------
117      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  ! maybe should be apply only once at the end
118      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
119      END IF
120
121      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
122      ! -------------------------------------------------------------
123      DO ib_bdy = 1, nb_bdy
124         idx => idx_bdy(ib_bdy)
125         !
126         jgrd = 2                               ! correct u component
127         DO jb = 1, idx%nblenrim(jgrd)
128               ii = idx%nbi(jb,jgrd)
129               ij = idx%nbj(jb,jgrd)
130               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
131               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij)
132         END DO
133         jgrd = 3                              ! correct v component
134         DO jb = 1, idx%nblenrim(jgrd)
135               ii = idx%nbi(jb,jgrd)
136               ij = idx%nbj(jb,jgrd)
137               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
138               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1)
139         END DO
140         !
141      END DO
142      !
143      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
144      ! ------------------------------------------------------
145      IF( MOD( kt, nwrite ) == 0 .AND. ( kc == 1 ) ) THEN
146         !
147         ! compute residual transport across boundary
148         ztranst = 0._wp
149         DO ib_bdy = 1, nb_bdy
150            idx => idx_bdy(ib_bdy)
151            !
152            jgrd = 2                               ! correct u component
153            DO jb = 1, idx%nblenrim(jgrd)
154                  ii = idx%nbi(jb,jgrd)
155                  ij = idx%nbj(jb,jgrd)
156                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
157                  ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
158            END DO
159            jgrd = 3                              ! correct v component
160            DO jb = 1, idx%nblenrim(jgrd)
161                  ii = idx%nbi(jb,jgrd)
162                  ij = idx%nbj(jb,jgrd)
163                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
164                  ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
165            END DO
166            !
167         END DO
168         IF( lk_mpp )   CALL mpp_sum('bdyvol', ztranst )   ! sum over the global domain
169
170
171         IF(lwp) WRITE(numout,*)
172         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
173         IF(lwp) WRITE(numout,*)'~~~~~~~ '
174         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
175         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)'
176         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
177         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
178      END IF 
179      !
180   END SUBROUTINE bdy_vol2d
181   !
182   REAL(wp) FUNCTION bdy_segs_surf(phu, phv)
183      !!----------------------------------------------------------------------
184      !!                 ***  ROUTINE bdy_ctl_seg  ***
185      !!
186      !! ** Purpose :   Compute total lateral surface for volume correction
187      !!
188      !!----------------------------------------------------------------------
189
190      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv ! water column thickness at U and V points
191      INTEGER            ::  igrd, ib_bdy, ib                ! loop indexes
192      INTEGER , POINTER  ::  nbi, nbj                        ! short cuts
193      REAL(wp), POINTER  ::  zflagu, zflagv                  !    -   -
194
195      ! Compute total lateral surface for volume correction:
196      ! ----------------------------------------------------
197      bdy_segs_surf = 0._wp
198      igrd = 2      ! Lateral surface at U-points
199      DO ib_bdy = 1, nb_bdy
200         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
201            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
202            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
203            IF( nbi == 1 .OR. nbi == jpi .OR. nbj == 1 .OR. nbj == jpj )  CYCLE
204            zflagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
205            bdy_segs_surf = bdy_segs_surf + phu(nbi, nbj)                              &
206               &                            * e2u(nbi, nbj) * ABS( zflagu )            &
207               &                            * tmask_i(nbi, nbj) * tmask_i(nbi+1, nbj)
208         END DO
209      END DO
210
211      igrd=3 ! Add lateral surface at V-points
212      DO ib_bdy = 1, nb_bdy
213         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
214            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
215            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
216            IF( nbi == 1 .OR. nbi == jpi .OR. nbj == 1 .OR. nbj == jpj )  CYCLE
217            zflagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
218            bdy_segs_surf = bdy_segs_surf + phv(nbi, nbj)                              &
219               &                            * e1v(nbi, nbj) * ABS( zflagv )            &
220               &                            * tmask_i(nbi, nbj) * tmask_i(nbi, nbj+1)
221         END DO
222      END DO
223      !
224      ! redirect the time to bdyvol as this variable is only used by bdyvol
225      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', bdy_segs_surf ) ! sum over the global domain
226      !
227   END FUNCTION bdy_segs_surf
228   !!======================================================================
229END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.