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 trunk/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMO/OPA_SRC/BDY/bdyvol.F90 @ 911

Last change on this file since 911 was 911, checked in by ctlod, 16 years ago

Implementation of the BDY package, see ticket: #126
  • Property svn:executable set to *
File size: 8.7 KB
Line 
1MODULE bdyvol
2   !!=================================================================================
3   !!                       ***  MODULE  bdyvol  ***
4   !! Ocean dynamic :  Volume constraint when unstructured boundary
5   !!                  and Free surface are used
6   !!=================================================================================
7#if   defined key_bdy   &&   defined key_dynspg_flt
8   !!---------------------------------------------------------------------------------
9   !!   'key_bdy'               and              unstructured open boundary conditions
10   !!   'key_dynspg_flt'                                  constant volume free surface
11   !!---------------------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE phycst          ! physical constants
16   USE bdy_oce         ! ocean open boundary conditions
17   USE lib_mpp         ! for mppsum
18   USE in_out_manager  ! I/O manager
19   USE sbc_oce         ! ocean surface boundary conditions
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Accessibility
25   PUBLIC bdy_vol        ! routine called by dynspg_flt.h90
26
27   !! * Substitutions
28#  include "domzgr_substitute.h90"
29   !!---------------------------------------------------------------------------------
30   !!   OPA 9.0 , LODYC-IPSL  (2003)
31   !!---------------------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE bdy_vol ( kt )
36      !!------------------------------------------------------------------------------
37      !!                      ***  ROUTINE bdyvol  ***
38      !!
39      !! ** Purpose :
40      !!      This routine is called in dynspg_flt to control
41      !!      the volume of the system. A correction velocity is calculated
42      !!      to correct the total transport through the unstructured OBC.
43      !!      The total depth used is constant (H0) to be consistent with the
44      !!      linear free surface coded in OPA 8.2
45      !!
46      !! ** Method : 
47      !!      The correction velocity (zubtpecor here) is defined calculating
48      !!      the total transport through all open boundaries (trans_bdy) minus
49      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral
50      !!      surface (bdysurftot) of the unstructured boundary.
51      !!
52      !!      zubtpecor = [trans_bdy - zCflxemp ]*(1./bdysurftot)
53      !!
54      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
55      !!                       over all the domain in m3/s at each time step.
56      !!
57      !!      zCflxemp < 0 when precipitation dominate
58      !!      zCflxemp > 0 when evaporation dominate
59      !!
60      !!      There are 2 options (user's desiderata):
61      !!
62      !!         1/ The volume changes according to E-P, this is the default
63      !!            option. In this case the cumulate E-P flux are setting to
64      !!            zero (zCflxemp=0) to calculate the correction velocity. So
65      !!            it will only balance the flux through open boundaries.
66      !!            (set volbdy to 0 in tne namelist for this option)
67      !!
68      !!         2/ The volume is constant even with E-P flux. In this case
69      !!            the correction velocity must balance both the flux
70      !!            through open boundaries and the ones through the free
71      !!            surface.
72      !!            (set volbdy to 1 in tne namelist for this option)
73      !!
74      !! History :
75      !!   8.5  !  05-01 (J. Chanut, A. Sellar) Original code
76      !!        !  06-01 (J. Chanut) Bug correction
77      !!----------------------------------------------------------------------------
78      !! * Arguments
79      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
80
81      !! * Local declarations
82      INTEGER ::   ji,jj,jb, jgrd, jk
83      REAL(wp) ::   zubtpecor
84      REAL(wp) ::   zCflxemp
85      REAL(wp) ::   ztranst
86      !!-----------------------------------------------------------------------------
87
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)'        '
90         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC'
91         IF(lwp) WRITE(numout,*)'~~~~~~~'
92         IF(lwp) WRITE(numout,*)'        '
93      END IF 
94
95      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
96      ! ---------------------------------------------------------------------------
97 
98      zCflxemp = 0.e0
99
100      DO jj = 1, jpj
101         DO ji = 1, jpi
102            zCflxemp = zCflxemp + ( (emp(ji,jj)*bdytmask(ji,jj)*tmask_i(ji,jj) )/rauw) &
103                                  *e1v(ji,jj)*e2u(ji,jj)
104         END DO
105      END DO
106      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain
107
108      ! 2. Barotropic velocity through the unstructured open boundary
109      ! -------------------------------------------------------------
110
111      zubtpecor = 0.e0
112
113      jgrd = 2 ! cumulate u component contribution first
114      DO jb = 1, nblenrim(jgrd)
115        DO jk = 1, jpkm1
116          zubtpecor = zubtpecor + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
117                                            * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &
118                                            * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk)
119        END DO
120      END DO
121
122      jgrd = 3 ! then add v component contribution
123       DO jb = 1, nblenrim(jgrd)
124        DO jk = 1, jpkm1
125          zubtpecor = zubtpecor + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
126                                            * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) &
127                                            * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
128        END DO
129      END DO
130
131      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
132
133
134      ! 3. The normal velocity correction
135      ! ---------------------------------
136
137      IF (volbdy==1) THEN
138        zubtpecor = (zubtpecor - zCflxemp)*(1./bdysurftot) 
139      ELSE
140        zubtpecor =  zubtpecor*(1./bdysurftot)
141      END IF
142
143      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
144         IF(lwp) WRITE(numout,*)'        '
145         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
146         IF(lwp) WRITE(numout,*)'~~~~~~~ '
147         IF(lwp) WRITE(numout,*)'        '
148         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
149         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',bdysurftot,'(m2)'
150         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
151         IF(lwp) WRITE(numout,*)'        '
152      END IF 
153
154      ! 4. Correction of the total velocity on the unstructured
155      !    boundary to respect the mass flux conservation
156      ! -------------------------------------------------------
157
158      ztranst = 0.e0
159
160      jgrd = 2 ! correct u component
161      DO jb = 1, nblenrim(jgrd)
162        DO jk = 1, jpkm1
163          ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) = ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) &
164                            -flagu(jb) * zubtpecor * umask(nbi(jb, jgrd), nbj(jb, jgrd), jk)
165          ztranst = ztranst + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
166                                        * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &
167                                        * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk)
168        END DO
169      END DO
170
171      jgrd = 3 ! correct v component
172       DO jb = 1, nblenrim(jgrd)
173        DO jk = 1, jpkm1
174          va(nbi(jb, jgrd), nbj(jb, jgrd), jk) = va(nbi(jb, jgrd), nbj(jb, jgrd), jk) &
175                            -flagv(jb) * zubtpecor * vmask(nbi(jb, jgrd), nbj(jb, jgrd), jk)
176          ztranst = ztranst + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
177                                        * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) &
178                                        * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk)
179        END DO
180      END DO
181     
182      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
183 
184      ! 5. Check the cumulate transport through unstructured OBC
185      !    once barotropic velocities corrected
186      ! --------------------------------------------------------
187
188      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
189         IF(lwp) WRITE(numout,*)'        '
190         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
191         IF(lwp) WRITE(numout,*)'        '
192      END IF
193
194   END SUBROUTINE bdy_vol
195
196#else
197   !!---------------------------------------------------------------------------------
198   !!  Default option :                                                   Empty module
199   !!---------------------------------------------------------------------------------
200CONTAINS
201   SUBROUTINE bdy_vol        ! Empty routine
202   END SUBROUTINE bdy_vol
203#endif
204
205   !!=================================================================================
206END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.