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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 16.6 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   !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   div_cur    : Compute the horizontal divergence and relative
24   !!                vorticity fields
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
29   USE sbcrnf          ! river runoff
30   USE sbcisf          ! ice shelf
31   USE cla             ! cross land advection             (cla_div routine)
32   USE in_out_manager  ! I/O manager
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE lib_mpp         ! MPP library
35   USE wrk_nemo        ! Memory Allocation
36   USE timing          ! Timing
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53#if defined key_noslip_accurate
54   !!----------------------------------------------------------------------
55   !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast
56   !!----------------------------------------------------------------------
57
58   SUBROUTINE div_cur( kt )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE div_cur  ***
61      !!
62      !! ** Purpose :   compute the horizontal divergence and the relative
63      !!              vorticity at before and now time-step
64      !!
65      !! ** Method  : I.  divergence :
66      !!         - save the divergence computed at the previous time-step
67      !!      (note that the Asselin filter has not been applied on hdivb)
68      !!         - compute the now divergence given by :
69      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
70      !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf)
71      !!      and cross land flow (div_cla)
72      !!              II. vorticity :
73      !!         - save the curl computed at the previous time-step
74      !!            rotb = rotn
75      !!      (note that the Asselin time filter has not been applied to rotb)
76      !!         - compute the now curl in tensorial formalism:
77      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
78      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
79      !!      the no-slip boundary condition is computed using Schchepetkin
80      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
81      !!      For example, along east coast, the one-sided finite difference
82      !!      approximation used for di[v] is:
83      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
84      !!
85      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
86      !!              - update rotb , rotn , the before & now rel. vorticity
87      !!----------------------------------------------------------------------
88      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
89      !
90      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices
91      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer
92      REAL(wp) ::  zraur, zdep              ! local scalar
93      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace
94      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace
95      !!----------------------------------------------------------------------
96      !
97      IF( nn_timing == 1 )  CALL timing_start('div_cur')
98      !
99      CALL wrk_alloc( jpi  , jpj+2, zwu  )
100      CALL wrk_alloc( jpi+2, jpj  , zwv  )
101      !
102      IF( kt == nit000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity'
105         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
106         IF(lwp .AND. lflush) CALL flush(numout)
107      ENDIF
108
109      !                                                ! ===============
110      DO jk = 1, jpkm1                                 ! Horizontal slab
111         !                                             ! ===============
112         !
113         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
114         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
115         !
116         !                                             ! --------
117         ! Horizontal divergence                       !   div
118         !                                             ! --------
119         DO jj = 2, jpjm1
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               hdivn(ji,jj,jk) =   &
122                  (  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)       &
123                   + 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)  )    &
124                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
125            END DO
126         END DO
127
128         IF( .NOT. AGRIF_Root() ) THEN
129            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
130            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
131            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
132            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
133         ENDIF
134
135         !                                             ! --------
136         ! relative vorticity                          !   rot
137         !                                             ! --------
138         ! contravariant velocity (extended for lateral b.c.)
139         ! inside the model domain
140         DO jj = 1, jpj
141            DO ji = 1, jpi
142               zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk)
143               zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk)
144            END DO 
145         END DO 
146 
147         ! East-West boundary conditions
148         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
149            zwv(  0  ,:) = zwv(jpi-2,:)
150            zwv( -1  ,:) = zwv(jpi-3,:)
151            zwv(jpi+1,:) = zwv(  3  ,:)
152            zwv(jpi+2,:) = zwv(  4  ,:)
153         ELSE
154            zwv(  0  ,:) = 0.e0
155            zwv( -1  ,:) = 0.e0
156            zwv(jpi+1,:) = 0.e0
157            zwv(jpi+2,:) = 0.e0
158         ENDIF
159
160         ! North-South boundary conditions
161         IF( nperio == 3 .OR. nperio == 4 ) THEN
162            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
163            zwu(jpi,jpj+1) = 0.e0
164            zwu(jpi,jpj+2) = 0.e0
165            DO ji = 1, jpi-1
166               iju = jpi - ji + 1
167               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
168               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
169            END DO
170         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
171            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\
172            zwu(jpi,jpj+1) = 0.e0
173            zwu(jpi,jpj+2) = 0.e0
174            DO ji = 1, jpi-1
175               iju = jpi - ji
176               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
177               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
178               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
179            END DO
180            DO ji = -1, jpi+2
181               ijt = jpi - ji + 1
182               zwv(ji,jpj) = - zwv(ijt,jpj-2)
183            END DO
184            DO ji = jpi/2+1, jpi+2
185               ijt = jpi - ji + 1
186               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
187            END DO
188         ELSE
189            ! closed
190            zwu(:,jpj+1) = 0.e0
191            zwu(:,jpj+2) = 0.e0
192         ENDIF
193
194         ! relative vorticity (vertical component of the velocity curl)
195         DO jj = 1, jpjm1
196            DO ji = 1, fs_jpim1   ! vector opt.
197               rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
198                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
199            END DO
200         END DO
201
202         ! second order accurate scheme along straight coast
203         DO jl = 1, npcoa(1,jk)
204            ii = nicoa(jl,1,jk)
205            ij = njcoa(jl,1,jk)
206            rotn(ii,ij,jk) = 1. / ( e1f(ii,ij) * e2f(ii,ij) )   &
207                           * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) )
208         END DO
209         DO jl = 1, npcoa(2,jk)
210            ii = nicoa(jl,2,jk)
211            ij = njcoa(jl,2,jk)
212            rotn(ii,ij,jk) = 1./(e1f(ii,ij)*e2f(ii,ij))   &
213               *(-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij))
214         END DO
215         DO jl = 1, npcoa(3,jk)
216            ii = nicoa(jl,3,jk)
217            ij = njcoa(jl,3,jk)
218            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
219               * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) )
220         END DO
221         DO jl = 1, npcoa(4,jk)
222            ii = nicoa(jl,4,jk)
223            ij = njcoa(jl,4,jk)
224            rotn(ii,ij,jk) = -1. / ( e1f(ii,ij)*e2f(ii,ij) )   &
225               * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
226         END DO
227         !                                             ! ===============
228      END DO                                           !   End of slab
229      !                                                ! ===============
230
231      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
232      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
233      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
234     
235      ! 4. Lateral boundary conditions on hdivn and rotn
236      ! ---------------------------------=======---======
237      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
238      !
239      CALL wrk_dealloc( jpi  , jpj+2, zwu )
240      CALL wrk_dealloc( jpi+2, jpj  , zwv )
241      !
242      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
243      !
244   END SUBROUTINE div_cur
245   
246#else
247   !!----------------------------------------------------------------------
248   !!   Default option                           2nd order centered schemes
249   !!----------------------------------------------------------------------
250
251   SUBROUTINE div_cur( kt )
252      !!----------------------------------------------------------------------
253      !!                  ***  ROUTINE div_cur  ***
254      !!                   
255      !! ** Purpose :   compute the horizontal divergence and the relative
256      !!      vorticity at before and now time-step
257      !!
258      !! ** Method  : - Divergence:
259      !!      - save the divergence computed at the previous time-step
260      !!      (note that the Asselin filter has not been applied on hdivb)
261      !!      - compute the now divergence given by :
262      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
263      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)
264      !!              - Relavtive Vorticity :
265      !!      - save the curl computed at the previous time-step (rotb = rotn)
266      !!      (note that the Asselin time filter has not been applied to rotb)
267      !!      - compute the now curl in tensorial formalism:
268      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
269      !!      Note: Coastal boundary condition: lateral friction set through
270      !!      the value of fmask along the coast (see dommsk.F90) and shlat
271      !!      (namelist parameter)
272      !!
273      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
274      !!              - update rotb , rotn , the before & now rel. vorticity
275      !!----------------------------------------------------------------------
276      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
277      !
278      INTEGER  ::   ji, jj, jk    ! dummy loop indices
279      REAL(wp) ::   zraur, zdep   ! local scalars
280      !!----------------------------------------------------------------------
281      !
282      IF( nn_timing == 1 )  CALL timing_start('div_cur')
283      !
284      IF( kt == nit000 ) THEN
285         IF(lwp) WRITE(numout,*)
286         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and'
287         IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity'
288         IF(lwp .AND. lflush) CALL flush(numout)
289      ENDIF
290
291      !                                                ! ===============
292      DO jk = 1, jpkm1                                 ! Horizontal slab
293         !                                             ! ===============
294         !
295         hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
296         rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
297         !
298         !                                             ! --------
299         ! Horizontal divergence                       !   div
300         !                                             ! --------
301         DO jj = 2, jpjm1
302            DO ji = fs_2, fs_jpim1   ! vector opt.
303               hdivn(ji,jj,jk) =   &
304                  (  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)       &
305                   + 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)  )    &
306                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
307            END DO 
308         END DO 
309
310         IF( .NOT. AGRIF_Root() ) THEN
311            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east
312            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west
313            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north
314            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south
315         ENDIF
316
317         !                                             ! --------
318         ! relative vorticity                          !   rot
319         !                                             ! --------
320         DO jj = 1, jpjm1
321            DO ji = 1, fs_jpim1   ! vector opt.
322               rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
323                  &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
324                  &           * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
325            END DO
326         END DO
327         !                                             ! ===============
328      END DO                                           !   End of slab
329      !                                                ! ===============
330
331      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
332      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
333      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
334      !
335      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
336      !
337      IF( nn_timing == 1 )  CALL timing_stop('div_cur')
338      !
339   END SUBROUTINE div_cur
340   
341#endif
342   !!======================================================================
343END MODULE divcur
Note: See TracBrowser for help on using the repository browser.