MODULE trcbdy !!====================================================================== !! *** MODULE bdytrc *** !! Ocean tracers: Apply boundary conditions for tracers in TOP component !!====================================================================== !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code !! 3.0 ! 2008-04 (NEMO team) add in the reference version !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Optimization of BDY communications !! 3.6 ! 2015 (T. Lovato) Adapt BDY for tracers in TOP component !!---------------------------------------------------------------------- #if defined key_bdy && key_top !!---------------------------------------------------------------------- !! 'key_bdy' Unstructured Open Boundary Conditions !!---------------------------------------------------------------------- !! trc_bdy : Apply open boundary conditions to T and S !! trc_bdy_frs : Apply Flow Relaxation Scheme !!---------------------------------------------------------------------- USE timing ! Timing USE oce_trc ! ocean dynamics and tracers variables USE par_trc USE trc ! ocean space and time domain variables USE bdylib ! for orlanski library routines USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE in_out_manager ! I/O manager USE bdy_oce, only: idx_bdy, OBC_INDEX, BDYTMASK, lk_bdy ! ocean open boundary conditions IMPLICIT NONE PRIVATE PUBLIC trc_bdy ! routine called in trcnxt.F90 PUBLIC trc_bdy_dmp ! routine called in trcstp.F90 !!---------------------------------------------------------------------- !! NEMO/OPA 3.6 , NEMO Consortium (2015) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trc_bdy( kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE trc_bdy *** !! !! ** Purpose : - Apply open boundary conditions for tracers in TOP component !! and scale the tracer data !! !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! Main time step counter !! INTEGER :: ib_bdy, jn ! Loop indeces !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('trc_bdy') ! DO jn = 1, jptra DO ib_bdy=1, nb_bdy SELECT CASE( trcdta_bdy(jn,ib_bdy)%cn_obc ) CASE('none') CYCLE CASE('frs') CALL bdy_trc_frs( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt ) CASE('specified') CALL bdy_trc_spe( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt ) CASE('neumann') CALL bdy_trc_nmn( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), kt ) CASE('orlanski') CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.false. ) CASE('orlanski_npo') CALL bdy_trc_orlanski( jn, idx_bdy(ib_bdy), trcdta_bdy(jn,ib_bdy), ll_npo=.true. ) CASE DEFAULT CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) END SELECT ! Boundary points should be updated CALL lbc_bdy_lnk( tra(:,:,:,jn), 'T', 1., ib_bdy ) ENDDO ENDDO ! IF( nn_timing == 1 ) CALL timing_stop('trc_bdy') END SUBROUTINE trc_bdy SUBROUTINE bdy_trc_frs( jn, idx, dta, kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE bdy_trc_frs *** !! !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries. !! !! Reference : Engedahl H., 1995, Tellus, 365-382. !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt INTEGER, INTENT(in) :: jn ! Tracer index TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data !! REAL(wp) :: zwgt ! boundary weight INTEGER :: ib, ik, igrd ! dummy loop indices INTEGER :: ii, ij ! 2D addresses !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('bdy_trc_frs') ! igrd = 1 ! Everything is at T-points here DO ib = 1, idx%nblen(igrd) DO ik = 1, jpkm1 ii = idx%nbi(ib,igrd) ij = idx%nbj(ib,igrd) zwgt = idx%nbw(ib,igrd) * rdttrc(ik) / 86400.d0 ! damping with a timescale of day tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) + zwgt * ( ( dta%trc(ib,ik) * dta%rn_fac) & & - tra(ii,ij,ik,jn) ) ) * tmask(ii,ij,ik) END DO END DO ! tra is overwritten at the boundary so damping doesn't work here - need neumann in addition ! to duplicate the internal value at the boundary CALL bdy_trc_nmn( jn, idx, dta, kt ) ! IF( kt .eq. nit000 ) CLOSE( unit = 102 ) ! IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_frs') ! END SUBROUTINE bdy_trc_frs SUBROUTINE bdy_trc_spe( jn, idx, dta, kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE bdy_trc_frs *** !! !! ** Purpose : Apply a specified value for tracers at open boundaries. !! !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt INTEGER, INTENT(in) :: jn ! Tracer index TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data !! REAL(wp) :: zwgt ! boundary weight INTEGER :: ib, ik, igrd ! dummy loop indices INTEGER :: ii, ij ! 2D addresses !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('bdy_trc_spe') ! igrd = 1 ! Everything is at T-points here DO ib = 1, idx%nblenrim(igrd) ii = idx%nbi(ib,igrd) ij = idx%nbj(ib,igrd) DO ik = 1, jpkm1 tra(ii,ij,ik,jn) = dta%trc(ib,ik) * dta%rn_fac * tmask(ii,ij,ik) END DO END DO ! IF( kt .eq. nit000 ) CLOSE( unit = 102 ) ! IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_spe') ! END SUBROUTINE bdy_trc_spe SUBROUTINE bdy_trc_nmn( jn, idx, dta, kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE bdy_trc_nmn *** !! !! ** Purpose : Duplicate the value for tracers at open boundaries. !! !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt INTEGER, INTENT(in) :: jn ! Tracer index TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data !! REAL(wp) :: zcoef, zcoef1, zcoef2 ! boundary weight INTEGER :: ib, ik, igrd ! dummy loop indices INTEGER :: ii, ij, ip, jp ! 2D addresses !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('bdy_trc_nmn') ! igrd = 1 ! Everything is at T-points here DO ib = 1, idx%nblenrim(igrd) ii = idx%nbi(ib,igrd) ij = idx%nbj(ib,igrd) DO ik = 1, jpkm1 ! search the sense of the gradient zcoef1 = bdytmask(ii-1,ij )*tmask(ii-1,ij,ik) + bdytmask(ii+1,ij )*tmask(ii+1,ij,ik) zcoef2 = bdytmask(ii ,ij-1)*tmask(ii,ij-1,ik) + bdytmask(ii ,ij+1)*tmask(ii,ij+1,ik) IF ( nint(zcoef1+zcoef2) == 0) THEN ! corner zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) + tmask(ii,ij-1,ik) + tmask(ii,ij+1,ik) IF (zcoef > .5_wp) THEN ! Only set not isolated points. tra(ii,ij,ik,jn) = tra(ii-1,ij ,ik,jn) * tmask(ii-1,ij ,ik) + & & tra(ii+1,ij ,ik,jn) * tmask(ii+1,ij ,ik) + & & tra(ii ,ij-1,ik,jn) * tmask(ii ,ij-1,ik) + & & tra(ii ,ij+1,ik,jn) * tmask(ii ,ij+1,ik) tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / zcoef ) * tmask(ii,ij,ik) ENDIF ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN ! oblique corner zcoef = tmask(ii-1,ij,ik)*bdytmask(ii-1,ij ) + tmask(ii+1,ij,ik)*bdytmask(ii+1,ij ) + & & tmask(ii,ij-1,ik)*bdytmask(ii,ij -1 ) + tmask(ii,ij+1,ik)*bdytmask(ii,ij+1 ) tra(ii,ij,ik,jn) = tra(ii-1,ij ,ik,jn) * tmask(ii-1,ij ,ik)*bdytmask(ii-1,ij ) + & & tra(ii+1,ij ,ik,jn) * tmask(ii+1,ij ,ik)*bdytmask(ii+1,ij ) + & & tra(ii ,ij-1,ik,jn) * tmask(ii ,ij-1,ik)*bdytmask(ii,ij -1 ) + & & tra(ii ,ij+1,ik,jn) * tmask(ii ,ij+1,ik)*bdytmask(ii,ij+1 ) tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) / MAX(1._wp, zcoef) ) * tmask(ii,ij,ik) ELSE ip = nint(bdytmask(ii+1,ij )*tmask(ii+1,ij,ik) - bdytmask(ii-1,ij )*tmask(ii-1,ij,ik)) jp = nint(bdytmask(ii ,ij+1)*tmask(ii,ij+1,ik) - bdytmask(ii ,ij-1)*tmask(ii,ij-1,ik)) tra(ii,ij,ik,jn) = tra(ii+ip,ij+jp,ik,jn) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik) ENDIF END DO END DO ! IF( kt .eq. nit000 ) CLOSE( unit = 102 ) ! IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_nmn') ! END SUBROUTINE bdy_trc_nmn SUBROUTINE bdy_trc_orlanski( jn, idx, dta, ll_npo ) !!---------------------------------------------------------------------- !! *** SUBROUTINE bdy_trc_orlanski *** !! !! - Apply Orlanski radiation to tracers of TOP component. !! - Wrapper routine for bdy_orlanski_3d !! !! !! References: Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001) !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: jn ! Tracer index TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data LOGICAL, INTENT(in) :: ll_npo ! switch for NPO version INTEGER :: igrd ! grid index !!---------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('bdy_trc_orlanski') ! igrd = 1 ! Orlanski bc on tracers; ! CALL bdy_orlanski_3d( idx, igrd, trb(:,:,:,jn), tra(:,:,:,jn), (dta%trc * dta%rn_fac), ll_npo ) ! IF( nn_timing == 1 ) CALL timing_stop('bdy_trc_orlanski') ! END SUBROUTINE bdy_trc_orlanski SUBROUTINE trc_bdy_dmp( kt ) !!---------------------------------------------------------------------- !! *** SUBROUTINE trc_bdy_dmp *** !! !! ** Purpose : Apply damping for tracers at open boundaries. !! It currently applies the damping to all tracers!!! !! !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt !! INTEGER :: jn ! Tracer index REAL(wp) :: zwgt ! boundary weight REAL(wp) :: zta, zsa, ztime INTEGER :: ib, ik, igrd ! dummy loop indices INTEGER :: ii, ij ! 2D addresses INTEGER :: ib_bdy ! Loop index !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('trc_bdy_dmp') ! DO jn = 1, jptra DO ib_bdy=1, nb_bdy IF ( trcdta_bdy(jn, ib_bdy)%dmp ) THEN igrd = 1 ! Everything is at T-points here DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ii = idx_bdy(ib_bdy)%nbi(ib,igrd) ij = idx_bdy(ib_bdy)%nbj(ib,igrd) zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd) DO ik = 1, jpkm1 zta = zwgt * ( trcdta_bdy(jn, ib_bdy)%trc(ib,ik) - trb(ii,ij,ik,jn) ) * tmask(ii,ij,ik) tra(ii,ij,ik,jn) = tra(ii,ij,ik,jn) + zta END DO END DO ENDIF ENDDO ENDDO ! IF( nn_timing == 1 ) CALL timing_stop('trc_bdy_dmp') ! END SUBROUTINE trc_bdy_dmp #else !!---------------------------------------------------------------------- !! Dummy module NO Unstruct Open Boundary Conditions !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trc_bdy(kt) ! Empty routine WRITE(*,*) 'trc_bdy: You should not have seen this print! error?', kt END SUBROUTINE trc_bdy SUBROUTINE trc_bdy_dmp(kt) ! Empty routine WRITE(*,*) 'trc_bdy_dmp: You should not have seen this print! error?', kt END SUBROUTINE trc_bdy_dmp #endif !!====================================================================== END MODULE trcbdy