Line data Source code
1 : module viscosity_mod
2 : !
3 : ! This module should be renamed "global_deriv_mod.F90"
4 : !
5 : ! It is a collection of derivative operators that must be applied to the field
6 : ! over the sphere (as opposed to derivative operators that can be applied element
7 : ! by element)
8 : !
9 : !
10 : use shr_kind_mod, only: r8=>shr_kind_r8
11 : use thread_mod, only: max_num_threads, omp_get_num_threads
12 : use dimensions_mod, only: np, nc, nlev,nlevp, qsize,nelemd
13 : use hybrid_mod, only: hybrid_t, get_loop_ranges, config_thread_region
14 : use parallel_mod, only: parallel_t
15 : use element_mod, only: element_t
16 : use derivative_mod, only: derivative_t, laplace_sphere_wk, vlaplace_sphere_wk, vorticity_sphere, derivinit, divergence_sphere
17 : use edgetype_mod, only: EdgeBuffer_t, EdgeDescriptor_t
18 : use edge_mod, only: edgevpack, edgevunpack, edgeVunpackmin, edgeSunpackmin, &
19 : edgeVunpackmax, initEdgeBuffer, FreeEdgeBuffer, edgeSunpackmax, edgeSpack
20 : use bndry_mod, only: bndry_exchange, bndry_exchange_start,bndry_exchange_finish
21 : use control_mod, only: hypervis_scaling, nu, nu_div
22 : use thread_mod, only: vert_num_threads
23 :
24 : implicit none
25 : save
26 :
27 : public :: biharmonic_wk_scalar
28 : public :: biharmonic_wk_omega
29 : public :: neighbor_minmax, neighbor_minmax_start,neighbor_minmax_finish
30 :
31 : !
32 : ! compute vorticity/divergence and then project to make continious
33 : ! high-level routines uses only for I/O
34 : public :: compute_zeta_C0
35 : public :: compute_div_C0
36 :
37 : interface compute_zeta_C0
38 : module procedure compute_zeta_C0_hybrid ! hybrid version
39 : module procedure compute_zeta_C0_par ! single threaded
40 : end interface compute_zeta_C0
41 : interface compute_div_C0
42 : module procedure compute_div_C0_hybrid
43 : module procedure compute_div_C0_par
44 : end interface compute_div_C0
45 :
46 : public :: compute_zeta_C0_contra ! for older versions of sweq which carry
47 : public :: compute_div_C0_contra ! velocity around in contra-coordinates
48 :
49 : type (EdgeBuffer_t) :: edge1
50 :
51 : CONTAINS
52 :
53 6649344 : subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend)
54 : use derivative_mod, only : subcell_Laplace_fluxes
55 : use dimensions_mod, only : use_cslam, nu_div_lev,nu_lev
56 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
57 : ! compute weak biharmonic operator
58 : ! input: h,v (stored in elem()%, in lat-lon coordinates
59 : ! output: ttens,vtens overwritten with weak biharmonic of h,v (output in lat-lon coordinates)
60 : !
61 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
62 : type (hybrid_t) , intent(in) :: hybrid
63 : type (element_t) , intent(inout), target :: elem(:)
64 : integer , intent(in) :: nt,nets,nete
65 : integer , intent(in) :: kbeg, kend
66 : real (kind=r8), intent(out), dimension(nc,nc,4,nlev,nets:nete) :: dpflux
67 : real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: vtens
68 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: ttens,dptens
69 : type (EdgeBuffer_t) , intent(inout) :: edge3
70 : type (derivative_t) , intent(in) :: deriv
71 : ! local
72 : integer :: i,j,k,kptr,ie,kblk
73 : ! real (kind=r8), dimension(:,:), pointer :: rspheremv
74 : real (kind=r8), dimension(np,np) :: tmp
75 : real (kind=r8), dimension(np,np) :: tmp2
76 : real (kind=r8), dimension(np,np,2) :: v
77 :
78 : real (kind=r8), dimension(np,np,nlev) :: lap_p_wk
79 : real (kind=r8), dimension(np,np,nlevp) :: T_i
80 :
81 :
82 : real (kind=r8) :: nu_ratio1, nu_ratio2, dp_thresh
83 : logical var_coef1
84 :
85 6649344 : kblk = kend - kbeg + 1
86 :
87 >23049*10^7 : if (use_cslam) dpflux = 0
88 : !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
89 : !so tensor is only used on second call to laplace_sphere_wk
90 6649344 : var_coef1 = .true.
91 6649344 : if(hypervis_scaling > 0) var_coef1 = .false.
92 6649344 : dp_thresh=.025_r8 ! tunable coefficient
93 53402544 : do ie=nets,nete
94 : !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
95 4394800800 : do k=kbeg,kend
96 4348047600 : nu_ratio1=1
97 4348047600 : nu_ratio2=1
98 4348047600 : if (nu_div_lev(k)/=nu_lev(k)) then
99 4348047600 : if(hypervis_scaling /= 0) then
100 : ! we have a problem with the tensor in that we cant seperate
101 : ! div and curl components. So we do, with tensor V:
102 : ! nu * (del V del ) * ( nu_ratio * grad(div) - curl(curl))
103 0 : nu_ratio1=nu_div_lev(k)/nu_lev(k)
104 : nu_ratio2=1
105 : else
106 4348047600 : nu_ratio1=sqrt(nu_div_lev(k)/nu_lev(k))
107 4348047600 : nu_ratio2=sqrt(nu_div_lev(k)/nu_lev(k))
108 : endif
109 : endif
110 :
111 91308999600 : tmp=elem(ie)%state%T(:,:,k,nt)-elem(ie)%derived%T_ref(:,:,k)
112 4348047600 : call laplace_sphere_wk(tmp,deriv,elem(ie),ttens(:,:,k,ie),var_coef=var_coef1)
113 :
114 91308999600 : tmp=elem(ie)%state%dp3d(:,:,k,nt)-elem(ie)%derived%dp_ref(:,:,k)
115 4348047600 : call laplace_sphere_wk(tmp,deriv,elem(ie),dptens(:,:,k,ie),var_coef=var_coef1)
116 :
117 4348047600 : call vlaplace_sphere_wk(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),.true.,vtens(:,:,:,k,ie), &
118 8742848400 : var_coef=var_coef1,nu_ratio=nu_ratio1)
119 : enddo
120 :
121 46753200 : kptr = kbeg - 1
122 46753200 : call edgeVpack(edge3,ttens(:,:,kbeg:kend,ie),kblk,kptr,ie)
123 :
124 46753200 : kptr = kbeg - 1 + nlev
125 91355752800 : call edgeVpack(edge3,vtens(:,:,1,kbeg:kend,ie),kblk,kptr,ie)
126 :
127 46753200 : kptr = kbeg - 1 + 2*nlev
128 91355752800 : call edgeVpack(edge3,vtens(:,:,2,kbeg:kend,ie),kblk,kptr,ie)
129 :
130 46753200 : kptr = kbeg - 1 + 3*nlev
131 53402544 : call edgeVpack(edge3,dptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
132 : enddo
133 :
134 6649344 : call bndry_exchange(hybrid,edge3,location='biharmonic_wk_dp3d')
135 :
136 53402544 : do ie=nets,nete
137 : !CLEAN rspheremv => elem(ie)%rspheremp(:,:)
138 :
139 46753200 : kptr = kbeg - 1
140 46753200 : call edgeVunpack(edge3,ttens(:,:,kbeg:kend,ie),kblk,kptr,ie)
141 :
142 46753200 : kptr = kbeg - 1 + nlev
143 >18266*10^7 : call edgeVunpack(edge3,vtens(:,:,1,kbeg:kend,ie),kblk,kptr,ie)
144 :
145 46753200 : kptr = kbeg - 1 + 2*nlev
146 >18266*10^7 : call edgeVunpack(edge3,vtens(:,:,2,kbeg:kend,ie),kblk,kptr,ie)
147 :
148 46753200 : kptr = kbeg - 1 + 3*nlev
149 46753200 : call edgeVunpack(edge3,dptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
150 :
151 46753200 : if (use_cslam) then
152 4394800800 : do k=1,nlev
153 : !CLEAN tmp(:,:)= rspheremv(:,:)*dptens(:,:,k,ie)
154 91308999600 : tmp(:,:)= elem(ie)%rspheremp(:,:)*dptens(:,:,k,ie)
155 4394800800 : call subcell_Laplace_fluxes(tmp, deriv, elem(ie), np, nc,dpflux(:,:,:,k,ie))
156 : enddo
157 : endif
158 :
159 : ! apply inverse mass matrix, then apply laplace again
160 : !$omp parallel do num_threads(vert_num_threads) private(k,v,tmp,tmp2)
161 4401450144 : do k=kbeg,kend
162 : !CLEAN tmp(:,:)=rspheremv(:,:)*ttens(:,:,k,ie)
163 91308999600 : tmp(:,:)=elem(ie)%rspheremp(:,:)*ttens(:,:,k,ie)
164 4348047600 : call laplace_sphere_wk(tmp,deriv,elem(ie),ttens(:,:,k,ie),var_coef=.true.)
165 : !CLEAN tmp2(:,:)=rspheremv(:,:)*dptens(:,:,k,ie)
166 91308999600 : tmp2(:,:)=elem(ie)%rspheremp(:,:)*dptens(:,:,k,ie)
167 4348047600 : call laplace_sphere_wk(tmp2,deriv,elem(ie),dptens(:,:,k,ie),var_coef=.true.)
168 : !CLEAN v(:,:,1)=rspheremv(:,:)*vtens(:,:,1,k,ie)
169 : !CLEAN v(:,:,2)=rspheremv(:,:)*vtens(:,:,2,k,ie)
170 :
171 91308999600 : v(:,:,1)=elem(ie)%rspheremp(:,:)*vtens(:,:,1,k,ie)
172 91308999600 : v(:,:,2)=elem(ie)%rspheremp(:,:)*vtens(:,:,2,k,ie)
173 : call vlaplace_sphere_wk(v(:,:,:),deriv,elem(ie),.true.,vtens(:,:,:,k,ie), &
174 4394800800 : var_coef=.true.,nu_ratio=nu_ratio2)
175 :
176 : enddo
177 : enddo
178 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
179 6649344 : end subroutine biharmonic_wk_dp3d
180 :
181 :
182 1112832 : subroutine biharmonic_wk_omega(elem,ptens,deriv,edge3,hybrid,nets,nete,kbeg,kend)
183 : type (hybrid_t) , intent(in) :: hybrid
184 : type (element_t) , intent(inout), target :: elem(:)
185 : integer , intent(in) :: nets,nete
186 : integer , intent(in) :: kbeg, kend
187 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: ptens
188 : type (EdgeBuffer_t) , intent(inout) :: edge3
189 : type (derivative_t) , intent(in) :: deriv
190 :
191 : ! local
192 : integer :: i,j,k,kptr,ie,kblk
193 1112832 : real (kind=r8), dimension(:,:), pointer :: rspheremv
194 : real (kind=r8), dimension(np,np) :: tmp
195 : real (kind=r8), dimension(np,np) :: tmp2
196 : real (kind=r8), dimension(np,np,2) :: v
197 : real (kind=r8) :: nu_ratio1, nu_ratio2
198 : logical var_coef1
199 :
200 1112832 : kblk = kend - kbeg + 1
201 :
202 : !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
203 : !so tensor is only used on second call to laplace_sphere_wk
204 1112832 : var_coef1 = .true.
205 0 : if(hypervis_scaling > 0) var_coef1 = .false.
206 :
207 1112832 : nu_ratio1=1
208 1112832 : nu_ratio2=1
209 :
210 8937432 : do ie=nets,nete
211 :
212 : !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
213 735512400 : do k=kbeg,kend
214 15281443800 : tmp=elem(ie)%derived%omega(:,:,k)
215 735512400 : call laplace_sphere_wk(tmp,deriv,elem(ie),ptens(:,:,k,ie),var_coef=var_coef1)
216 : enddo
217 :
218 7824600 : kptr = kbeg - 1
219 8937432 : call edgeVpack(edge3,ptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
220 : enddo
221 :
222 1112832 : call bndry_exchange(hybrid,edge3,location='biharmonic_wk_omega')
223 :
224 8937432 : do ie=nets,nete
225 7824600 : rspheremv => elem(ie)%rspheremp(:,:)
226 :
227 7824600 : kptr = kbeg - 1
228 7824600 : call edgeVunpack(edge3,ptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
229 :
230 : ! apply inverse mass matrix, then apply laplace again
231 : !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
232 736625232 : do k=kbeg,kend
233 15281443800 : tmp(:,:)=rspheremv(:,:)*ptens(:,:,k,ie)
234 735512400 : call laplace_sphere_wk(tmp,deriv,elem(ie),ptens(:,:,k,ie),var_coef=.true.)
235 : enddo
236 : enddo
237 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238 7762176 : end subroutine biharmonic_wk_omega
239 :
240 :
241 0 : subroutine biharmonic_wk_scalar(elem,qtens,deriv,edgeq,hybrid,nets,nete)
242 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
243 : ! compute weak biharmonic operator
244 : ! input: qtens = Q
245 : ! output: qtens = weak biharmonic of Q
246 : !
247 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
248 : type (hybrid_t) , intent(in) :: hybrid
249 : type (element_t) , intent(inout), target :: elem(:)
250 : integer :: nets,nete
251 : real (kind=r8), dimension(np,np,nlev,qsize,nets:nete) :: qtens
252 : type (EdgeBuffer_t) , intent(inout) :: edgeq
253 : type (derivative_t) , intent(in) :: deriv
254 :
255 : ! local
256 : integer :: k,kptr,i,j,ie,ic,q
257 : integer :: kbeg,kend,qbeg,qend
258 : real (kind=r8), dimension(np,np) :: lap_p
259 : logical var_coef1
260 : integer :: kblk,qblk ! The per thead size of the vertical and tracers
261 :
262 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
263 :
264 : !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
265 : !so tensor is only used on second call to laplace_sphere_wk
266 0 : var_coef1 = .true.
267 0 : if(hypervis_scaling > 0) var_coef1 = .false.
268 :
269 :
270 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
271 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
272 :
273 0 : do ie=nets,nete
274 0 : do q=qbeg,qend
275 0 : do k=kbeg,kend
276 0 : lap_p(:,:)=qtens(:,:,k,q,ie)
277 0 : call laplace_sphere_wk(lap_p,deriv,elem(ie),qtens(:,:,k,q,ie),var_coef=var_coef1)
278 : enddo
279 0 : kptr = nlev*(q-1) + kbeg - 1
280 0 : call edgeVpack(edgeq, qtens(:,:,kbeg:kend,q,ie),kblk,kptr,ie)
281 : enddo
282 : enddo
283 :
284 :
285 0 : call bndry_exchange(hybrid,edgeq,location='biharmonic_wk_scalar')
286 :
287 0 : do ie=nets,nete
288 :
289 : ! apply inverse mass matrix, then apply laplace again
290 0 : do q=qbeg,qend
291 0 : kptr = nlev*(q-1) + kbeg - 1
292 0 : call edgeVunpack(edgeq, qtens(:,:,kbeg:kend,q,ie),kblk,kptr,ie)
293 0 : do k=kbeg,kend
294 0 : lap_p(:,:)=elem(ie)%rspheremp(:,:)*qtens(:,:,k,q,ie)
295 0 : call laplace_sphere_wk(lap_p,deriv,elem(ie),qtens(:,:,k,q,ie),var_coef=.true.)
296 : enddo
297 : enddo
298 : enddo
299 :
300 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
301 0 : end subroutine biharmonic_wk_scalar
302 :
303 :
304 0 : subroutine make_C0(zeta,elem,hybrid,nets,nete)
305 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306 : ! apply DSS (aka assembly procedure) to zeta.
307 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308 :
309 : type (hybrid_t) , intent(in) :: hybrid
310 : type (element_t) , intent(in), target :: elem(:)
311 : integer :: nets,nete
312 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
313 :
314 : ! local
315 : integer :: k,i,j,ie,ic,kptr,nthread_save
316 :
317 :
318 0 : call initEdgeBuffer(hybrid%par,edge1,elem,nlev)
319 :
320 0 : do ie=nets,nete
321 : #if (defined COLUMN_OPENMP)
322 : !$omp parallel do num_threads(vert_num_threads) private(k)
323 : #endif
324 0 : do k=1,nlev
325 0 : zeta(:,:,k,ie)=zeta(:,:,k,ie)*elem(ie)%spheremp(:,:)
326 : enddo
327 0 : kptr=0
328 0 : call edgeVpack(edge1, zeta(1,1,1,ie),nlev,kptr,ie)
329 : enddo
330 0 : call bndry_exchange(hybrid,edge1,location='make_C0')
331 0 : do ie=nets,nete
332 0 : kptr=0
333 0 : call edgeVunpack(edge1, zeta(1,1,1,ie),nlev,kptr, ie)
334 : #if (defined COLUMN_OPENMP)
335 : !$omp parallel do num_threads(vert_num_threads) private(k)
336 : #endif
337 0 : do k=1,nlev
338 0 : zeta(:,:,k,ie)=zeta(:,:,k,ie)*elem(ie)%rspheremp(:,:)
339 : enddo
340 : enddo
341 :
342 0 : call FreeEdgeBuffer(edge1)
343 :
344 0 : end subroutine
345 :
346 :
347 0 : subroutine make_C0_vector(v,elem,hybrid,nets,nete)
348 : #if 1
349 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350 : ! apply DSS to a velocity vector
351 : ! this is a low-performance routine used for I/O and analysis.
352 : ! no need to optimize
353 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
354 : type (hybrid_t) , intent(in) :: hybrid
355 : type (element_t) , intent(in), target :: elem(:)
356 : integer :: nets,nete
357 : real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: v
358 :
359 : ! local
360 : integer :: k,i,j,ie,ic,kptr
361 : type (EdgeBuffer_t) :: edge2
362 0 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: v1
363 :
364 0 : v1(:,:,:,:) = v(:,:,1,:,:)
365 0 : call make_C0(v1,elem,hybrid,nets,nete)
366 0 : v(:,:,1,:,:) = v1(:,:,:,:)
367 :
368 0 : v1(:,:,:,:) = v(:,:,2,:,:)
369 0 : call make_C0(v1,elem,hybrid,nets,nete)
370 0 : v(:,:,2,:,:) = v1(:,:,:,:)
371 : #else
372 : type (hybrid_t) , intent(in) :: hybrid
373 : type (element_t) , intent(in), target :: elem(:)
374 : integer :: nets,nete
375 : real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: v
376 :
377 : ! local
378 : integer :: k,i,j,ie,ic,kptr
379 : type (EdgeBuffer_t) :: edge2
380 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: v1
381 :
382 :
383 :
384 : call initEdgeBuffer(hybrid%par,edge2,elem,2*nlev)
385 :
386 : do ie=nets,nete
387 : #if (defined COLUMN_OPENMP)
388 : !$omp parallel do num_threads(vert_num_threads) private(k)
389 : #endif
390 : do k=1,nlev
391 : v(:,:,1,k,ie)=v(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
392 : v(:,:,2,k,ie)=v(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
393 : enddo
394 : kptr=0
395 : call edgeVpack(edge2, v(1,1,1,1,ie),2*nlev,kptr,ie)
396 : enddo
397 : call bndry_exchange(hybrid,edge2,location='make_C0_vector')
398 : do ie=nets,nete
399 : kptr=0
400 : call edgeVunpack(edge2, v(1,1,1,1,ie),2*nlev,kptr,ie)
401 : #if (defined COLUMN_OPENMP)
402 : !$omp parallel do num_threads(vert_num_threads) private(k)
403 : #endif
404 : do k=1,nlev
405 : v(:,:,1,k,ie)=v(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
406 : v(:,:,2,k,ie)=v(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
407 : enddo
408 : enddo
409 :
410 : call FreeEdgeBuffer(edge2)
411 : #endif
412 0 : end subroutine
413 :
414 :
415 :
416 :
417 :
418 :
419 0 : subroutine compute_zeta_C0_contra(zeta,elem,hybrid,nets,nete,nt)
420 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
421 : ! compute C0 vorticity. That is, solve:
422 : ! < PHI, zeta > = <PHI, curl(elem%state%v >
423 : !
424 : ! input: v (stored in elem()%, in contra-variant coordinates)
425 : ! output: zeta(:,:,:,:)
426 : !
427 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428 :
429 : type (hybrid_t) , intent(in) :: hybrid
430 : type (element_t) , intent(in), target :: elem(:)
431 : integer :: nt,nets,nete
432 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
433 : real (kind=r8), dimension(np,np,2) :: ulatlon
434 : real (kind=r8), dimension(np,np) :: v1,v2
435 :
436 : ! local
437 : integer :: k,ie
438 : type (derivative_t) :: deriv
439 :
440 0 : call derivinit(deriv)
441 :
442 0 : do k=1,nlev
443 0 : do ie=nets,nete
444 0 : v1 = elem(ie)%state%v(:,:,1,k,nt)
445 0 : v2 = elem(ie)%state%v(:,:,2,k,nt)
446 0 : ulatlon(:,:,1) = elem(ie)%D(:,:,1,1)*v1 + elem(ie)%D(:,:,1,2)*v2
447 0 : ulatlon(:,:,2) = elem(ie)%D(:,:,2,1)*v1 + elem(ie)%D(:,:,2,2)*v2
448 0 : call vorticity_sphere(ulatlon,deriv,elem(ie),zeta(:,:,k,ie))
449 : enddo
450 : enddo
451 :
452 0 : call make_C0(zeta,elem,hybrid,nets,nete)
453 :
454 0 : end subroutine
455 :
456 :
457 :
458 0 : subroutine compute_div_C0_contra(zeta,elem,hybrid,nets,nete,nt)
459 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
460 : ! compute C0 divergence. That is, solve:
461 : ! < PHI, zeta > = <PHI, div(elem%state%v >
462 : !
463 : ! input: v (stored in elem()%, in contra-variant coordinates)
464 : ! output: zeta(:,:,:,:)
465 : !
466 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
467 :
468 : type (hybrid_t) , intent(in) :: hybrid
469 : type (element_t) , intent(in), target :: elem(:)
470 : integer :: nt,nets,nete
471 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
472 : real (kind=r8), dimension(np,np,2) :: ulatlon
473 : real (kind=r8), dimension(np,np) :: v1,v2
474 :
475 : ! local
476 : integer :: k,ie
477 : type (derivative_t) :: deriv
478 :
479 0 : call derivinit(deriv)
480 :
481 0 : do k=1,nlev
482 0 : do ie=nets,nete
483 0 : v1 = elem(ie)%state%v(:,:,1,k,nt)
484 0 : v2 = elem(ie)%state%v(:,:,2,k,nt)
485 0 : ulatlon(:,:,1) = elem(ie)%D(:,:,1,1)*v1 + elem(ie)%D(:,:,1,2)*v2
486 0 : ulatlon(:,:,2) = elem(ie)%D(:,:,2,1)*v1 + elem(ie)%D(:,:,2,2)*v2
487 0 : call divergence_sphere(ulatlon,deriv,elem(ie),zeta(:,:,k,ie))
488 : enddo
489 : enddo
490 :
491 0 : call make_C0(zeta,elem,hybrid,nets,nete)
492 :
493 0 : end subroutine
494 :
495 0 : subroutine compute_zeta_C0_par(zeta,elem,par,nt)
496 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
497 : ! compute C0 vorticity. That is, solve:
498 : ! < PHI, zeta > = <PHI, curl(elem%state%v >
499 : !
500 : ! input: v (stored in elem()%, in lat-lon coordinates)
501 : ! output: zeta(:,:,:,:)
502 : !
503 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
504 : type (parallel_t) :: par
505 : type (element_t) , intent(in), target :: elem(:)
506 : real (kind=r8), dimension(np,np,nlev,nelemd) :: zeta
507 : integer :: nt
508 :
509 : ! local
510 : type (hybrid_t) :: hybrid
511 : integer :: k,i,j,ie,ic
512 : type (derivative_t) :: deriv
513 :
514 : ! single thread
515 0 : hybrid = config_thread_region(par,'serial')
516 :
517 0 : call compute_zeta_C0_hybrid(zeta,elem,hybrid,1,nelemd,nt)
518 :
519 0 : end subroutine
520 :
521 :
522 0 : subroutine compute_div_C0_par(zeta,elem,par,nt)
523 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
524 : ! compute C0 divergence. That is, solve:
525 : ! < PHI, zeta > = <PHI, div(elem%state%v >
526 : !
527 : ! input: v (stored in elem()%, in lat-lon coordinates)
528 : ! output: zeta(:,:,:,:)
529 : !
530 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
531 :
532 : type (parallel_t) :: par
533 : type (element_t) , intent(in), target :: elem(:)
534 : real (kind=r8), dimension(np,np,nlev,nelemd) :: zeta
535 : integer :: nt
536 :
537 : ! local
538 : type (hybrid_t) :: hybrid
539 : integer :: k,i,j,ie,ic
540 : type (derivative_t) :: deriv
541 :
542 : ! single thread
543 0 : hybrid = config_thread_region(par,'serial')
544 :
545 0 : call compute_div_C0_hybrid(zeta,elem,hybrid,1,nelemd,nt)
546 :
547 0 : end subroutine
548 :
549 :
550 :
551 0 : subroutine compute_zeta_C0_hybrid(zeta,elem,hybrid,nets,nete,nt)
552 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
553 : ! compute C0 vorticity. That is, solve:
554 : ! < PHI, zeta > = <PHI, curl(elem%state%v >
555 : !
556 : ! input: v (stored in elem()%, in lat-lon coordinates)
557 : ! output: zeta(:,:,:,:)
558 : !
559 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
560 :
561 : type (hybrid_t) , intent(in) :: hybrid
562 : type (element_t) , intent(in), target :: elem(:)
563 : integer :: nt,nets,nete
564 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
565 :
566 : ! local
567 : integer :: k,i,j,ie,ic
568 : type (derivative_t) :: deriv
569 :
570 0 : call derivinit(deriv)
571 :
572 0 : do ie=nets,nete
573 : #if (defined COLUMN_OPENMP)
574 : !$omp parallel do num_threads(vert_num_threads) private(k)
575 : #endif
576 0 : do k=1,nlev
577 0 : call vorticity_sphere(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),zeta(:,:,k,ie))
578 : enddo
579 : enddo
580 :
581 0 : call make_C0(zeta,elem,hybrid,nets,nete)
582 :
583 0 : end subroutine
584 :
585 :
586 0 : subroutine compute_div_C0_hybrid(zeta,elem,hybrid,nets,nete,nt)
587 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
588 : ! compute C0 divergence. That is, solve:
589 : ! < PHI, zeta > = <PHI, div(elem%state%v >
590 : !
591 : ! input: v (stored in elem()%, in lat-lon coordinates)
592 : ! output: zeta(:,:,:,:)
593 : !
594 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
595 :
596 : type (hybrid_t) , intent(in) :: hybrid
597 : type (element_t) , intent(in), target :: elem(:)
598 : integer :: nt,nets,nete
599 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
600 :
601 : ! local
602 : integer :: k,i,j,ie,ic
603 : type (derivative_t) :: deriv
604 :
605 0 : call derivinit(deriv)
606 :
607 0 : do ie=nets,nete
608 : #if (defined COLUMN_OPENMP)
609 : !$omp parallel do num_threads(vert_num_threads) private(k)
610 : #endif
611 0 : do k=1,nlev
612 0 : call divergence_sphere(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),zeta(:,:,k,ie))
613 : enddo
614 : enddo
615 :
616 0 : call make_C0(zeta,elem,hybrid,nets,nete)
617 :
618 0 : end subroutine
619 :
620 :
621 :
622 :
623 :
624 :
625 :
626 :
627 0 : subroutine neighbor_minmax(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
628 :
629 : type (hybrid_t) , intent(in) :: hybrid
630 : type (EdgeBuffer_t) , intent(inout) :: edgeMinMax
631 : integer :: nets,nete
632 : real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
633 : real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
634 : integer :: kblk, qblk
635 : ! local
636 : integer:: ie, q, k, kptr
637 : integer:: kbeg, kend, qbeg, qend
638 :
639 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
640 :
641 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
642 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
643 :
644 0 : do ie=nets,nete
645 0 : do q = qbeg, qend
646 0 : kptr = nlev*(q - 1) + kbeg - 1
647 0 : call edgeSpack(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
648 0 : kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
649 0 : call edgeSpack(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
650 : enddo
651 : enddo
652 :
653 0 : call bndry_exchange(hybrid,edgeMinMax,location='neighbor_minmax')
654 :
655 0 : do ie=nets,nete
656 0 : do q=qbeg,qend
657 0 : kptr = nlev*(q - 1) + kbeg - 1
658 0 : call edgeSunpackMIN(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
659 0 : kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
660 0 : call edgeSunpackMAX(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
661 0 : do k=kbeg,kend
662 0 : min_neigh(k,q,ie) = max(min_neigh(k,q,ie),0.0_r8)
663 : enddo
664 : enddo
665 : enddo
666 :
667 0 : end subroutine neighbor_minmax
668 :
669 :
670 0 : subroutine neighbor_minmax_start(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
671 :
672 : type (hybrid_t) , intent(in) :: hybrid
673 : type (EdgeBuffer_t) , intent(inout) :: edgeMinMax
674 : integer :: nets,nete
675 : real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
676 : real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
677 : integer :: kblk, qblk
678 : integer :: kbeg, kend, qbeg, qend
679 :
680 : ! local
681 : integer :: ie,q, k,kptr
682 :
683 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
684 :
685 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
686 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
687 :
688 0 : do ie=nets,nete
689 0 : do q=qbeg, qend
690 0 : kptr = nlev*(q - 1) + kbeg - 1
691 0 : call edgeSpack(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
692 0 : kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
693 0 : call edgeSpack(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
694 : enddo
695 : enddo
696 :
697 0 : call bndry_exchange_start(hybrid,edgeMinMax,location='neighbor_minmax_start')
698 :
699 0 : end subroutine neighbor_minmax_start
700 :
701 0 : subroutine neighbor_minmax_finish(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
702 :
703 : type (hybrid_t) , intent(in) :: hybrid
704 : type (EdgeBuffer_t) , intent(inout) :: edgeMinMax
705 : integer :: nets,nete
706 : real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
707 : real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
708 : integer :: kblk, qblk
709 : integer :: ie,q, k,kptr
710 : integer :: kbeg, kend, qbeg, qend
711 :
712 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
713 :
714 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
715 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
716 :
717 0 : call bndry_exchange_finish(hybrid,edgeMinMax,location='neighbor_minmax_finish')
718 :
719 0 : do ie=nets,nete
720 0 : do q=qbeg, qend
721 0 : kptr = nlev*(q - 1) + kbeg - 1
722 0 : call edgeSunpackMIN(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
723 0 : kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
724 0 : call edgeSunpackMAX(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
725 0 : do k=kbeg,kend
726 0 : min_neigh(k,q,ie) = max(min_neigh(k,q,ie),0.0_r8)
727 : enddo
728 : enddo
729 : enddo
730 :
731 0 : end subroutine neighbor_minmax_finish
732 :
733 : end module viscosity_mod
|
