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