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/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90 @ 11979

Last change on this file since 11979 was 11978, checked in by dford, 4 years ago

Update BDY scheme to remove dead zone.

File size: 14.0 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) * rdttrc(ik) / 86400.d0  ! damping with a timescale of day
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! tra is overwritten at the boundary so damping doesn't work here - need neumann in addition
117! to duplicate the internal value at the boundary
118       CALL bdy_trc_nmn( jn,  idx, dta, kt )
119
120
121      !
122      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
123      !
124      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_frs')
125      !
126   END SUBROUTINE bdy_trc_frs
127 
128   SUBROUTINE bdy_trc_spe( jn, idx, dta, kt )
129      !!----------------------------------------------------------------------
130      !!                 ***  SUBROUTINE bdy_trc_frs  ***
131      !!                   
132      !! ** Purpose : Apply a specified value for tracers at open boundaries.
133      !!
134      !!----------------------------------------------------------------------
135      INTEGER,         INTENT(in) ::   kt
136      INTEGER,         INTENT(in) ::   jn   ! Tracer index
137      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
138      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
139      !!
140      REAL(wp) ::   zwgt           ! boundary weight
141      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
142      INTEGER  ::   ii, ij         ! 2D addresses
143      !!----------------------------------------------------------------------
144      !
145      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_spe')
146      !
147      igrd = 1                       ! Everything is at T-points here
148      DO ib = 1, idx%nblenrim(igrd)
149         ii = idx%nbi(ib,igrd)
150         ij = idx%nbj(ib,igrd)
151         DO ik = 1, jpkm1
152            tra(ii,ij,ik,jn) = dta%trc(ib,ik) * dta%rn_fac * tmask(ii,ij,ik)
153         END DO
154      END DO
155      !
156      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
157      !
158      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_spe')
159      !
160   END SUBROUTINE bdy_trc_spe
161
162   SUBROUTINE bdy_trc_nmn( jn, idx, dta, kt )
163      !!----------------------------------------------------------------------
164      !!                 ***  SUBROUTINE bdy_trc_nmn  ***
165      !!                   
166      !! ** Purpose : Duplicate the value for tracers at open boundaries.
167      !!
168      !!----------------------------------------------------------------------
169      INTEGER,         INTENT(in) ::   kt
170      INTEGER,         INTENT(in) ::   jn   ! Tracer index
171      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
172      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
173      !!
174      REAL(wp) ::   zcoef, zcoef1, zcoef2           ! boundary weight
175      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
176      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
177      !!----------------------------------------------------------------------
178      !
179      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_nmn')
180      !
181      igrd = 1                       ! Everything is at T-points here
182      DO ib = 1, idx%nblenrim(igrd)
183         ii = idx%nbi(ib,igrd)
184         ij = idx%nbj(ib,igrd)
185         DO ik = 1, jpkm1
186            ! search the sense of the gradient
187            zcoef1 = bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik) +  bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik)
188            zcoef2 = bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik) +  bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik)
189            IF ( nint(zcoef1+zcoef2) == 0) THEN
190               ! corner
191               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
192               IF (zcoef > .5_wp) THEN ! Only set not isolated points.
193                 tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik) + &
194                   &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik) + &
195                   &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik) + &
196                   &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)
197                 tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / zcoef ) * tmask(ii,ij,ik)
198               ENDIF
199            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN
200               ! oblique corner
201               zcoef = tmask(ii-1,ij,ik)*bdytmask(ii-1,ij  ) + tmask(ii+1,ij,ik)*bdytmask(ii+1,ij  ) + &
202                  &  tmask(ii,ij-1,ik)*bdytmask(ii,ij -1 ) +  tmask(ii,ij+1,ik)*bdytmask(ii,ij+1  )
203               tra(ii,ij,ik,jn) = tra(ii-1,ij  ,ik,jn) * tmask(ii-1,ij  ,ik)*bdytmask(ii-1,ij  ) + &
204                  &              tra(ii+1,ij  ,ik,jn) * tmask(ii+1,ij  ,ik)*bdytmask(ii+1,ij  )  + &
205                  &              tra(ii  ,ij-1,ik,jn) * tmask(ii  ,ij-1,ik)*bdytmask(ii,ij -1 ) + &
206                  &              tra(ii  ,ij+1,ik,jn) * tmask(ii  ,ij+1,ik)*bdytmask(ii,ij+1  )
207 
208               tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / MAX(1._wp, zcoef) ) * tmask(ii,ij,ik)
209            ELSE
210               ip = nint(bdytmask(ii+1,ij  )*tmask(ii+1,ij,ik) - bdytmask(ii-1,ij  )*tmask(ii-1,ij,ik))
211               jp = nint(bdytmask(ii  ,ij+1)*tmask(ii,ij+1,ik) - bdytmask(ii  ,ij-1)*tmask(ii,ij-1,ik))
212               tra(ii,ij,ik,jn) = tra(ii+ip,ij+jp,ik,jn) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
213            ENDIF
214         END DO
215      END DO
216      !
217      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
218      !
219      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_nmn')
220      !
221   END SUBROUTINE bdy_trc_nmn
222 
223
224   SUBROUTINE bdy_trc_orlanski( jn, idx, dta, ll_npo )
225      !!----------------------------------------------------------------------
226      !!                 ***  SUBROUTINE bdy_trc_orlanski  ***
227      !!             
228      !!              - Apply Orlanski radiation to tracers of TOP component.
229      !!              - Wrapper routine for bdy_orlanski_3d
230      !!
231      !!
232      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
233      !!----------------------------------------------------------------------
234      INTEGER,                      INTENT(in) ::   jn      ! Tracer index
235      TYPE(OBC_INDEX),              INTENT(in) ::   idx     ! OBC indices
236      TYPE(OBC_DATA),               INTENT(in) ::   dta     ! OBC external data
237      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
238
239      INTEGER  ::   igrd                                    ! grid index
240      !!----------------------------------------------------------------------
241
242      IF( nn_timing == 1 ) CALL timing_start('bdy_trc_orlanski')
243      !
244      igrd = 1      ! Orlanski bc on tracers;
245      !           
246      CALL bdy_orlanski_3d( idx, igrd, trb(:,:,:,jn), tra(:,:,:,jn), (dta%trc * dta%rn_fac), ll_npo )
247      !
248      IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_orlanski')
249      !
250
251   END SUBROUTINE bdy_trc_orlanski
252
253   SUBROUTINE trc_bdy_dmp( kt )
254      !!----------------------------------------------------------------------
255      !!                 ***  SUBROUTINE trc_bdy_dmp  ***
256      !!                   
257      !! ** Purpose : Apply damping for tracers at open boundaries.
258      !!             It currently applies the damping to all tracers!!!
259      !!
260      !!----------------------------------------------------------------------
261      INTEGER,         INTENT(in) ::   kt
262      !!
263      INTEGER  ::   jn             ! Tracer index
264      REAL(wp) ::   zwgt           ! boundary weight
265      REAL(wp) ::   zta, zsa, ztime
266      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
267      INTEGER  ::   ii, ij         ! 2D addresses
268      INTEGER  ::   ib_bdy         ! Loop index
269      !!----------------------------------------------------------------------
270      !
271      IF( nn_timing == 1 ) CALL timing_start('trc_bdy_dmp')
272      !
273      DO jn = 1, jptra
274         DO ib_bdy=1, nb_bdy
275            IF ( trcdta_bdy(jn, ib_bdy)%dmp ) THEN
276               igrd = 1                       ! Everything is at T-points here
277               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
278                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
279                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
280                  zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
281                  DO ik = 1, jpkm1
282                     zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - trb(ii,ij,ik,jn) ) * tmask(ii,ij,ik)
283                     tra(ii,ij,ik,jn) = tra(ii,ij,ik,jn) + zta
284                  END DO
285               END DO
286            ENDIF
287         ENDDO
288      ENDDO
289      !
290      IF( nn_timing == 1 ) CALL timing_stop('trc_bdy_dmp')
291      !
292   END SUBROUTINE trc_bdy_dmp
293 
294#else
295   !!----------------------------------------------------------------------
296   !!   Dummy module                   NO Unstruct Open Boundary Conditions
297   !!----------------------------------------------------------------------
298CONTAINS
299   SUBROUTINE trc_bdy(kt)      ! Empty routine
300      WRITE(*,*) 'trc_bdy: You should not have seen this print! error?', kt
301   END SUBROUTINE trc_bdy
302
303   SUBROUTINE trc_bdy_dmp(kt)      ! Empty routine
304      WRITE(*,*) 'trc_bdy_dmp: You should not have seen this print! error?', kt
305   END SUBROUTINE trc_bdy_dmp
306
307#endif
308
309   !!======================================================================
310END MODULE trcbdy
Note: See TracBrowser for help on using the repository browser.