--- 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