Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch
--- MITgcm_contrib/plumes/plume2dyn.F 2004/05/13 22:21:45 1.1
+++ MITgcm_contrib/plumes/plume2dyn.F 2004/05/25 18:08:51 1.2
@@ -1,83 +1,91 @@
- subroutine plume2dyn(qplume,idimin,jdimin,Lmplume,
- . idim1,idim2,jdim1,jdim2,Lmout,Nsx,Nsy,bi,bj,qdyn)
+ subroutine plume2dyn(qplume,Nxplume,Lmplume,uref,vref,flag,
+ . idim1,idim2,jdim1,jdim2,Lmout,Nsx,Nsy,bi,bj,qdyn1,qdyn2)
C***********************************************************************
C Purpose:
-C To interpolate an arbitrary quantity from higher resolution plumes
+C To interpolate an arbitrary quantity from higher resolution plume
C grid to the model's dynamics grid
C Algorithm:
-C Plumes -> Dynamics computes the plumes are mean value
+C Plumes -> Dynamics computes the plumes mean value, and in the case
+C of a vector field, preserves the direction of a vector
+C given in (uref,vref)
C
C Input:
-C qplume... [im,jm,Lmplume] Arbitrary Quantity on Input Grid
-C pephy.... [im,jm,Lmplume+1] Pressures at bottom edges of input levels
+C qplume... [idim2,jdim2,im,Lmplume,bi] Quantity on Input Grid
C idimin... Longitude Dimension of Input
-C jdimin... Latitude Dimension of Input
C Lmplume.. Vertical Dimension of Input
+C uref .... [im,jm,Lmout,bi,bj] Reference u-component of velocity
+C vref .... [im,jm,Lmout,bi,bj] Reference v-component of velocity
+C flag .... Flag to indicate vector (1) or scalar (0) interpolation
+C idim1,2.. Beginning and ending i-values of output grid
+C jdim1,2.. Beginning and ending j-values of output grid
+C Lmout.... Vertical Dimension of Output
C Nsx...... Number of processes in x-direction
C Nsy...... Number of processes in y-direction
-C idim1,2.. Beginning and ending i-values to calculate
-C jdim1,2.. Beginning and ending j-values to calculate
C bi....... Index of process number in x-direction
C bj....... Index of process number in x-direction
-C pedyn.... [im,jm,Lmout+1] Pressures at bottom edges of output levels
-C Lmout.... Vertical Dimension of Output
-C nlperdyn. Mapping Array-Highest Physics level in each dynmics level
C
C Output:
-C qdyn..... [im,jm,Lmout] Quantity at output grid (physics grid)
+C qdyn1..... [im,jm,Lmout,bi,bj] Field at output grid (dynamics)
+C qdyn2..... [im,jm,Lmout,bi,bj] Field at output grid (dynamics)
C
C Notes:
-C 1) This algorithm assumes that the output (physics) grid levels
-C fit exactly into the input (dynamics) grid levels
+C 1) Assume (for now) that the number of vertical levels is the
+C same on both the input and output grids
C***********************************************************************
implicit none
#include "CPP_OPTIONS.h"
- integer idimin, jdimin, Lmout, Lmplume, Nsx, Nsy
- integer idim1, idim2, jdim1, jdim2, bi, bj
- _RL qplume(idimin,jdimin,Lmplume,Nsx,Nsy)
- _RL pedyn(idimin,jdimin,Lmout+1,Nsx,Nsy)
- _RL pephy(idimin,jdimin,Lmplume+1,Nsx,Nsy)
- integer nlperdyn(idimin,jdimin,Lmout,Nsx,Nsy)
- _RL qdyn(idimin,jdimin,Lmout,Nsx,Nsy)
- integer Lbot(idimin,jdimin,Nsx,Nsy)
-
- integer i,j,L,Lout1,Lout1p1,Lout2,Lphy
- _RL getcon, kappa, dpkephy, dpkedyn, sum
-
- kappa = getcon('KAPPA')
-
-c do loop for all dynamics (output) levels
- do L = 1,Lmout
-c do loop for all grid points
- do j = jdim1,jdim2
- do i = idim1,idim2
- qdyn(i,j,L,bi,bj) = 0.
-c Check to make sure we are above ground - otherwise do nothing
- if(L.ge.Lbot(i,j,bi,bj))then
- if(L.eq.Lbot(i,j,bi,bj)) then
- Lout1 = 0
- else
- Lout1 = nlperdyn(i,j,L-1,bi,bj)
- endif
- Lout2 = nlperdyn(i,j,L,bi,bj)
-c do loop for all physics levels contained in this dynamics level
-cinterp1 dpkedyn = (pedyn(i,j,L,bi,bj)**kappa)-
-cinterp1 (pedyn(i,j,L+1,bi,bj)**kappa)
- dpkedyn = pedyn(i,j,L,bi,bj)-pedyn(i,j,L+1,bi,bj)
- sum = 0.
- Lout1p1 = Lout1+1
- do Lphy = Lout1p1,Lout2
-cinterp1 dpkephy = (pephy(i,j,Lphy,bi,bj)**kappa)-
-cinterp1 (pephy(i,j,Lphy+1,bi,bj)**kappa)
- dpkephy = pephy(i,j,Lphy,bi,bj)-pephy(i,j,Lphy+1,bi,bj)
- sum=sum+qplume(i,j,Lphy,bi,bj)*(dpkephy/dpkedyn)
- enddo
- qdyn(i,j,L,bi,bj) = sum
- endif
+ integer Nxplume, Lmplume, Lmout, Nsx, Nsy
+ integer idim1, idim2, jdim1, jdim2, bi, bj, flag
+ _RL qplume(idim2,jdim2,Nxplume,Lmplume,Nsx)
+ _RL uref(idim1:idim2,jdim1:jdim2,Lmout,Nsx,Nsy)
+ _RL vref(idim1:idim2,jdim1:jdim2,Lmout,Nsx,Nsy)
+ _RL qdyn1(idim1:idim2,jdim1:jdim2,Lmout,Nsx,Nsy)
+ _RL qdyn2(idim1:idim2,jdim1:jdim2,Lmout,Nsx,Nsy)
+
+ integer i,j,L,iplume
+ _RL qplumeav(idim1,jdim2,Lmplume)
+ _RL sqrtarg
+
+C First step - compute the average of qplume over Nxplume
+ do j = jdim1,jdim2
+ do i = idim1,idim2
+ do L = 1,Lmplume
+ qplumeav(i,j,L) = 0.
+ do iplume = 1,Nxplume
+ qplumeav(i,j,L)=qplumeav(i,j,L)+qplume(i,j,iplume,L,bi)/Nxplume
enddo
enddo
enddo
+ enddo
+
+C Now check the flag -- if a scalar, we are done - just assign
+C the average to all the i and j points of the output grid.
+C If a vector, there is some more work to do in order to preserve
+C the angle given by uref and vref
+
+ if (flag.eq.0) then
+ do j = jdim1,jdim2
+ do i = idim1,idim2
+ do L = 1,Lmplume
+ qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L)
+ enddo
+ enddo
+ enddo
+ elseif (flag.eq.1) then
+ do j = jdim1,jdim2
+ do i = idim1,idim2
+ do L = 1,Lmplume
+ sqrtarg = (qplumeav(i,j,L)*qplumeav(i,j,L)) /
+ . ( ( (uref(i,j,L,bi,bj)*uref(i,j,L,bi,bj)) /
+ . (vref(i,j,L,bi,bj)*vref(i,j,L,bi,bj)) ) + 1. )
+ qdyn2(i,j,L,bi,bj) = sqrt(sqrtarg)
+ qdyn1(i,j,L,bi,bj) = qdyn2(i,j,L,bi,bj) *
+ . (uref(i,j,L,bi,bj)/vref(i,j,L,bi,bj))
+ enddo
+ enddo
+ enddo
+ endif
return
end
| ViewVC Help | |
| Powered by ViewVC 1.1.22 |