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.
trcbdy.F90 in branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90 @ 11137

Last change on this file since 11137 was 11137, checked in by jcastill, 5 years ago

Add missing files
File size: 13.8 KB
Line 
1MODULE trcbdy
2   !!======================================================================
3   !!                       ***  MODULE  bdytrc  ***
4   !! Ocean tracers:   Apply boundary conditions for tracers in TOP component
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   !!            3.6  !  2015     (T. Lovato) Adapt BDY for tracers in TOP component
11   !!----------------------------------------------------------------------
12#if defined key_bdy && key_top
13   !!----------------------------------------------------------------------
14   !!   'key_bdy'                     Unstructured Open Boundary Conditions
15   !!----------------------------------------------------------------------
16   !!   trc_bdy            : Apply open boundary conditions to T and S
17   !!   trc_bdy_frs        : Apply Flow Relaxation Scheme
18   !!----------------------------------------------------------------------
19   USE timing                       ! Timing
20   USE oce_trc                      ! ocean dynamics and tracers variables
21   USE par_trc
22   USE trc                          ! ocean space and time domain variables
23   USE bdylib                       ! for orlanski library routines
24   USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
25   USE in_out_manager               ! I/O manager
26   USE bdy_oce, only: idx_bdy, OBC_INDEX, BDYTMASK, lk_bdy       ! ocean open boundary conditions
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC trc_bdy      ! routine called in trcnxt.F90
32   PUBLIC trc_bdy_dmp  ! routine called in trcstp.F90
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
36   !! $Id$
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE trc_bdy( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  SUBROUTINE trc_bdy  ***
44      !!
45      !! ** Purpose : - Apply open boundary conditions for tracers in TOP component
46      !!                and scale the tracer data
47      !!
48      !!----------------------------------------------------------------------
49      INTEGER, INTENT( in ) :: kt     ! Main time step counter
50      !!
51      INTEGER               :: ib_bdy, jn ! Loop indeces
52      !!----------------------------------------------------------------------
53      !
54      IF( nn_timing == 1 ) CALL timing_start('trc_bdy')
55      !
56      DO jn = 1, jptra
57         DO ib_bdy=1, nb_bdy
58
59            SELECT CASE( trcdta_bdy(jn,ib_bdy)%cn_obc )
60            CASE('none')
61               CYCLE
62            CASE('frs')
63               CALL bdy_trc_frs( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
64            CASE('specified')
65               CALL bdy_trc_spe( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
66            CASE('neumann')
67               CALL bdy_trc_nmn( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt )
68            CASE('orlanski')
69               CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.false. )
70            CASE('orlanski_npo')
71               CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.true. )
72            CASE DEFAULT
73               CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
74            END SELECT
75
76            ! Boundary points should be updated
77            CALL lbc_bdy_lnk( tra(:,:,:,jn), 'T', 1., ib_bdy )
78
79         ENDDO
80      ENDDO
81      !
82      IF( nn_timing == 1 ) CALL timing_stop('trc_bdy')
83
84   END SUBROUTINE trc_bdy
85
86   SUBROUTINE bdy_trc_frs( jn, idx, dta, kt )
87      !!----------------------------------------------------------------------
88      !!                 ***  SUBROUTINE bdy_trc_frs  ***
89      !!                   
90      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
91      !!
92      !! Reference : Engedahl H., 1995, Tellus, 365-382.
93      !!----------------------------------------------------------------------
94      INTEGER,         INTENT(in) ::   kt
95      INTEGER,         INTENT(in) ::   jn   ! Tracer index
96      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
97      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
98      !!
99      REAL(wp) ::   zwgt           ! boundary weight
100      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
101      INTEGER  ::   ii, ij         ! 2D addresses
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_frs')
105      !
106      igrd = 1                       ! Everything is at T-points here
107      DO ib = 1, idx%nblen(igrd)
108         DO ik = 1, jpkm1
109            ii = idx%nbi(ib,igrd)
110            ij = idx%nbj(ib,igrd)
111            zwgt = idx%nbw(ib,igrd)
112            tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) + zwgt * ( ( dta%trc(ib,ik) * dta%rn_fac)  & 
113                        &  - tra(ii,ij,ik,jn) ) ) * tmask(ii,ij,ik)
114         END DO
115      END DO 
116      !
117      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
118      !
119      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_frs')
120      !
121   END SUBROUTINE bdy_trc_frs
122 
123   SUBROUTINE bdy_trc_spe( jn, idx, dta, kt )
124      !!----------------------------------------------------------------------
125      !!                 ***  SUBROUTINE bdy_trc_frs  ***
126      !!                   
127      !! ** Purpose : Apply a specified value for tracers at open boundaries.
128      !!
129      !!----------------------------------------------------------------------
130      INTEGER,         INTENT(in) ::   kt
131      INTEGER,         INTENT(in) ::   jn   ! Tracer index
132      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
133      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
134      !!
135      REAL(wp) ::   zwgt           ! boundary weight
136      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
137      INTEGER  ::   ii, ij         ! 2D addresses
138      !!----------------------------------------------------------------------
139      !
140      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_spe')
141      !
142      igrd = 1                       ! Everything is at T-points here
143      DO ib = 1, idx%nblenrim(igrd)
144         ii = idx%nbi(ib,igrd)
145         ij = idx%nbj(ib,igrd)
146         DO ik = 1, jpkm1
147            tra(ii,ij,ik,jn) = dta%trc(ib,ik) * dta%rn_fac * tmask(ii,ij,ik)
148         END DO
149      END DO
150      !
151      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
152      !
153      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_spe')
154      !
155   END SUBROUTINE bdy_trc_spe
156
157   SUBROUTINE bdy_trc_nmn( jn, idx, dta, kt )
158      !!----------------------------------------------------------------------
159      !!                 ***  SUBROUTINE bdy_trc_nmn  ***
160      !!                   
161      !! ** Purpose : Duplicate the value for tracers at open boundaries.
162      !!
163      !!----------------------------------------------------------------------
164      INTEGER,         INTENT(in) ::   kt
165      INTEGER,         INTENT(in) ::   jn   ! Tracer index
166      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
167      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
168      !!
169      REAL(wp) ::   zcoef, zcoef1, zcoef2           ! boundary weight
170      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
171      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
172      !!----------------------------------------------------------------------
173      !
174      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_nmn')
175      !
176      igrd = 1                       ! Everything is at T-points here
177      DO ib = 1, idx%nblenrim(igrd)
178         ii = idx%nbi(ib,igrd)
179         ij = idx%nbj(ib,igrd)
180         DO ik = 1, jpkm1
181            ! search the sense of the gradient
182            zcoef1 = bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik) +  bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik)
183            zcoef2 = bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik) +  bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik)
184            IF ( nint(zcoef1+zcoef2) == 0) THEN
185               ! corner
186               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
187               IF (zcoef > .5_wp) THEN ! Only set not isolated points.
188                 tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik) + &
189                   &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik) + &
190                   &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik) + &
191                   &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)
192                 tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / zcoef ) * tmask(ii,ij,ik)
193               ENDIF
194            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN
195               ! oblique corner
196               zcoef = tmask(ii-1,ij,ik)*bdytmask(ii-1,ij  ) + tmask(ii+1,ij,ik)*bdytmask(ii+1,ij  ) + &
197                  &  tmask(ii,ij-1,ik)*bdytmask(ii,ij -1 ) +  tmask(ii,ij+1,ik)*bdytmask(ii,ij+1  )
198               tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik)*bdytmask(ii-1,ij  ) + &
199                  &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik)*bdytmask(ii+1,ij  )  + &
200                  &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik)*bdytmask(ii,ij -1 ) + &
201                  &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)*bdytmask(ii,ij+1  )
202 
203               tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / MAX(1._wp, zcoef) ) * tmask(ii,ij,ik)
204            ELSE
205               ip = nint(bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik) - bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik))
206               jp = nint(bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik) - bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik))
207               tra(ii,ij,ik,jn) = tra(ii+ip,ij+jp,ik,jn) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
208            ENDIF
209         END DO
210      END DO
211      !
212      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
213      !
214      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_nmn')
215      !
216   END SUBROUTINE bdy_trc_nmn
217 
218
219   SUBROUTINE bdy_trc_orlanski( jn, idx, dta, ll_npo )
220      !!----------------------------------------------------------------------
221      !!                 ***  SUBROUTINE bdy_trc_orlanski  ***
222      !!             
223      !!              - Apply Orlanski radiation to tracers of TOP component.
224      !!              - Wrapper routine for bdy_orlanski_3d
225      !!
226      !!
227      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
228      !!----------------------------------------------------------------------
229      INTEGER,                      INTENT(in) ::   jn      ! Tracer index
230      TYPE(OBC_INDEX),              INTENT(in) ::   idx     ! OBC indices
231      TYPE(OBC_DATA),               INTENT(in) ::   dta     ! OBC external data
232      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
233
234      INTEGER  ::   igrd                                    ! grid index
235      !!----------------------------------------------------------------------
236
237      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_orlanski')
238      !
239      igrd = 1      ! Orlanski bc on tracers;
240      !           
241      CALL bdy_orlanski_3d( idx, igrd, trb(:,:,:,jn), tra(:,:,:,jn), (dta%trc * dta%rn_fac), ll_npo )
242      !
243      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_orlanski')
244      !
245
246   END SUBROUTINE bdy_trc_orlanski
247
248   SUBROUTINE trc_bdy_dmp( kt )
249      !!----------------------------------------------------------------------
250      !!                 ***  SUBROUTINE trc_bdy_dmp  ***
251      !!                   
252      !! ** Purpose : Apply damping for tracers at open boundaries.
253      !!             It currently applies the damping to all tracers!!!
254      !!
255      !!----------------------------------------------------------------------
256      INTEGER,         INTENT(in) ::   kt
257      !!
258      INTEGER  ::   jn             ! Tracer index
259      REAL(wp) ::   zwgt           ! boundary weight
260      REAL(wp) ::   zta, zsa, ztime
261      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
262      INTEGER  ::   ii, ij         ! 2D addresses
263      INTEGER  ::   ib_bdy         ! Loop index
264      !!----------------------------------------------------------------------
265      !
266      IF( nn_timing == 1 ) CALL timing_start('trc_bdy_dmp')
267      !
268      DO jn = 1, jptra
269         DO ib_bdy=1, nb_bdy
270            IF ( trcdta_bdy(jn, ib_bdy)%dmp ) THEN
271               igrd = 1                       ! Everything is at T-points here
272               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
273                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
274                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
275                  zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
276                  DO ik = 1, jpkm1
277                     zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - trb(ii,ij,ik,jn) ) * tmask(ii,ij,ik)
278                     tra(ii,ij,ik,jn) = tra(ii,ij,ik,jn) + zta
279                  END DO
280               END DO
281            ENDIF
282         ENDDO
283      ENDDO
284      !
285      IF( nn_timing == 1 ) CALL timing_stop('trc_bdy_dmp')
286      !
287   END SUBROUTINE trc_bdy_dmp
288 
289#else
290   !!----------------------------------------------------------------------
291   !!   Dummy module                   NO Unstruct Open Boundary Conditions
292   !!----------------------------------------------------------------------
293CONTAINS
294   SUBROUTINE trc_bdy(kt)      ! Empty routine
295      WRITE(*,*) 'trc_bdy: You should not have seen this print! error?', kt
296   END SUBROUTINE trc_bdy
297
298   SUBROUTINE trc_bdy_dmp(kt)      ! Empty routine
299      WRITE(*,*) 'trc_bdy_dmp: You should not have seen this print! error?', kt
300   END SUBROUTINE trc_bdy_dmp
301
302#endif
303
304   !!======================================================================
305END MODULE trcbdy
Note: See TracBrowser for help on using the repository browser.