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.
divcur.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

  1. Initial restructuring.
  2. Use fldread to read open boundary data.
  • Property svn:keywords set to Id
File size: 15.8 KB
Line 
1MODULE divcur
2   !!==============================================================================
3   !!                       ***  MODULE  divcur  ***
4   !! Ocean diagnostic variable : horizontal divergence and relative vorticity
5   !!==============================================================================
6   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
9   !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
10   !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
11   !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
12   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
13   !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
14   !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
15   !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
16   !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
17   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
18   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   div_cur    : Compute the horizontal divergence and relative
23   !!                vorticity fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
28   USE sbcrnf          ! river runoff
29   USE cla             ! cross land advection             (cla_div routine)
30   USE in_out_manager  ! I/O manager
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE lib_mpp         ! MPP library
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49#if defined key_noslip_accurate
50   !!----------------------------------------------------------------------
51   !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast
52   !!----------------------------------------------------------------------
53
54   SUBROUTINE div_cur( kt )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE div_cur  ***
57      !!
58      !! ** Purpose :   compute the horizontal divergence and the relative
59      !!              vorticity at before and now time-step
60      !!
61      !! ** Method  : I.  divergence :
62      !!         - save the divergence computed at the previous time-step
63      !!      (note that the Asselin filter has not been applied on hdivb)
64      !!         - compute the now divergence given by :
65      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
66      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
67      !!              II. vorticity :
68      !!         - save the curl computed at the previous time-step
69      !!            rotb = rotn
70      !!      (note that the Asselin time filter has not been applied to rotb)
71      !!         - compute the now curl in tensorial formalism:
72      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
73      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
74      !!      the no-slip boundary condition is computed using Schchepetkin
75      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
76      !!      For example, along east coast, the one-sided finite difference
77      !!      approximation used for di[v] is:
78      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
79      !!
80      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
81      !!              - update rotb , rotn , the before & now rel. vorticity
82      !!----------------------------------------------------------------------
83      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
84      !
85      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwu   ! specific 2D workspace
86      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwv   ! specific 2D workspace
87      !
88      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
89      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
90      REAL(wp) ::  zraur, zdep              ! local scalar
91      !!----------------------------------------------------------------------
92
93      IF( kt == nit000 ) THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
96         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
97         !
98         ALLOCATE( zwu( jpi, 1:jpj+2) , zwv(-1:jpi+2, jpj) , STAT=ierr )
99         IF( lk_mpp    )   CALL mpp_sum( ierr )
100         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'div_cur : unable to allocate arrays' )
101      ENDIF
102
103      !                                                ! ===============
104      DO jk = 1, jpkm1                                 ! Horizontal slab
105         !                                             ! ===============
106         !
107         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
108         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
109         !
110         !                                             ! --------
111         ! Horizontal divergence                       !   div
112         !                                             ! --------
113         DO jj = 2, jpjm1
114            DO ji = fs_2, fs_jpim1   ! vector opt.
115               hdivn(ji,jj,jk) =   &
116                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
117                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
118                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
119            END DO
120         END DO
121
122         IF( .NOT. AGRIF_Root() ) THEN
123            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
124            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
125            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
126            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
127         ENDIF
128
129         !                                             ! --------
130         ! relative vorticity                          !   rot
131         !                                             ! --------
132         ! contravariant velocity (extended for lateral b.c.)
133         ! inside the model domain
134         DO jj = 1, jpj
135            DO ji = 1, jpi
136               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
137               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
138            END DO 
139         END DO 
140 
141         ! East-West boundary conditions
142         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
143            zwv(  0  ,:) = zwv(jpi-2,:)
144            zwv( -1  ,:) = zwv(jpi-3,:)
145            zwv(jpi+1,:) = zwv(  3  ,:)
146            zwv(jpi+2,:) = zwv(  4  ,:)
147         ELSE
148            zwv(  0  ,:) = 0.e0
149            zwv( -1  ,:) = 0.e0
150            zwv(jpi+1,:) = 0.e0
151            zwv(jpi+2,:) = 0.e0
152         ENDIF
153
154         ! North-South boundary conditions
155         IF( nperio == 3 .OR. nperio == 4 ) THEN
156            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
157            zwu(jpi,jpj+1) = 0.e0
158            zwu(jpi,jpj+2) = 0.e0
159            DO ji = 1, jpi-1
160               iju = jpi - ji + 1
161               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
162               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
163            END DO
164         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
165            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
166            zwu(jpi,jpj+1) = 0.e0
167            zwu(jpi,jpj+2) = 0.e0
168            DO ji = 1, jpi-1
169               iju = jpi - ji
170               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
171               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
172               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
173            END DO
174            DO ji = -1, jpi+2
175               ijt = jpi - ji + 1
176               zwv(ji,jpj) = - zwv(ijt,jpj-2)
177            END DO
178            DO ji = jpi/2+1, jpi+2
179               ijt = jpi - ji + 1
180               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
181            END DO
182         ELSE
183            ! closed
184            zwu(:,jpj+1) = 0.e0
185            zwu(:,jpj+2) = 0.e0
186         ENDIF
187
188         ! relative vorticity (vertical component of the velocity curl)
189         DO jj = 1, jpjm1
190            DO ji = 1, fs_jpim1   ! vector opt.
191               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
192                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
193            END DO
194         END DO
195
196         ! second order accurate scheme along straight coast
197         DO jl = 1, npcoa(1,jk)
198            ii = nicoa(jl,1,jk)
199            ij = njcoa(jl,1,jk)
200            rotn(ii,ij,jk) = 1. / ( e1f(ii,ij) * e2f(ii,ij) )   &
201                           * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) )
202         END DO
203         DO jl = 1, npcoa(2,jk)
204            ii = nicoa(jl,2,jk)
205            ij = njcoa(jl,2,jk)
206            rotn(ii,ij,jk) = 1./(e1f(ii,ij)*e2f(ii,ij))   &
207               *(-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
208         END DO
209         DO jl = 1, npcoa(3,jk)
210            ii = nicoa(jl,3,jk)
211            ij = njcoa(jl,3,jk)
212            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
213               * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
214         END DO
215         DO jl = 1, npcoa(4,jk)
216            ii = nicoa(jl,4,jk)
217            ij = njcoa(jl,4,jk)
218            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
219               * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
220         END DO
221         !                                             ! ===============
222      END DO                                           !   End of slab
223      !                                                ! ===============
224
225      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
226      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
227     
228      ! 4. Lateral boundary conditions on hdivn and rotn
229      ! ---------------------------------=======---======
230      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
231      !
232   END SUBROUTINE div_cur
233   
234#else
235   !!----------------------------------------------------------------------
236   !!   Default option                           2nd order centered schemes
237   !!----------------------------------------------------------------------
238
239   SUBROUTINE div_cur( kt )
240      !!----------------------------------------------------------------------
241      !!                  ***  ROUTINE div_cur  ***
242      !!                   
243      !! ** Purpose :   compute the horizontal divergence and the relative
244      !!      vorticity at before and now time-step
245      !!
246      !! ** Method  : - Divergence:
247      !!      - save the divergence computed at the previous time-step
248      !!      (note that the Asselin filter has not been applied on hdivb)
249      !!      - compute the now divergence given by :
250      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
251      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
252      !!              - Relavtive Vorticity :
253      !!      - save the curl computed at the previous time-step (rotb = rotn)
254      !!      (note that the Asselin time filter has not been applied to rotb)
255      !!      - compute the now curl in tensorial formalism:
256      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
257      !!      Note: Coastal boundary condition: lateral friction set through
258      !!      the value of fmask along the coast (see dommsk.F90) and shlat
259      !!      (namelist parameter)
260      !!
261      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
262      !!              - update rotb , rotn , the before & now rel. vorticity
263      !!----------------------------------------------------------------------
264      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
265      !
266      INTEGER  ::   ji, jj, jk    ! dummy loop indices
267      REAL(wp) ::   zraur, zdep   ! local scalars
268      !!----------------------------------------------------------------------
269
270      IF( kt == nit000 ) THEN
271         IF(lwp) WRITE(numout,*)
272         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and'
273         IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity'
274      ENDIF
275
276      !                                                ! ===============
277      DO jk = 1, jpkm1                                 ! Horizontal slab
278         !                                             ! ===============
279         !
280         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
281         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
282         !
283         !                                             ! --------
284         ! Horizontal divergence                       !   div
285         !                                             ! --------
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1   ! vector opt.
288               hdivn(ji,jj,jk) =   &
289                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       &
290                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    &
291                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
292            END DO 
293         END DO 
294
295         IF( .NOT. AGRIF_Root() ) THEN
296            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
297            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
298            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
299            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
300         ENDIF
301
302         !                                             ! --------
303         ! relative vorticity                          !   rot
304         !                                             ! --------
305         DO jj = 1, jpjm1
306            DO ji = 1, fs_jpim1   ! vector opt.
307               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
308                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
309                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
310            END DO
311         END DO
312         !                                             ! ===============
313      END DO                                           !   End of slab
314      !                                                ! ===============
315
316      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
317      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
318      !
319      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
320      !
321   END SUBROUTINE div_cur
322   
323#endif
324   !!======================================================================
325END MODULE divcur
Note: See TracBrowser for help on using the repository browser.