1 |
subroutine plume2dyn(qplume,Nxplume,Lmplume,uref,vref,flag, |
2 |
. idim1,idim2,jdim1,jdim2,i1,i2,j1,j2,Nsx,Nsy,bi,bj,qdyn1,qdyn2) |
3 |
C*********************************************************************** |
4 |
C Purpose: |
5 |
C To interpolate an arbitrary quantity from higher resolution plume |
6 |
C grid to the model's dynamics grid |
7 |
C Algorithm: |
8 |
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 |
C |
12 |
C Input: |
13 |
C qplume... [idim2,jdim2,im,Lmplume,bi] Quantity on Input Grid |
14 |
C Nxplume . Longitude Dimension of Input |
15 |
C Lmplume.. Vertical Dimension of Input |
16 |
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 |
C flag .... Flag to indicate vector (1) or scalar (0) interpolation |
19 |
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 |
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 |
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 |
C |
32 |
C Notes: |
33 |
C 1) Assume (for now) that the number of vertical levels is the |
34 |
C same on both the input and output grids |
35 |
C*********************************************************************** |
36 |
implicit none |
37 |
#include "CPP_OPTIONS.h" |
38 |
|
39 |
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 |
|
48 |
integer i,j,L,iplume |
49 |
_RL qplumeav(i2,j2,Lmplume) |
50 |
_RL sqrtarg |
51 |
|
52 |
C First step - compute the average of qplume over Nxplume |
53 |
do j = j1,j2 |
54 |
do i = i1,i2 |
55 |
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 |
|
64 |
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 |
|
69 |
if (flag.eq.0) then |
70 |
do j = j1,j2 |
71 |
do i = i1,i2 |
72 |
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 |
do j = j1,j2 |
79 |
do i = i1,i2 |
80 |
do L = 1,Lmplume |
81 |
if(vref(i,j,L,bi,bj).ne.0.) then |
82 |
sqrtarg = (qplumeav(i,j,L)*qplumeav(i,j,L)) / |
83 |
. ( ( (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 |
qdyn2(i,j,L,bi,bj) = sqrt(sqrtarg) |
86 |
qdyn1(i,j,L,bi,bj) = qdyn2(i,j,L,bi,bj) * |
87 |
. (uref(i,j,L,bi,bj)/vref(i,j,L,bi,bj)) |
88 |
else |
89 |
qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L) |
90 |
qdyn2(i,j,L,bi,bj) = 0. |
91 |
endif |
92 |
enddo |
93 |
enddo |
94 |
enddo |
95 |
endif |
96 |
|
97 |
return |
98 |
end |