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 branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

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

Last change on this file since 3994 was 3994, checked in by davestorkey, 11 years ago

Bug fixes: add masking in bdy_orlanski routines and avoid division by
zero error. Also need TARGET attributes in dom_oce.F90.
File size: 12.2 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, mask, 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      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   mask     ! land/sea mask
53      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
54
55      INTEGER  ::   jb                                     ! dummy loop indices
56      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
57      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
58      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
59      INTEGER  ::   flagu, flagv                           ! short cuts
60      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
61      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
62      !!----------------------------------------------------------------------
63
64      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_2d')
65
66      ! ----------------------------------!
67      ! Orlanski boundary conditions     :!
68      ! ----------------------------------!
69     
70      !
71      DO jb = 1, idx%nblenrim(igrd)
72         ii  = idx%nbi(jb,igrd)
73         ij  = idx%nbj(jb,igrd) 
74         flagu = int( idx%flagu(jb,igrd) )
75         flagv = int( idx%flagv(jb,igrd) )
76         !
77         ! calculate positions of b-1 and b-2 points for this rim point
78         ! also (b-1,j-1) and (b-1,j+1) points
79         iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
80         ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
81         !
82         iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
83         ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
84         !
85         iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
86         ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
87         !
88         ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
89         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
90         zdx = phia(iibm1,ijbm1) - phia(iibm2,ijbm2)
91         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
92         ! upstream differencing for tangential derivatives
93         zsign_ups = sign( 1., zdt * zdy_centred )
94         zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
95         zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) &
96       &        + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) )
97         znor2 = zdx * zdx + zdy * zdy
98         znor2 = max(znor2,rsmall)
99         zcx = zdt * zdx / znor2
100         zcy = zdt * zdy / znor2
101         !
102         ! update boundary value:
103         zout = sign( 1., zcx )
104         zout = 0.5*( zout + abs(zout) )
105         zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
106         ! only apply radiation on outflow points
107         if( ll_npo ) then     !! NPO version !!
108            phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   &
109           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1) ) / ( 1. + zcx ) 
110         else                  !! full oblique radiation !!
111            zsign_ups = sign( 1., zcy )
112            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
113            phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   &
114           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1)                         &
115           &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ) - phib(iijm1,ijjm1 ) ) &
116           &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) ) / ( 1. + zcx ) 
117         end if
118         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phia(ii,ij) ) 
119         phia(ii,ij) = phia(ii,ij) * mask(ii,ij)
120      END DO
121      !
122      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_2d')
123
124   END SUBROUTINE bdy_orlanski_2d
125
126
127   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, mask, ll_npo )
128      !!----------------------------------------------------------------------
129      !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
130      !!             
131      !!              - Apply Orlanski radiation condition adaptively to 3D fields:
132      !!                  - radiation plus weak nudging at outflow points
133      !!                  - no radiation and strong nudging at inflow points
134      !!
135      !!
136      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
137      !!----------------------------------------------------------------------
138      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
139      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
140      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phib     ! model before 3D field
141      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
142      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phi_ext  ! external forcing data
143      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   mask     ! land/sea mask
144      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
145
146      INTEGER  ::   jb, jk                                 ! dummy loop indices
147      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
148      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
149      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
150      INTEGER  ::   flagu, flagv                           ! short cuts
151      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
152      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
153      !!----------------------------------------------------------------------
154
155      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_3d')
156
157      ! ----------------------------------!
158      ! Orlanski boundary conditions     :!
159      ! ----------------------------------!
160     
161      DO jk = 1, jpk
162         !           
163         DO jb = 1, idx%nblenrim(igrd)
164            ii  = idx%nbi(jb,igrd)
165            ij  = idx%nbj(jb,igrd) 
166            flagu = int( idx%flagu(jb,igrd) )
167            flagv = int( idx%flagv(jb,igrd) )
168            !
169            ! calculate positions of b-1 and b-2 points for this rim point
170            ! also (b-1,j-1) and (b-1,j+1) points
171            iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
172            ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
173            !
174            iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
175            ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
176            !
177            iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
178            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
179            !
180            ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
181            zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)
182            zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk)
183            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk)
184            ! upstream differencing for tangential derivatives
185            zsign_ups = sign( 1., zdt * zdy_centred )
186            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
187            zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) &
188           &    + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) )
189            znor2 = zdx * zdx + zdy * zdy
190            znor2 = max(znor2,rsmall)
191            zcx = zdt * zdx / znor2
192            zcy = zdt * zdy / znor2
193            !
194            ! update boundary value:
195            zout = sign( 1., zcx )
196            zout = 0.5*( zout + abs(zout) )
197            zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
198            ! only apply radiation on outflow points
199            if( ll_npo ) then     !! NPO version !!
200               phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   &
201              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk) ) / ( 1. + zcx ) 
202            else                  !! full oblique radiation !!
203               zsign_ups = sign( 1., zcy )
204               zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
205               phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   &
206              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk)                      &
207              &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ,jk) - phib(iijm1,ijjm1,jk ) ) &
208              &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk ) ) ) / ( 1. + zcx ) 
209            end if
210            phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phia(ii,ij,jk) ) 
211            phia(ii,ij,jk) = phia(ii,ij,jk) * mask(ii,ij,jk)
212         END DO
213         !
214      END DO
215
216      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_3d')
217
218   END SUBROUTINE bdy_orlanski_3d
219
220
221#else
222   !!----------------------------------------------------------------------
223   !!   Dummy module                   NO Unstruct Open Boundary Conditions
224   !!----------------------------------------------------------------------
225CONTAINS
226   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
227      WRITE(*,*) 'bdy_orlanski_2d: You should not have seen this print! error?', kt
228   END SUBROUTINE bdy_orlanski_2d
229   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
230      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
231   END SUBROUTINE bdy_orlanski_3d
232#endif
233
234   !!======================================================================
235END MODULE bdylib
Note: See TracBrowser for help on using the repository browser.