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

Last change on this file since 4008 was 4008, checked in by davestorkey, 8 years ago

Uncomment important lines in bdylib.F90!

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.