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.
bdytra.F90 in branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90 @ 4007

Last change on this file since 4007 was 4007, checked in by davestorkey, 11 years ago
  1. Bug fixes for flagu/flagv calculation in bdyini.F90.
  2. Introduce masking of derivatives in radiation velocity calculation in bdylib.F90
  3. Change relaxation term in Orlanski implementation to explicit timestepping in bdylib.F90.
  4. Remove bdyfmask (redundant).
  • Property svn:keywords set to Id
File size: 14.5 KB
Line 
1MODULE bdytra
2   !!======================================================================
3   !!                       ***  MODULE  bdytra  ***
4   !! Ocean tracers:   Apply boundary conditions for tracers
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
8   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
9   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
10   !!----------------------------------------------------------------------
11#if defined key_bdy
12   !!----------------------------------------------------------------------
13   !!   'key_bdy'                     Unstructured Open Boundary Conditions
14   !!----------------------------------------------------------------------
15   !!   bdy_tra            : Apply open boundary conditions to T and S
16   !!   bdy_tra_frs        : Apply Flow Relaxation Scheme
17   !!----------------------------------------------------------------------
18   USE timing          ! Timing
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
21   USE bdy_oce         ! ocean open boundary conditions
22   USE bdylib          ! for orlanski library routines
23   USE bdydta, ONLY:   bf
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE in_out_manager  ! I/O manager
26
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC bdy_tra      ! routine called in tranxt.F90
32   PUBLIC bdy_tra_dmp  ! routine called in step.F90
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE bdy_tra( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  SUBROUTINE bdy_tra  ***
44      !!
45      !! ** Purpose : - Apply open boundary conditions for temperature and salinity
46      !!
47      !!----------------------------------------------------------------------
48      INTEGER, INTENT( in ) :: kt     ! Main time step counter
49      !!
50      INTEGER               :: ib_bdy ! Loop index
51
52      DO ib_bdy=1, nb_bdy
53
54         SELECT CASE( cn_tra(ib_bdy) )
55         CASE('none')
56            CYCLE
57         CASE('frs')
58            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
59         CASE('specified')
60            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
61         CASE('neumann')
62            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
63         CASE('orlanski')
64            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.false. )
65         CASE('orlanski_npo')
66            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.true. )
67         CASE('runoff')
68            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
69         CASE DEFAULT
70            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
71         END SELECT
72         ! Boundary points should be updated
73         CALL lbc_bdy_lnk( tsa(:,:,:,jp_tem), 'T', 1., ib_bdy )
74         CALL lbc_bdy_lnk( tsa(:,:,:,jp_sal), 'T', 1., ib_bdy )
75      ENDDO
76      !
77
78   END SUBROUTINE bdy_tra
79
80   SUBROUTINE bdy_tra_frs( idx, dta, kt )
81      !!----------------------------------------------------------------------
82      !!                 ***  SUBROUTINE bdy_tra_frs  ***
83      !!                   
84      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
85      !!
86      !! Reference : Engedahl H., 1995, Tellus, 365-382.
87      !!----------------------------------------------------------------------
88      INTEGER,         INTENT(in) ::   kt
89      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
90      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
91      !!
92      REAL(wp) ::   zwgt           ! boundary weight
93      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
94      INTEGER  ::   ii, ij         ! 2D addresses
95      !!----------------------------------------------------------------------
96      !
97      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_frs')
98      !
99      igrd = 1                       ! Everything is at T-points here
100      DO ib = 1, idx%nblen(igrd)
101         DO ik = 1, jpkm1
102            ii = idx%nbi(ib,igrd)
103            ij = idx%nbj(ib,igrd)
104            zwgt = idx%nbw(ib,igrd)
105            tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta%tem(ib,ik) - tsa(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik)         
106            tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta%sal(ib,ik) - tsa(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik)
107         END DO
108      END DO 
109      !
110      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
111      !
112      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_frs')
113      !
114   END SUBROUTINE bdy_tra_frs
115 
116   SUBROUTINE bdy_tra_spe( idx, dta, kt )
117      !!----------------------------------------------------------------------
118      !!                 ***  SUBROUTINE bdy_tra_frs  ***
119      !!                   
120      !! ** Purpose : Apply a specified value for tracers at open boundaries.
121      !!
122      !!----------------------------------------------------------------------
123      INTEGER,         INTENT(in) ::   kt
124      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
125      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
126      !!
127      REAL(wp) ::   zwgt           ! boundary weight
128      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
129      INTEGER  ::   ii, ij         ! 2D addresses
130      !!----------------------------------------------------------------------
131      !
132      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_spe')
133      !
134      igrd = 1                       ! Everything is at T-points here
135      DO ib = 1, idx%nblenrim(igrd)
136         ii = idx%nbi(ib,igrd)
137         ij = idx%nbj(ib,igrd)
138         DO ik = 1, jpkm1
139            tsa(ii,ij,ik,jp_tem) = dta%tem(ib,ik) * tmask(ii,ij,ik)
140            tsa(ii,ij,ik,jp_sal) = dta%sal(ib,ik) * tmask(ii,ij,ik)
141         END DO
142      END DO
143      !
144      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
145      !
146      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe')
147      !
148   END SUBROUTINE bdy_tra_spe
149
150   SUBROUTINE bdy_tra_nmn( idx, dta, kt )
151      !!----------------------------------------------------------------------
152      !!                 ***  SUBROUTINE bdy_tra_nmn  ***
153      !!                   
154      !! ** Purpose : Duplicate the value for tracers at open boundaries.
155      !!
156      !!----------------------------------------------------------------------
157      INTEGER,         INTENT(in) ::   kt
158      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
159      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
160      !!
161      REAL(wp) ::   zwgt           ! boundary weight
162      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
163      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
164      !!----------------------------------------------------------------------
165      !
166      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn')
167      !
168      igrd = 1                       ! Everything is at T-points here
169      DO ib = 1, idx%nblenrim(igrd)
170         ii = idx%nbi(ib,igrd)
171         ij = idx%nbj(ib,igrd)
172         DO ik = 1, jpkm1
173            ! search the sense of the gradient
174            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
175            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
176            IF ( zcoef1+zcoef2 == 0) THEN
177               ! corner
178               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
179               tsa(ii,ij,ik,jp_tem) = tsa(ii-1,ij  ,ik,jp_tem) * tmask(ii-1,ij  ,ik) + &
180                 &                    tsa(ii+1,ij  ,ik,jp_tem) * tmask(ii+1,ij  ,ik) + &
181                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
182                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
183               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
184               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
185                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
186                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
187                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
188               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
189            ELSE
190               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
191               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
192               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
193               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
194            ENDIF
195         END DO
196      END DO
197      !
198      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
199      !
200      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn')
201      !
202   END SUBROUTINE bdy_tra_nmn
203 
204
205   SUBROUTINE bdy_tra_orlanski( idx, dta, ll_npo )
206      !!----------------------------------------------------------------------
207      !!                 ***  SUBROUTINE bdy_tra_orlanski  ***
208      !!             
209      !!              - Apply Orlanski radiation to temperature and salinity.
210      !!              - Wrapper routine for bdy_orlanski_3d
211      !!
212      !!
213      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
214      !!----------------------------------------------------------------------
215      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
216      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
217      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
218
219      INTEGER  ::   igrd                                    ! grid index
220      !!----------------------------------------------------------------------
221
222      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_orlanski')
223      !
224      igrd = 1      ! Orlanski bc on temperature;
225      !           
226      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_tem), tsa(:,:,:,jp_tem), dta%tem, ll_npo )
227
228      igrd = 1      ! Orlanski bc on salinity;
229      ! 
230      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_sal), tsa(:,:,:,jp_sal), dta%sal, ll_npo )
231      !
232      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_orlanski')
233      !
234
235   END SUBROUTINE bdy_tra_orlanski
236
237
238   SUBROUTINE bdy_tra_rnf( idx, dta, kt )
239      !!----------------------------------------------------------------------
240      !!                 ***  SUBROUTINE bdy_tra_rnf  ***
241      !!                   
242      !! ** Purpose : Apply the runoff values for tracers at open boundaries:
243      !!                  - specified to 0.1 PSU for the salinity
244      !!                  - duplicate the value for the temperature
245      !!
246      !!----------------------------------------------------------------------
247      INTEGER,         INTENT(in) ::   kt
248      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
249      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
250      !!
251      REAL(wp) ::   zwgt           ! boundary weight
252      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
253      INTEGER  ::   ii, ij, ip, jp ! 2D addresses
254      !!----------------------------------------------------------------------
255      !
256      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf')
257      !
258      igrd = 1                       ! Everything is at T-points here
259      DO ib = 1, idx%nblenrim(igrd)
260         ii = idx%nbi(ib,igrd)
261         ij = idx%nbj(ib,igrd)
262         DO ik = 1, jpkm1
263            ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
264            jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
265            tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii,ij,ik)
266            tsa(ii,ij,ik,jp_sal) =                        0.1 * tmask(ii,ij,ik)
267         END DO
268      END DO
269      !
270      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
271      !
272      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf')
273      !
274   END SUBROUTINE bdy_tra_rnf
275
276   SUBROUTINE bdy_tra_dmp( kt )
277      !!----------------------------------------------------------------------
278      !!                 ***  SUBROUTINE bdy_tra_dmp  ***
279      !!                   
280      !! ** Purpose : Apply damping for tracers at open boundaries.
281      !!
282      !!----------------------------------------------------------------------
283      INTEGER,         INTENT(in) ::   kt
284      !!
285      REAL(wp) ::   zwgt           ! boundary weight
286      REAL(wp) ::   zta, zsa, ztime
287      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
288      INTEGER  ::   ii, ij         ! 2D addresses
289      INTEGER  ::   ib_bdy         ! Loop index
290      !!----------------------------------------------------------------------
291      !
292      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp')
293      !
294      DO ib_bdy=1, nb_bdy
295         IF ( ln_tra_dmp(ib_bdy) ) THEN
296            igrd = 1                       ! Everything is at T-points here
297            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
298               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
299               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
300               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
301               DO ik = 1, jpkm1
302                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik)
303                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik)
304                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta
305                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa
306               END DO
307            END DO
308         ENDIF
309      ENDDO
310      !
311      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp')
312      !
313   END SUBROUTINE bdy_tra_dmp
314 
315#else
316   !!----------------------------------------------------------------------
317   !!   Dummy module                   NO Unstruct Open Boundary Conditions
318   !!----------------------------------------------------------------------
319CONTAINS
320   SUBROUTINE bdy_tra(kt)      ! Empty routine
321      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt
322   END SUBROUTINE bdy_tra
323
324   SUBROUTINE bdy_tra_dmp(kt)      ! Empty routine
325      WRITE(*,*) 'bdy_tra_dmp: You should not have seen this print! error?', kt
326   END SUBROUTINE bdy_tra_dmp
327
328#endif
329
330   !!======================================================================
331END MODULE bdytra
Note: See TracBrowser for help on using the repository browser.