source: branches/2013/dev_r3891_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90 @ 3902

Last change on this file since 3902 was 3902, checked in by davestorkey, 9 years ago

bug fixes

File size: 8.7 KB
Line 
1MODULE bdydyn2d
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 )
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(in/out) ::   phia     ! model after 2D field (to be updated)
51      REAL(wp), DIMENSION(:),     INTENT(in)     ::   phi_ext  ! external forcing data
52
53      INTEGER  ::   jb                                     ! dummy loop indices
54      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
55      INTEGER, POINTER :: flagu, flagv                     ! short cuts
56      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zout, zwgt  ! intermediate calculations
57      !!----------------------------------------------------------------------
58
59      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_2d')
60
61      ! ----------------------------------!
62      ! Orlanski boundary conditions     :!
63      ! ----------------------------------!
64     
65      !
66      DO jb = 1, idx%nblenrim(igrd)
67         ii  = idx%nbi(jb,igrd)
68         ij  = idx%nbj(jb,igrd) 
69         flagu =>  idx%flagu(jb,igrd)
70         flagv =>  idx%flagu(jb,igrd)
71         !
72         ! calculate positions of b-1 and b-2 points for this rim point
73         ! also (b-1,j-1) and (b-1,j+1) points
74         iibm1 = ii + int(flagu) ; iibm2 = ii + 2*int(flagu) 
75         ijbm1 = ij + int(flagv) ; ijbm2 = ij + 2*int(flagv) 
76         !
77         iibm1jm1 = ii + int(flagu) - abs(int(flagv)) ; iibm1jp1 = ii + int(flagu) + abs(int(flagv)) 
78         ijbm1jm1 = ij + int(flagv) - abs(int(flagu)) ; ijbm1jp1 = ij + int(flagv) + abs(int(flagu)) 
79         !
80         ! normal radiation velocity:
81         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
82         zdx = phia(iibm1,ijbm1) - phia(iibm2,ijbm2)
83         zdy = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
84         znor2 = zdx * zdx + zdy * zdy
85         zcx = zdt * zdx / znor2
86         !
87         ! update boundary value:
88         zout = sign( 1., zcx )
89         zout = 0.5*( zout + abs(zout) )
90         zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
91         ! only apply radiation on outflow points
92         phia(ii,ij) = (1.-zout) * phia(ii,ij) + zout * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1) ) / ( 1. + zcx ) 
93         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phia(ii,ij) ) 
94      END DO
95      !
96
97      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_2d')
98
99   END SUBROUTINE bdy_orlanski_2d
100
101
102   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext )
103      !!----------------------------------------------------------------------
104      !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
105      !!             
106      !!              - Apply Orlanski radiation condition adaptively to 3D fields:
107      !!                  - radiation plus weak nudging at outflow points
108      !!                  - no radiation and strong nudging at inflow points
109      !!
110      !!
111      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)   
112      !!----------------------------------------------------------------------
113      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
114      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
115      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phib     ! model before 3D field
116      REAL(wp), DIMENSION(:,:,:), INTENT(in/out) ::   phia     ! model after 3D field (to be updated)
117      REAL(wp), DIMENSION(:,:),   INTENT(in/out) ::   phi_ext  ! external forcing data
118
119      INTEGER  ::   jb                                     ! dummy loop indices
120      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
121      INTEGER, POINTER :: flagu, flagv                     ! short cuts
122      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zout, zwgt  ! intermediate calculations
123      !!----------------------------------------------------------------------
124
125      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_3d')
126
127      ! ----------------------------------!
128      ! Orlanski boundary conditions     :!
129      ! ----------------------------------!
130     
131      DO ik = 1, jpk
132         !           
133         DO jb = 1, idx%nblenrim(igrd)
134            ii  = idx%nbi(jb,igrd)
135            ij  = idx%nbj(jb,igrd) 
136            flagu =>  idx%flagu(jb,igrd)
137            flagv =>  idx%flagu(jb,igrd)
138            !
139            ! calculate positions of b-1 and b-2 points for this rim point
140            ! also (b-1,j-1) and (b-1,j+1) points
141            iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
142            ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv 
143            !
144            iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
145            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
146            !
147            ! normal radiation velocity:
148            zdt = phia(iibm1,ijbm1,ik) - phib(iibm1,ijbm1,ik)
149            zdx = phia(iibm1,ijbm1,ik) - phia(iibm2,ijbm2,ik)
150            zdy = phib(iibm1jp1,ijbm1jp1,ik) - phib(iibm1jm1,ijbm1jm1,ik)
151            znor2 = zdx * zdx + zdy * zdy
152            zcx = zdt * zdx / znor2
153            !
154            ! update boundary value:
155            zout = sign( 1., zcx )
156            zout = 0.5*( zout + abs(zout) )
157            zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
158            ! only apply radiation on outflow points
159            phia(ii,ij,ik) = (1.-zout) * phia(ii,ij,ik) + zout * ( phib(ii,ij,ik) + zcx*phia(iibm1,ijbm1,ik) ) &
160           &                 / ( 1. + zcx ) 
161            phia(ii,ij,ik) = phia(ii,ij,ik) + zwgt * ( phi_ext(jb,ik) - phia(ii,ij,ik) ) 
162         END DO
163         !
164      END DO
165
166      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_3d')
167
168   END SUBROUTINE bdy_orlanski_3d
169
170
171#else
172   !!----------------------------------------------------------------------
173   !!   Dummy module                   NO Unstruct Open Boundary Conditions
174   !!----------------------------------------------------------------------
175CONTAINS
176   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
177      WRITE(*,*) 'bdy_orlanski_2d: You should not have seen this print! error?', kt
178   END SUBROUTINE bdy_orlanski_2d
179   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
180      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
181   END SUBROUTINE bdy_orlanski_3d
182#endif
183
184   !!======================================================================
185END MODULE bdylib
Note: See TracBrowser for help on using the repository browser.