source: NEMO/trunk/src/OCE/BDY/bdyvol.F90 @ 10425

Last change on this file since 10425 was 10425, checked in by smasson, 22 months ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 7.9 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   !!----------------------------------------------------------------------
12   USE oce            ! ocean dynamics and tracers
13   USE bdy_oce        ! ocean open boundary conditions
14   USE sbc_oce        ! ocean surface boundary conditions
15   USE dom_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   USE sbcisf         ! ice shelf
18   !
19   USE in_out_manager ! I/O manager
20   USE lib_mpp        ! for mppsum
21   USE lib_fortran    ! Fortran routines library
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   bdy_vol    ! called by ???
27
28   !!----------------------------------------------------------------------
29   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
30   !! $Id$
31   !! Software governed by the CeCILL license (see ./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   SUBROUTINE bdy_vol( kt )
36      !!----------------------------------------------------------------------
37      !!                      ***  ROUTINE bdyvol  ***
38      !!
39      !! ** Purpose :   This routine controls the volume of the system.
40      !!      A correction velocity is calculated to correct the total transport
41      !!      through the unstructured OBC.
42      !!      The total depth used is constant (H0) to be consistent with the
43      !!      linear free surface coded in OPA 8.2    <<<=== !!gm  ???? true ????
44      !!
45      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
46      !!      the total transport through all open boundaries (trans_bdy) minus
47      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral
48      !!      surface (bdysurftot) of the unstructured boundary.
49      !!         zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot)
50      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
51      !!                       over all the domain in m3/s at each time step.
52      !!      z_cflxemp < 0 when precipitation dominate
53      !!      z_cflxemp > 0 when evaporation dominate
54      !!
55      !!      There are 2 options (user's desiderata):
56      !!         1/ The volume changes according to E-P, this is the default
57      !!            option. In this case the cumulate E-P flux are setting to
58      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
59      !!            it will only balance the flux through open boundaries.
60      !!            (set nn_volctl to 0 in tne namelist for this option)
61      !!         2/ The volume is constant even with E-P flux. In this case
62      !!            the correction velocity must balance both the flux
63      !!            through open boundaries and the ones through the free
64      !!            surface.
65      !!            (set nn_volctl to 1 in tne namelist for this option)
66      !!----------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
68      !
69      INTEGER  ::   ji, jj, jk, jb, jgrd
70      INTEGER  ::   ib_bdy, ii, ij
71      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
72      TYPE(OBC_INDEX), POINTER :: idx
73      !!-----------------------------------------------------------------------------
74      !
75      IF( ln_vol ) THEN
76      !
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC'
80         IF(lwp) WRITE(numout,*)'~~~~~~~'
81      END IF 
82
83      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
84      ! -----------------------------------------------------------------------
85!!gm replace these lines :
86      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
87      CALL mpp_sum( 'bdyvol', z_cflxemp )     ! sum over the global domain
88!!gm   by :
89!!gm      z_cflxemp = glob_sum(  'bdyvol', ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
90!!gm
91
92      ! Transport through the unstructured open boundary
93      ! ------------------------------------------------
94      zubtpecor = 0._wp
95      DO ib_bdy = 1, nb_bdy
96         idx => idx_bdy(ib_bdy)
97         !
98         jgrd = 2                               ! cumulate u component contribution first
99         DO jb = 1, idx%nblenrim(jgrd)
100            DO jk = 1, jpkm1
101               ii = idx%nbi(jb,jgrd)
102               ij = idx%nbj(jb,jgrd)
103               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
104            END DO
105         END DO
106         jgrd = 3                               ! then add v component contribution
107         DO jb = 1, idx%nblenrim(jgrd)
108            DO jk = 1, jpkm1
109               ii = idx%nbi(jb,jgrd)
110               ij = idx%nbj(jb,jgrd)
111               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * e3v_n(ii,ij,jk) 
112            END DO
113         END DO
114         !
115      END DO
116      CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
117
118      ! The normal velocity correction
119      ! ------------------------------
120      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot 
121      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
122      END IF
123
124      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
125      ! -------------------------------------------------------------
126      ztranst = 0._wp
127      DO ib_bdy = 1, nb_bdy
128         idx => idx_bdy(ib_bdy)
129         !
130         jgrd = 2                               ! correct u component
131         DO jb = 1, idx%nblenrim(jgrd)
132            DO jk = 1, jpkm1
133               ii = idx%nbi(jb,jgrd)
134               ij = idx%nbj(jb,jgrd)
135               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk)
136               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
137            END DO
138         END DO
139         jgrd = 3                              ! correct v component
140         DO jb = 1, idx%nblenrim(jgrd)
141            DO jk = 1, jpkm1
142               ii = idx%nbi(jb,jgrd)
143               ij = idx%nbj(jb,jgrd)
144               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk)
145               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * e3v_n(ii,ij,jk)
146            END DO
147         END DO
148         !
149      END DO
150      CALL mpp_sum( 'bdyvol', ztranst )   ! sum over the global domain
151 
152      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
153      ! ------------------------------------------------------
154      IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN
155         IF(lwp) WRITE(numout,*)
156         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
157         IF(lwp) WRITE(numout,*)'~~~~~~~ '
158         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
159         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)'
160         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
161         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
162      END IF 
163      !
164      END IF ! ln_vol
165      !
166   END SUBROUTINE bdy_vol
167
168   !!======================================================================
169END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.