1 |
molod |
1.2 |
subroutine plume2dyn(qplume,Nxplume,Lmplume,uref,vref,flag, |
2 |
molod |
1.4 |
. idim1,idim2,jdim1,jdim2,i1,i2,j1,j2,Nsx,Nsy,bi,bj,qdyn1,qdyn2) |
3 |
molod |
1.1 |
C*********************************************************************** |
4 |
|
|
C Purpose: |
5 |
molod |
1.2 |
C To interpolate an arbitrary quantity from higher resolution plume |
6 |
molod |
1.1 |
C grid to the model's dynamics grid |
7 |
|
|
C Algorithm: |
8 |
molod |
1.2 |
C Plumes -> Dynamics computes the plumes mean value, and in the case |
9 |
|
|
C of a vector field, preserves the direction of a vector |
10 |
|
|
C given in (uref,vref) |
11 |
molod |
1.1 |
C |
12 |
|
|
C Input: |
13 |
molod |
1.2 |
C qplume... [idim2,jdim2,im,Lmplume,bi] Quantity on Input Grid |
14 |
molod |
1.4 |
C Nxplume . Longitude Dimension of Input |
15 |
molod |
1.1 |
C Lmplume.. Vertical Dimension of Input |
16 |
molod |
1.4 |
C uref .... [im,jm,Lmplume,bi,bj] Reference u-component of velocity |
17 |
|
|
C vref .... [im,jm,plume,bi,bj] Reference v-component of velocity |
18 |
molod |
1.2 |
C flag .... Flag to indicate vector (1) or scalar (0) interpolation |
19 |
molod |
1.4 |
C idim1,2.. Beginning and ending dimension of output grid |
20 |
|
|
C jdim1,2.. Beginning and ending dimension of output grid |
21 |
|
|
C i1,2..... Beginning and ending x-direction span |
22 |
|
|
C j1,2..... Beginning and ending y-direction span |
23 |
molod |
1.1 |
C Nsx...... Number of processes in x-direction |
24 |
|
|
C Nsy...... Number of processes in y-direction |
25 |
|
|
C bi....... Index of process number in x-direction |
26 |
|
|
C bj....... Index of process number in x-direction |
27 |
|
|
C |
28 |
|
|
C Output: |
29 |
molod |
1.4 |
C qdyn1..... [im,jm,plume,bi,bj] Field at output grid (dynamics) |
30 |
|
|
C qdyn2..... [im,jm,plume,bi,bj] Field at output grid (dynamics) |
31 |
molod |
1.1 |
C |
32 |
|
|
C Notes: |
33 |
molod |
1.2 |
C 1) Assume (for now) that the number of vertical levels is the |
34 |
|
|
C same on both the input and output grids |
35 |
molod |
1.1 |
C*********************************************************************** |
36 |
|
|
implicit none |
37 |
|
|
#include "CPP_OPTIONS.h" |
38 |
|
|
|
39 |
molod |
1.4 |
integer Nxplume, Lmplume, Nsx, Nsy |
40 |
|
|
integer idim1, idim2, jdim1, jdim2, i1, i2, j1, j2 |
41 |
|
|
integer bi, bj, flag |
42 |
|
|
_RL qplume(i2,j2,Nxplume,Lmplume,Nsx) |
43 |
|
|
_RL uref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy) |
44 |
|
|
_RL vref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy) |
45 |
|
|
_RL qdyn1(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy) |
46 |
|
|
_RL qdyn2(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy) |
47 |
molod |
1.1 |
|
48 |
molod |
1.2 |
integer i,j,L,iplume |
49 |
molod |
1.4 |
_RL qplumeav(i1,j2,Lmplume) |
50 |
molod |
1.2 |
_RL sqrtarg |
51 |
|
|
|
52 |
|
|
C First step - compute the average of qplume over Nxplume |
53 |
molod |
1.4 |
do j = j1,j2 |
54 |
|
|
do i = i1,i2 |
55 |
molod |
1.2 |
do L = 1,Lmplume |
56 |
|
|
qplumeav(i,j,L) = 0. |
57 |
|
|
do iplume = 1,Nxplume |
58 |
|
|
qplumeav(i,j,L)=qplumeav(i,j,L)+qplume(i,j,iplume,L,bi)/Nxplume |
59 |
|
|
enddo |
60 |
|
|
enddo |
61 |
|
|
enddo |
62 |
|
|
enddo |
63 |
molod |
1.1 |
|
64 |
molod |
1.2 |
C Now check the flag -- if a scalar, we are done - just assign |
65 |
|
|
C the average to all the i and j points of the output grid. |
66 |
|
|
C If a vector, there is some more work to do in order to preserve |
67 |
|
|
C the angle given by uref and vref |
68 |
molod |
1.1 |
|
69 |
molod |
1.2 |
if (flag.eq.0) then |
70 |
molod |
1.4 |
do j = j1,j2 |
71 |
|
|
do i = i1,i2 |
72 |
molod |
1.2 |
do L = 1,Lmplume |
73 |
|
|
qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L) |
74 |
|
|
enddo |
75 |
|
|
enddo |
76 |
|
|
enddo |
77 |
|
|
elseif (flag.eq.1) then |
78 |
molod |
1.4 |
do j = j1,j2 |
79 |
|
|
do i = i1,i2 |
80 |
molod |
1.2 |
do L = 1,Lmplume |
81 |
molod |
1.3 |
if(vref(i,j,L,bi,bj).ne.0.) then |
82 |
|
|
sqrtarg = (qplumeav(i,j,L)*qplumeav(i,j,L)) / |
83 |
molod |
1.2 |
. ( ( (uref(i,j,L,bi,bj)*uref(i,j,L,bi,bj)) / |
84 |
|
|
. (vref(i,j,L,bi,bj)*vref(i,j,L,bi,bj)) ) + 1. ) |
85 |
molod |
1.3 |
qdyn2(i,j,L,bi,bj) = sqrt(sqrtarg) |
86 |
|
|
qdyn1(i,j,L,bi,bj) = qdyn2(i,j,L,bi,bj) * |
87 |
molod |
1.2 |
. (uref(i,j,L,bi,bj)/vref(i,j,L,bi,bj)) |
88 |
molod |
1.3 |
else |
89 |
|
|
qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L) |
90 |
|
|
qdyn2(i,j,L,bi,bj) = 0. |
91 |
molod |
1.2 |
enddo |
92 |
|
|
enddo |
93 |
molod |
1.1 |
enddo |
94 |
molod |
1.2 |
endif |
95 |
molod |
1.1 |
|
96 |
|
|
return |
97 |
|
|
end |