/[MITgcm]/MITgcm_contrib/plumes/plume2dyn.F
ViewVC logotype

Diff of /MITgcm_contrib/plumes/plume2dyn.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by molod, Thu May 13 22:21:45 2004 UTC revision 1.6 by molod, Tue Jun 1 22:31:34 2004 UTC
# Line 1  Line 1 
1        subroutine plume2dyn(qplume,idimin,jdimin,Lmplume,        subroutine plume2dyn(qplume,Nxplume,Lmplume,uref,vref,flag,
2       .     idim1,idim2,jdim1,jdim2,Lmout,Nsx,Nsy,bi,bj,qdyn)       . idim1,idim2,jdim1,jdim2,i1,i2,j1,j2,Nsx,Nsy,bi,bj,qdyn1,qdyn2)
3  C***********************************************************************  C***********************************************************************
4  C Purpose:  C Purpose:
5  C   To interpolate an arbitrary quantity from higher resolution plumes  C   To interpolate an arbitrary quantity from higher resolution plume
6  C         grid to the model's dynamics grid  C         grid to the model's dynamics grid
7  C Algorithm:  C Algorithm:
8  C   Plumes -> Dynamics computes the plumes are mean value  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  C
12  C Input:  C Input:
13  C   qplume... [im,jm,Lmplume] Arbitrary Quantity on Input Grid  C   qplume... [idim2,jdim2,im,Lmplume,bi] Quantity on Input Grid
14  C   pephy.... [im,jm,Lmplume+1] Pressures at bottom edges of input levels  C   Nxplume . Longitude Dimension of Input
 C   idimin... Longitude Dimension of Input  
 C   jdimin... Latitude  Dimension of Input  
15  C   Lmplume.. Vertical  Dimension of Input  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  C   Nsx...... Number of processes in x-direction
24  C   Nsy...... Number of processes in y-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  
25  C   bi....... Index of process number in x-direction  C   bi....... Index of process number in x-direction
26  C   bj....... 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  
27  C  C
28  C Output:  C Output:
29  C   qdyn..... [im,jm,Lmout] Quantity at output grid (physics grid)  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  C
32  C Notes:  C Notes:
33  C   1) This algorithm assumes that the output (physics) grid levels  C   1)  Assume (for now) that the number of vertical levels is the
34  C      fit exactly into the input (dynamics) grid levels  C       same on both the input and output grids
35  C***********************************************************************  C***********************************************************************
36        implicit none        implicit none
37  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
38    
39        integer  idimin, jdimin, Lmout, Lmplume, Nsx, Nsy        integer  Nxplume, Lmplume, Nsx, Nsy
40        integer idim1, idim2, jdim1, jdim2, bi, bj        integer idim1, idim2, jdim1, jdim2, i1, i2, j1, j2
41        _RL qplume(idimin,jdimin,Lmplume,Nsx,Nsy)        integer bi, bj, flag
42        _RL pedyn(idimin,jdimin,Lmout+1,Nsx,Nsy)        _RL qplume(i2,j2,Nxplume,Lmplume,Nsx)
43        _RL pephy(idimin,jdimin,Lmplume+1,Nsx,Nsy)        _RL uref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
44        integer nlperdyn(idimin,jdimin,Lmout,Nsx,Nsy)        _RL vref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
45        _RL qdyn(idimin,jdimin,Lmout,Nsx,Nsy)        _RL qdyn1(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
46        integer Lbot(idimin,jdimin,Nsx,Nsy)        _RL qdyn2(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
47    
48        integer  i,j,L,Lout1,Lout1p1,Lout2,Lphy        integer i,j,L,iplume
49        _RL getcon, kappa, dpkephy, dpkedyn, sum        _RL qplumeav(i2,j2,Lmplume)
50          _RL sqrtarg
51        kappa = getcon('KAPPA')  
52    C First step - compute the average of qplume over Nxplume
53  c do loop for all dynamics (output) levels        do j = j1,j2
54        do L = 1,Lmout        do i = i1,i2
55  c do loop for all grid points         do L = 1,Lmplume
56         do j = jdim1,jdim2          qplumeav(i,j,L) = 0.
57          do i = idim1,idim2          do iplume = 1,Nxplume
58           qdyn(i,j,L,bi,bj) = 0.           qplumeav(i,j,L)=qplumeav(i,j,L)+qplume(i,j,iplume,L,bi)/Nxplume
 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  
59          enddo          enddo
60         enddo         enddo
61        enddo        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        return
98        end        end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22