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.
bdylib.F90 in NEMO/trunk/src/OCE/BDY – NEMO

source: NEMO/trunk/src/OCE/BDY/bdylib.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 27.7 KB
Line 
1MODULE bdylib
2   !!======================================================================
3   !!                       ***  MODULE  bdylib  ***
4   !! Unstructured Open Boundary Cond. :  Library module of generic boundary algorithms.
5   !!======================================================================
6   !! History :  3.6  !  2013     (D. Storkey) original code
7   !!            4.0  !  2014     (T. Lovato) Generalize OBC structure
8   !!----------------------------------------------------------------------
9   !!----------------------------------------------------------------------
10   !!   bdy_orlanski_2d
11   !!   bdy_orlanski_3d
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE bdy_oce        ! ocean open boundary conditions
16   USE phycst         ! physical constants
17   USE bdyini
18   !
19   USE in_out_manager !
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21   USE lib_mpp, ONLY: ctl_stop
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   bdy_frs, bdy_spe, bdy_nmn, bdy_orl
27   PUBLIC   bdy_orlanski_2d
28   PUBLIC   bdy_orlanski_3d
29
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
33   !! Software governed by the CeCILL license (see ./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE bdy_frs( idx, phia, dta )
38      !!----------------------------------------------------------------------
39      !!                 ***  SUBROUTINE bdy_frs  ***
40      !!
41      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
42      !!
43      !! Reference : Engedahl H., 1995, Tellus, 365-382.
44      !!----------------------------------------------------------------------
45      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
46      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
47      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
48      !!
49      REAL(wp) ::   zwgt           ! boundary weight
50      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
51      INTEGER  ::   ii, ij         ! 2D addresses
52      !!----------------------------------------------------------------------
53      !
54      igrd = 1                       ! Everything is at T-points here
55      DO ib = 1, idx%nblen(igrd)
56         DO ik = 1, jpkm1
57            ii = idx%nbi(ib,igrd) 
58            ij = idx%nbj(ib,igrd)
59            zwgt = idx%nbw(ib,igrd)
60            phia(ii,ij,ik) = ( phia(ii,ij,ik) + zwgt * (dta(ib,ik) - phia(ii,ij,ik) ) ) * tmask(ii,ij,ik)
61         END DO
62      END DO
63      !
64   END SUBROUTINE bdy_frs
65
66
67   SUBROUTINE bdy_spe( idx, phia, dta )
68      !!----------------------------------------------------------------------
69      !!                 ***  SUBROUTINE bdy_spe  ***
70      !!
71      !! ** Purpose : Apply a specified value for tracers at open boundaries.
72      !!
73      !!----------------------------------------------------------------------
74      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
75      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
76      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
77      !!
78      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
79      INTEGER  ::   ii, ij         ! 2D addresses
80      !!----------------------------------------------------------------------
81      !
82      igrd = 1                       ! Everything is at T-points here
83      DO ib = 1, idx%nblenrim(igrd)
84         ii = idx%nbi(ib,igrd)
85         ij = idx%nbj(ib,igrd)
86         DO ik = 1, jpkm1
87            phia(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
88         END DO
89      END DO
90      !
91   END SUBROUTINE bdy_spe
92
93
94   SUBROUTINE bdy_orl( idx, phib, phia, dta, lrim0, ll_npo )
95      !!----------------------------------------------------------------------
96      !!                 ***  SUBROUTINE bdy_orl  ***
97      !!
98      !! ** Purpose : Apply Orlanski radiation for tracers at open boundaries.
99      !!              This is a wrapper routine for bdy_orlanski_3d below
100      !!
101      !!----------------------------------------------------------------------
102      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
103      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
104      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phib  ! before tracer field
105      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
106      LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated
107      LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version
108      !!
109      INTEGER  ::   igrd                                    ! grid index
110      !!----------------------------------------------------------------------
111      !
112      igrd = 1                       ! Everything is at T-points here
113      !
114      CALL bdy_orlanski_3d( idx, igrd, phib(:,:,:), phia(:,:,:), dta, lrim0, ll_npo )
115      !
116   END SUBROUTINE bdy_orl
117
118
119   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, lrim0, ll_npo )
120      !!----------------------------------------------------------------------
121      !!                 ***  SUBROUTINE bdy_orlanski_2d  ***
122      !!             
123      !!              - Apply Orlanski radiation condition adaptively to 2D fields:
124      !!                  - radiation plus weak nudging at outflow points
125      !!                  - no radiation and strong nudging at inflow points
126      !!
127      !!
128      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
129      !!----------------------------------------------------------------------
130      TYPE(OBC_INDEX),          INTENT(in   ) ::   idx      ! BDY indices
131      INTEGER ,                 INTENT(in   ) ::   igrd     ! grid index
132      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phib     ! model before 2D field
133      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   phia     ! model after 2D field (to be updated)
134      REAL(wp), DIMENSION(:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
135      LOGICAL, OPTIONAL,        INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
136      LOGICAL ,                 INTENT(in   ) ::   ll_npo   ! switch for NPO version
137      !
138      INTEGER  ::   jb                                     ! dummy loop indices
139      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
140      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
141      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
142      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
143      INTEGER  ::   flagu, flagv                           ! short cuts
144      INTEGER  ::   ibeg, iend                             ! length of rim to be treated (rim 0 or rim 1 or both)
145      REAL(wp) ::   zmask_x, zmask_y1, zmask_y2
146      REAL(wp) ::   zex1, zex2, zey, zey1, zey2
147      REAL(wp) ::   zdt, zdx, zdy, znor2, zrx, zry         ! intermediate calculations
148      REAL(wp) ::   zout, zwgt, zdy_centred
149      REAL(wp) ::   zdy_1, zdy_2, zsign_ups
150      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
151      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask      ! land/sea mask for field
152      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_xdif ! land/sea mask for x-derivatives
153      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_ydif ! land/sea mask for y-derivatives
154      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives
155      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
156      !!----------------------------------------------------------------------
157      !
158      ! ----------------------------------!
159      ! Orlanski boundary conditions     :!
160      ! ----------------------------------!
161     
162      SELECT CASE(igrd)
163         CASE(1)
164            zmask      => tmask(:,:,1)
165            zmask_xdif => umask(:,:,1)
166            zmask_ydif => vmask(:,:,1)
167            pe_xdif    => e1u(:,:)
168            pe_ydif    => e2v(:,:)
169            ii_offset = 0
170            ij_offset = 0
171         CASE(2)
172            zmask      => umask(:,:,1)
173            zmask_xdif => tmask(:,:,1)
174            zmask_ydif => fmask(:,:,1)
175            pe_xdif    => e1t(:,:)
176            pe_ydif    => e2f(:,:)
177            ii_offset = 1
178            ij_offset = 0
179         CASE(3)
180            zmask      => vmask(:,:,1)
181            zmask_xdif => fmask(:,:,1)
182            zmask_ydif => tmask(:,:,1)
183            pe_xdif    => e1f(:,:)
184            pe_ydif    => e2t(:,:)
185            ii_offset = 0
186            ij_offset = 1
187         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
188      END SELECT
189      !
190      IF( PRESENT(lrim0) ) THEN
191         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
192         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
193         END IF
194      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
195      END IF
196      !
197      DO jb = ibeg, iend
198         ii  = idx%nbi(jb,igrd)
199         ij  = idx%nbj(jb,igrd) 
200         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
201         flagu = int( idx%flagu(jb,igrd) )
202         flagv = int( idx%flagv(jb,igrd) )
203         !
204         ! Calculate positions of b-1 and b-2 points for this rim point
205         ! also (b-1,j-1) and (b-1,j+1) points
206         iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
207         ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
208          !
209         iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
210         ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
211         !
212         iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
213         ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
214         !
215         ! Calculate scale factors for calculation of spatial derivatives.       
216         zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
217        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1   +ij_offset) ) 
218         zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   &
219        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2   +ij_offset) ) 
220         zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   & 
221        &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
222         zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
223        &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1   +ij_offset) ) 
224         ! make sure scale factors are nonzero
225         if( zey1 .lt. rsmall ) zey1 = zey2
226         if( zey2 .lt. rsmall ) zey2 = zey1
227         zex1 = max(zex1,rsmall); zex2 = max(zex2,rsmall)
228         zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall); 
229         !
230         ! Calculate masks for calculation of spatial derivatives.
231         zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2               )   &
232        &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset) ) 
233         zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1            )   & 
234        &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
235         zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1               )   &
236        &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset) ) 
237
238         ! Calculation of terms required for both versions of the scheme.
239         ! Mask derivatives to ensure correct land boundary conditions for each variable.
240         ! Centred derivative is calculated as average of "left" and "right" derivatives for
241         ! this reason.
242         ! Note no rdt factor in expression for zdt because it cancels in the expressions for
243         ! zrx and zry.
244         zdt   =     phia(iibm1   ,ijbm1   ) - phib(iibm1   ,ijbm1   )
245         zdx   = ( ( phia(iibm1   ,ijbm1   ) - phia(iibm2   ,ijbm2   ) ) / zex2 ) * zmask_x 
246         zdy_1 = ( ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1   
247         zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) ) / zey2 ) * zmask_y2 
248         zdy_centred = 0.5 * ( zdy_1 + zdy_2 )
249!!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
250         ! upstream differencing for tangential derivatives
251         zsign_ups = sign( 1., zdt * zdy_centred )
252         zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
253         zdy = zsign_ups * zdy_1 + (1. - zsign_ups) * zdy_2
254         znor2 = zdx * zdx + zdy * zdy
255         znor2 = max(znor2,zepsilon)
256         !
257         zrx = zdt * zdx / ( zex1 * znor2 ) 
258!!$         zrx = min(zrx,2.0_wp)
259         zout = sign( 1., zrx )
260         zout = 0.5*( zout + abs(zout) )
261         zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
262         ! only apply radiation on outflow points
263         if( ll_npo ) then     !! NPO version !!
264            phia(ii,ij) = (1.-zout) * ( phib(ii,ij) + zwgt * ( phi_ext(jb) - phib(ii,ij) ) )        &
265           &            + zout      * ( phib(ii,ij) + zrx*phia(iibm1,ijbm1)                         &
266           &                            + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 
267         else                  !! full oblique radiation !!
268            zsign_ups = sign( 1., zdt * zdy )
269            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
270            zey = zsign_ups * zey1 + (1.-zsign_ups) * zey2 
271            zry = zdt * zdy / ( zey * znor2 ) 
272            phia(ii,ij) = (1.-zout) * ( phib(ii,ij) + zwgt * ( phi_ext(jb) - phib(ii,ij) ) )        &
273           &            + zout      * ( phib(ii,ij) + zrx*phia(iibm1,ijbm1)                         &
274           &                    - zsign_ups      * zry * ( phib(ii   ,ij   ) - phib(iijm1,ijjm1 ) ) &
275           &                    - (1.-zsign_ups) * zry * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) &
276           &                    + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 
277         end if
278         phia(ii,ij) = phia(ii,ij) * zmask(ii,ij)
279      END DO
280      !
281   END SUBROUTINE bdy_orlanski_2d
282
283
284   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, lrim0, ll_npo )
285      !!----------------------------------------------------------------------
286      !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
287      !!             
288      !!              - Apply Orlanski radiation condition adaptively to 3D fields:
289      !!                  - radiation plus weak nudging at outflow points
290      !!                  - no radiation and strong nudging at inflow points
291      !!
292      !!
293      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
294      !!----------------------------------------------------------------------
295      TYPE(OBC_INDEX),            INTENT(in   ) ::   idx      ! BDY indices
296      INTEGER ,                   INTENT(in   ) ::   igrd     ! grid index
297      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phib     ! model before 3D field
298      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phia     ! model after 3D field (to be updated)
299      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
300      LOGICAL, OPTIONAL,          INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
301      LOGICAL ,                   INTENT(in   ) ::   ll_npo   ! switch for NPO version
302      !
303      INTEGER  ::   jb, jk                                 ! dummy loop indices
304      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
305      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
306      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
307      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
308      INTEGER  ::   flagu, flagv                           ! short cuts
309      INTEGER  ::   ibeg, iend                             ! length of rim to be treated (rim 0 or rim 1 or both)
310      REAL(wp) ::   zmask_x, zmask_y1, zmask_y2
311      REAL(wp) ::   zex1, zex2, zey, zey1, zey2
312      REAL(wp) ::   zdt, zdx, zdy, znor2, zrx, zry         ! intermediate calculations
313      REAL(wp) ::   zout, zwgt, zdy_centred
314      REAL(wp) ::   zdy_1, zdy_2,  zsign_ups
315      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
316      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask      ! land/sea mask for field
317      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_xdif ! land/sea mask for x-derivatives
318      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_ydif ! land/sea mask for y-derivatives
319      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives
320      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
321      !!----------------------------------------------------------------------
322      !
323      ! ----------------------------------!
324      ! Orlanski boundary conditions     :!
325      ! ----------------------------------!
326      !
327      SELECT CASE(igrd)
328         CASE(1)
329            zmask      => tmask(:,:,:)
330            zmask_xdif => umask(:,:,:)
331            zmask_ydif => vmask(:,:,:)
332            pe_xdif    => e1u(:,:)
333            pe_ydif    => e2v(:,:)
334            ii_offset = 0
335            ij_offset = 0
336         CASE(2)
337            zmask      => umask(:,:,:)
338            zmask_xdif => tmask(:,:,:)
339            zmask_ydif => fmask(:,:,:)
340            pe_xdif    => e1t(:,:)
341            pe_ydif    => e2f(:,:)
342            ii_offset = 1
343            ij_offset = 0
344         CASE(3)
345            zmask      => vmask(:,:,:)
346            zmask_xdif => fmask(:,:,:)
347            zmask_ydif => tmask(:,:,:)
348            pe_xdif    => e1f(:,:)
349            pe_ydif    => e2t(:,:)
350            ii_offset = 0
351            ij_offset = 1
352         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
353      END SELECT
354      !
355      IF( PRESENT(lrim0) ) THEN
356         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
357         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
358         END IF
359      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
360      END IF
361      !
362      DO jk = 1, jpk
363         !           
364         DO jb = ibeg, iend
365            ii  = idx%nbi(jb,igrd)
366            ij  = idx%nbj(jb,igrd) 
367            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
368            flagu = int( idx%flagu(jb,igrd) )
369            flagv = int( idx%flagv(jb,igrd) )
370            !
371            ! calculate positions of b-1 and b-2 points for this rim point
372            ! also (b-1,j-1) and (b-1,j+1) points
373            iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
374            ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
375            !
376            iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
377            ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
378            !
379            iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
380            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
381            !
382            ! Calculate scale factors for calculation of spatial derivatives.       
383            zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
384           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) ) 
385            zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   &
386           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2+ij_offset   ) ) 
387            zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   & 
388           &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
389            zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
390           &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) ) 
391            ! make sure scale factors are nonzero
392            if( zey1 .lt. rsmall ) zey1 = zey2
393            if( zey2 .lt. rsmall ) zey2 = zey1
394            zex1 = max(zex1,rsmall); zex2 = max(zex2,rsmall); 
395            zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall); 
396            !
397            ! Calculate masks for calculation of spatial derivatives.       
398            zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2             ,jk)   &
399           &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset,jk) ) 
400            zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   & 
401           &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) ) 
402            zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1             ,jk)   &
403           &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset,jk) ) 
404            !
405            ! Calculate normal (zrx) and tangential (zry) components of radiation velocities.
406            ! Mask derivatives to ensure correct land boundary conditions for each variable.
407            ! Centred derivative is calculated as average of "left" and "right" derivatives for
408            ! this reason.
409            zdt   =     phia(iibm1   ,ijbm1   ,jk) - phib(iibm1   ,ijbm1   ,jk)
410            zdx   = ( ( phia(iibm1   ,ijbm1   ,jk) - phia(iibm2   ,ijbm2   ,jk) ) / zex2 ) * zmask_x                 
411            zdy_1 = ( ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1 
412            zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) ) / zey2 ) * zmask_y2     
413            zdy_centred = 0.5 * ( zdy_1 + zdy_2 )
414!!$            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk)
415            ! upstream differencing for tangential derivatives
416            zsign_ups = sign( 1., zdt * zdy_centred )
417            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
418            zdy = zsign_ups * zdy_1 + (1. - zsign_ups) * zdy_2
419            znor2 = zdx * zdx + zdy * zdy
420            znor2 = max(znor2,zepsilon)
421            !
422            ! update boundary value:
423            zrx = zdt * zdx / ( zex1 * znor2 )
424!!$            zrx = min(zrx,2.0_wp)
425            zout = sign( 1., zrx )
426            zout = 0.5*( zout + abs(zout) )
427            zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
428            ! only apply radiation on outflow points
429            if( ll_npo ) then     !! NPO version !!
430               phia(ii,ij,jk) = (1.-zout) * ( phib(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) &
431              &               + zout      * ( phib(ii,ij,jk) + zrx*phia(iibm1,ijbm1,jk)                     &
432              &                            + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 
433            else                  !! full oblique radiation !!
434               zsign_ups = sign( 1., zdt * zdy )
435               zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
436               zey = zsign_ups * zey1 + (1.-zsign_ups) * zey2 
437               zry = zdt * zdy / ( zey * znor2 ) 
438               phia(ii,ij,jk) = (1.-zout) * ( phib(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) )    &
439              &               + zout      * ( phib(ii,ij,jk) + zrx*phia(iibm1,ijbm1,jk)                        &
440              &                       - zsign_ups      * zry * ( phib(ii   ,ij   ,jk) - phib(iijm1,ijjm1,jk) ) &
441              &                       - (1.-zsign_ups) * zry * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk) ) &
442              &                       + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 
443            end if
444            phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk)
445         END DO
446         !
447      END DO
448      !
449   END SUBROUTINE bdy_orlanski_3d
450
451   SUBROUTINE bdy_nmn( idx, igrd, phia, lrim0 )
452      !!----------------------------------------------------------------------
453      !!                 ***  SUBROUTINE bdy_nmn  ***
454      !!                   
455      !! ** Purpose : Duplicate the value at open boundaries, zero gradient.
456      !!
457      !!
458      !! ** Method  : - take the average of free ocean neighbours
459      !!
460      !!      ___   !   |_____|   !   ___|    !   __|x o   !   |_   _|     ! |     
461      !!   __|x     !      x      !     x o   !      o     !     |_|       ! |x o   
462      !!      o     !      o      !     o     !            !    o x o      ! |x_x_
463      !!                                                   !      o     
464      !!----------------------------------------------------------------------
465      INTEGER,                    INTENT(in   )  ::   igrd     ! grid index
466      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated), must be masked
467      TYPE(OBC_INDEX),            INTENT(in   )  ::   idx      ! OBC indices
468      LOGICAL, OPTIONAL,          INTENT(in   )  ::   lrim0    ! indicate if rim 0 is treated
469      !!
470      REAL(wp) ::   zweight
471      REAL(wp), POINTER, DIMENSION(:,:,:)      :: zmask         ! land/sea mask for field
472      INTEGER  ::   ib, ik   ! dummy loop indices
473      INTEGER  ::   ii, ij   ! 2D addresses
474      INTEGER  ::   ipkm1    ! size of phia third dimension minus 1
475      INTEGER  ::   ibeg, iend                          ! length of rim to be treated (rim 0 or rim 1 or both)
476      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, itreat
477      !!----------------------------------------------------------------------
478      !
479      ipkm1 = MAX( SIZE(phia,3) - 1, 1 ) 
480      !
481      SELECT CASE(igrd)
482         CASE(1)   ;   zmask => tmask(:,:,:)
483         CASE(2)   ;   zmask => umask(:,:,:)
484         CASE(3)   ;   zmask => vmask(:,:,:)
485         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )
486      END SELECT
487      !
488      IF( PRESENT(lrim0) ) THEN
489         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
490         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
491         END IF
492      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
493      END IF
494      !
495      DO ib = ibeg, iend
496         ii = idx%nbi(ib,igrd)
497         ij = idx%nbj(ib,igrd)
498         itreat = idx%ntreat(ib,igrd)
499         CALL find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )   ! find free ocean neighbours
500         SELECT CASE( itreat )
501         CASE( 1:8 )
502            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
503            DO ik = 1, ipkm1
504               IF( zmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik) 
505            END DO
506         CASE( 9:12 )
507            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
508            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
509            DO ik = 1, ipkm1
510               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik)
511               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight
512            END DO
513         CASE( 13:16 )
514            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
515            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
516            IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE
517            DO ik = 1, ipkm1
518               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) + zmask(ii3,ij3,ik)
519               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight
520            END DO
521         END SELECT
522      END DO
523      !
524   END SUBROUTINE bdy_nmn
525
526   !!======================================================================
527END MODULE bdylib
Note: See TracBrowser for help on using the repository browser.