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.
bdyvol.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

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, 5 years 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
RevLine 
[911]1MODULE bdyvol
[1125]2   !!======================================================================
[911]3   !!                       ***  MODULE  bdyvol  ***
4   !! Ocean dynamic :  Volume constraint when unstructured boundary
[3294]5   !!                  and filtered free surface are used
[1125]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
[3294]10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
[10481]11   !!            4.0  !  2019-01  (P. Mathiot) adapted to time splitting     
[1125]12   !!----------------------------------------------------------------------
[6140]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
[5836]19   !
[6140]20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! for mppsum
22   USE lib_fortran    ! Fortran routines library
[911]23
24   IMPLICIT NONE
25   PRIVATE
26
[10481]27   PUBLIC   bdy_vol2d    ! called by dynspg_ts
[911]28
[1125]29   !!----------------------------------------------------------------------
[9598]30   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1146]31   !! $Id$
[10068]32   !! Software governed by the CeCILL license (see ./LICENSE)
[1125]33   !!----------------------------------------------------------------------
[911]34CONTAINS
35
[10481]36   SUBROUTINE bdy_vol2d( kt, kc, pua2d, pva2d, phu, phv )
[1125]37      !!----------------------------------------------------------------------
[911]38      !!                      ***  ROUTINE bdyvol  ***
39      !!
[5930]40      !! ** Purpose :   This routine controls the volume of the system.
[6140]41      !!      A correction velocity is calculated to correct the total transport
42      !!      through the unstructured OBC.
[911]43      !!
[1125]44      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
[911]45      !!      the total transport through all open boundaries (trans_bdy) minus
[1125]46      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral
[911]47      !!      surface (bdysurftot) of the unstructured boundary.
[1125]48      !!         zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot)
49      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
[911]50      !!                       over all the domain in m3/s at each time step.
[1125]51      !!      z_cflxemp < 0 when precipitation dominate
52      !!      z_cflxemp > 0 when evaporation dominate
[911]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
[1125]57      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
[911]58      !!            it will only balance the flux through open boundaries.
[2528]59      !!            (set nn_volctl to 0 in tne namelist for this option)
[911]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.
[2528]64      !!            (set nn_volctl to 1 in tne namelist for this option)
[1125]65      !!----------------------------------------------------------------------
[10481]66      INTEGER, INTENT(in) ::   kt, kc   ! ocean time-step index, cycle time-step
[6140]67      !
[1125]68      INTEGER  ::   ji, jj, jk, jb, jgrd
[3294]69      INTEGER  ::   ib_bdy, ii, ij
[10481]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
[3294]74      TYPE(OBC_INDEX), POINTER :: idx
[911]75      !!-----------------------------------------------------------------------------
[5836]76      !
[1125]77      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
78      ! -----------------------------------------------------------------------
[10481]79      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
[911]80
[10481]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
[2528]91      ! Transport through the unstructured open boundary
92      ! ------------------------------------------------
[5836]93      zubtpecor = 0._wp
[3294]94      DO ib_bdy = 1, nb_bdy
95         idx => idx_bdy(ib_bdy)
[5836]96         !
[3294]97         jgrd = 2                               ! cumulate u component contribution first
98         DO jb = 1, idx%nblenrim(jgrd)
[10481]99            ii = idx%nbi(jb,jgrd)
100            ij = idx%nbj(jb,jgrd)
[11048]101            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? check tmask_i definition...
[10481]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)
[1125]103         END DO
[3294]104         jgrd = 3                               ! then add v component contribution
105         DO jb = 1, idx%nblenrim(jgrd)
[10481]106            ii = idx%nbi(jb,jgrd)
107            ij = idx%nbj(jb,jgrd)
[11048]108            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
[10481]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)
[1125]110         END DO
[5836]111         !
[911]112      END DO
[10481]113      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
[911]114
[1125]115      ! The normal velocity correction
116      ! ------------------------------
[10481]117      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  ! maybe should be apply only once at the end
[6140]118      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
[911]119      END IF
120
[1125]121      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
122      ! -------------------------------------------------------------
[3294]123      DO ib_bdy = 1, nb_bdy
124         idx => idx_bdy(ib_bdy)
[5836]125         !
[3294]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)
[11048]130               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
[10481]131               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij)
[1125]132         END DO
[3294]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)
[11048]137               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
[10481]138               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1)
[1125]139         END DO
[5836]140         !
[911]141      END DO
[10481]142      !
[1125]143      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
144      ! ------------------------------------------------------
[10481]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)
[11048]156                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
[10481]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)
[11048]163                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
[10481]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
[1125]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)'
[911]178      END IF 
[1125]179      !
[10481]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)
[11048]203            IF( nbi == 1 .OR. nbi == jpi .OR. nbj == 1 .OR. nbj == jpj )  CYCLE
[10481]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)
[11048]216            IF( nbi == 1 .OR. nbi == jpi .OR. nbj == 1 .OR. nbj == jpj )  CYCLE
[10481]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
[5836]223      !
[10481]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
[1125]228   !!======================================================================
[911]229END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.