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

Last change on this file since 4007 was 4007, checked in by davestorkey, 8 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).
File size: 15.9 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) new module
7   !!----------------------------------------------------------------------
8#if defined key_bdy 
9   !!----------------------------------------------------------------------
10   !!   'key_bdy' :                    Unstructured Open Boundary Condition
11   !!----------------------------------------------------------------------
12   !!   bdy_orlanski_2d
13   !!   bdy_orlanski_3d
14   !!----------------------------------------------------------------------
15   USE timing          ! Timing
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE bdy_oce         ! ocean open boundary conditions
19   USE phycst          ! physical constants
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE in_out_manager  !
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   bdy_orlanski_2d     ! routine called where?
27   PUBLIC   bdy_orlanski_3d     ! routine called where?
28
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $
32   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo )
37      !!----------------------------------------------------------------------
38      !!                 ***  SUBROUTINE bdy_orlanski_2d  ***
39      !!             
40      !!              - Apply Orlanski radiation condition adaptively to 2D fields:
41      !!                  - radiation plus weak nudging at outflow points
42      !!                  - no radiation and strong nudging at inflow points
43      !!
44      !!
45      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
46      !!----------------------------------------------------------------------
47      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
48      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
49      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phib     ! model before 2D field
50      REAL(wp), DIMENSION(:,:),   INTENT(inout)  ::   phia     ! model after 2D field (to be updated)
51      REAL(wp), DIMENSION(:),     INTENT(in)     ::   phi_ext  ! external forcing data
52      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
53
54      INTEGER  ::   jb                                     ! dummy loop indices
55      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
56      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
57      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
58      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
59      INTEGER  ::   flagu, flagv                           ! short cuts
60      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
61      REAL(wp) ::   zout, zwgt, zdy_centred
62      REAL(wp) ::   zdy_left, zdy_right, zsign_ups
63      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field
64      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdiv ! land/sea mask for x-derivatives
65      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydiv ! land/sea mask for y-derivatives
66      !!----------------------------------------------------------------------
67
68      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_2d')
69
70      ! ----------------------------------!
71      ! Orlanski boundary conditions     :!
72      ! ----------------------------------!
73     
74      SELECT CASE(igrd)
75         CASE(1)
76            pmask => tmask(:,:,1)
77            pmask_xdiv => umask(:,:,1)
78            ii_offset = 0
79            pmask_ydiv => vmask(:,:,1)
80            ij_offset = 0
81         CASE(2)
82            pmask => umask(:,:,1)
83            pmask_xdiv => tmask(:,:,1)
84            ii_offset = 1
85            pmask_ydiv => fmask(:,:,1)
86            ij_offset = 0
87         CASE(3)
88            pmask => vmask(:,:,1)
89            pmask_xdiv => fmask(:,:,1)
90            ii_offset = 0
91            pmask_ydiv => tmask(:,:,1)
92            ij_offset = 1
93         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
94      END SELECT
95      !
96      DO jb = 1, idx%nblenrim(igrd)
97         ii  = idx%nbi(jb,igrd)
98         ij  = idx%nbj(jb,igrd) 
99         flagu = int( idx%flagu(jb,igrd) )
100         flagv = int( idx%flagv(jb,igrd) )
101         !
102         ! calculate positions of b-1 and b-2 points for this rim point
103         ! also (b-1,j-1) and (b-1,j+1) points
104         iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
105         ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
106         !
107         iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
108         ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
109         !
110         iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
111         ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
112         !
113         ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities.
114         ! Mask derivatives to ensure correct land boundary conditions for each variable.
115         ! Centred derivative is calculated as average of "left" and "right" derivatives for
116         ! this reason.
117         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
118         zdx = ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) )                          &
119        &    * ( abs(iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          )   &
120        &      + abs(ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset) ) 
121         zdy_left  = phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1)                  &
122        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          )   & 
123        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
124         zdy_right = phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)                     &
125        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1)                   &
126        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset) ) 
127         zdy_centred = 0.5 * ( zdy_left + zdy_right )
128!!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
129         ! upstream differencing for tangential derivatives
130         zsign_ups = sign( 1., zdt * zdy_centred )
131         zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
132         zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right
133         znor2 = zdx * zdx + zdy * zdy
134         znor2 = max(znor2,rsmall)
135         zcx = zdt * zdx / znor2
136         zcy = zdt * zdy / znor2
137         !
138         ! update boundary value:
139         zout = sign( 1., zcx )
140         zout = 0.5*( zout + abs(zout) )
141         zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
142         ! only apply radiation on outflow points
143         if( ll_npo ) then     !! NPO version !!
144            phia(ii,ij) = (1.-zout) * phib(ii,ij)                                                   &
145           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1) ) / ( 1. + zcx ) 
146         else                  !! full oblique radiation !!
147            zsign_ups = sign( 1., zcy )
148            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
149            phia(ii,ij) = (1.-zout) * phib(ii,ij)                                                   &
150           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1)                         &
151           &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ) - phib(iijm1,ijjm1 ) ) &
152           &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) ) / ( 1. + zcx ) 
153         end if
154!!$         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phib(ii,ij) )
155         phia(ii,ij) = phia(ii,ij) * pmask(ii,ij)
156      END DO
157      !
158      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_2d')
159
160   END SUBROUTINE bdy_orlanski_2d
161
162
163   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, ll_npo )
164      !!----------------------------------------------------------------------
165      !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
166      !!             
167      !!              - Apply Orlanski radiation condition adaptively to 3D fields:
168      !!                  - radiation plus weak nudging at outflow points
169      !!                  - no radiation and strong nudging at inflow points
170      !!
171      !!
172      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
173      !!----------------------------------------------------------------------
174      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
175      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
176      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phib     ! model before 3D field
177      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
178      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phi_ext  ! external forcing data
179      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
180
181      INTEGER  ::   jb, jk                                 ! dummy loop indices
182      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
183      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
184      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
185      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
186      INTEGER  ::   flagu, flagv                           ! short cuts
187      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
188      REAL(wp) ::   zout, zwgt, zdy_centred
189      REAL(wp) ::   zdy_left, zdy_right,  zsign_ups
190      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field
191      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdiv ! land/sea mask for x-derivatives
192      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydiv ! land/sea mask for y-derivatives
193      !!----------------------------------------------------------------------
194
195      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_3d')
196
197      ! ----------------------------------!
198      ! Orlanski boundary conditions     :!
199      ! ----------------------------------!
200     
201      SELECT CASE(igrd)
202         CASE(1)
203            pmask => tmask(:,:,:)
204            pmask_xdiv => umask(:,:,:)
205            ii_offset = 0
206            pmask_ydiv => vmask(:,:,:)
207            ij_offset = 0
208         CASE(2)
209            pmask => umask(:,:,:)
210            pmask_xdiv => tmask(:,:,:)
211            ii_offset = 1
212            pmask_ydiv => fmask(:,:,:)
213            ij_offset = 0
214         CASE(3)
215            pmask => vmask(:,:,:)
216            pmask_xdiv => fmask(:,:,:)
217            ii_offset = 0
218            pmask_ydiv => tmask(:,:,:)
219            ij_offset = 1
220         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
221      END SELECT
222
223      DO jk = 1, jpk
224         !           
225         DO jb = 1, idx%nblenrim(igrd)
226            ii  = idx%nbi(jb,igrd)
227            ij  = idx%nbj(jb,igrd) 
228            flagu = int( idx%flagu(jb,igrd) )
229            flagv = int( idx%flagv(jb,igrd) )
230            !
231            ! calculate positions of b-1 and b-2 points for this rim point
232            ! also (b-1,j-1) and (b-1,j+1) points
233            iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
234            ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
235            !
236            iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
237            ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
238            !
239            iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
240            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
241            !
242            ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities.
243            ! Mask derivatives to ensure correct land boundary conditions for each variable.
244            ! Centred derivative is calculated as average of "left" and "right" derivatives for
245            ! this reason.
246            zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)
247            zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk)                     &
248        &    * ( (iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          ,jk)   &
249        &      + (ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset,jk) ) 
250            zdy_left  = phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk)            &
251        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   & 
252        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset,jk) ) 
253            zdy_right = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk)      &
254        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1          ,jk)   &
255        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset,jk) ) 
256            zdy_centred = 0.5 * ( zdy_left + zdy_right )
257!!$            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk)
258            ! upstream differencing for tangential derivatives
259            zsign_ups = sign( 1., zdt * zdy_centred )
260            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
261            zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right
262            znor2 = zdx * zdx + zdy * zdy
263            znor2 = max(znor2,rsmall)
264            zcx = zdt * zdx / znor2
265            zcy = zdt * zdy / znor2
266            !
267            ! update boundary value:
268            zout = sign( 1., zcx )
269            zout = 0.5*( zout + abs(zout) )
270            zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
271            ! only apply radiation on outflow points
272            if( ll_npo ) then     !! NPO version !!
273               phia(ii,ij,jk) = (1.-zout) * phib(ii,ij,jk)                                                   &
274              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk) ) / ( 1. + zcx ) 
275            else                  !! full oblique radiation !!
276               zsign_ups = sign( 1., zcy )
277               zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
278               phia(ii,ij,jk) = (1.-zout) * phib(ii,ij,jk)                                                   &
279              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk)                      &
280              &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ,jk) - phib(iijm1,ijjm1,jk ) ) &
281              &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk ) ) ) / ( 1. + zcx ) 
282            end if
283!!$            phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) )
284            phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk)
285         END DO
286         !
287      END DO
288
289      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_3d')
290
291   END SUBROUTINE bdy_orlanski_3d
292
293
294#else
295   !!----------------------------------------------------------------------
296   !!   Dummy module                   NO Unstruct Open Boundary Conditions
297   !!----------------------------------------------------------------------
298CONTAINS
299   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
300      WRITE(*,*) 'bdy_orlanski_2d: You should not have seen this print! error?', kt
301   END SUBROUTINE bdy_orlanski_2d
302   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
303      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
304   END SUBROUTINE bdy_orlanski_3d
305#endif
306
307   !!======================================================================
308END MODULE bdylib
Note: See TracBrowser for help on using the repository browser.