/[MITgcm]/MITgcm/pkg/diagnostics/diagnostics_interp_vert.F
ViewVC logotype

Annotation of /MITgcm/pkg/diagnostics/diagnostics_interp_vert.F

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


Revision 1.2 - (hide annotations) (download)
Tue Oct 25 16:11:20 2005 UTC (18 years, 7 months ago) by molod
Branch: MAIN
Changes since 1.1: +39 -14 lines
Add capability for vertical interpolation on the fly to APE levels

1 molod 1.2 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_interp_vert.F,v 1.1 2005/08/25 21:54:54 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "DIAG_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7     CBOP 0
8     C !ROUTINE: DIAGNOSTICS_INTERP_VERT
9    
10     C !INTERFACE:
11     SUBROUTINE DIAGNOSTICS_INTERP_VERT(
12     I listId, md, ndId, ip, im,
13     U nlevsout,
14     U qtmp1,
15     I undef,
16     I myTime, myIter, myThid )
17    
18     C !DESCRIPTION:
19     C Interpolate vertically a diagnostics field before writing to file.
20     C presently implemented (for Atmospheric fields only):
21     C Interpolation (linear in p^kappa) to standard pressure levels
22     C
23    
24     C !USES:
25     IMPLICIT NONE
26     #include "SIZE.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29     #include "GRID.h"
30     #include "DIAGNOSTICS_SIZE.h"
31     #include "DIAGNOSTICS.h"
32    
33     #ifdef ALLOW_FIZHI
34     #include "fizhi_SIZE.h"
35     #else
36     INTEGER Nrphys
37     PARAMETER (Nrphys=0)
38     #endif
39    
40    
41     C !INPUT PARAMETERS:
42     C listId :: Diagnostics list number being written
43     C md :: field number in the list "listId".
44     C ndId :: diagnostics Id number (in available diagnostics list)
45     C ip :: diagnostics pointer to storage array
46     C im :: counter-mate pointer to storage array
47     C nlevsout:: number of levels
48     C qtmp1 :: diagnostics field output array
49     C undef ::
50     C myTime :: current time of simulation (s)
51     C myIter :: current iteration number
52     C myThid :: my Thread Id number
53     INTEGER listId, md, ndId, ip, im
54     INTEGER nlevsout
55     _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)
56     _RL undef
57     _RL myTime
58     INTEGER myIter, myThid
59     CEOP
60    
61     C !LOCAL VARIABLES:
62     C i,j,k :: loop indices
63     INTEGER i, j, k
64     INTEGER bi, bj
65     _RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66     _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)
67     _RL getcon
68     EXTERNAL getcon
69 molod 1.2 integer nplevs1, nplevs2, nplevs3
70     parameter(nplevs1 = 16)
71     parameter(nplevs2 = 16)
72     parameter(nplevs3 = 17)
73     _RL plevs1(nplevs1)
74 jmc 1.1 data plevs1/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,
75     . 600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,
76     . 250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,
77     . 70.0 _d 2, 50.0 _d 2, 30.0 _d 2, 20.0 _d 2/
78 molod 1.2 _RL plevs2(nplevs2)
79 jmc 1.1 data plevs2/ 1000.0 _d 2, 950.0 _d 2, 900.0 _d 2, 850.0 _d 2,
80     . 800.0 _d 2, 750.0 _d 2, 700.0 _d 2, 600.0 _d 2,
81     . 500.0 _d 2, 400.0 _d 2, 300.0 _d 2, 250.0 _d 2,
82     . 200.0 _d 2, 150.0 _d 2, 100.0 _d 2, 50.0 _d 2/
83 molod 1.2 _RL plevs3(nplevs3)
84     data plevs3/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,
85     . 600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,
86     . 250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,
87     . 70.0 _d 2, 50.0 _d 2, 30.0 _d 2, 20.0 _d 2,
88     . 10.0 _d 2/
89     C Use the biggest of nplevs 1-3 (or any others) for the size of qprs
90     _RL qprs(sNx,sNy,nplevs3)
91 jmc 1.1 _RL qinp(sNx,sNy,Nr+Nrphys)
92     _RL pkz(sNx,sNy,Nr+Nrphys)
93     _RL pksrf(sNx,sNy)
94     _RL p
95     _RL kappa
96     #ifdef NONLIN_FRSURF
97     INTEGER jpoint1,ipoint1
98     INTEGER jpoint2,ipoint2
99     logical foundp
100     data foundp /.false./
101     #endif
102     CHARACTER*(MAX_LEN_MBUF) msgBuf
103    
104     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105    
106     IF(fflags(listId)(2:2).eq.'P') then
107     kappa = getcon('KAPPA')
108    
109     c If nonlinear free surf is active, need averaged pressures
110     #ifdef NONLIN_FRSURF
111     if(select_rStar.GT.0)then
112     call diagnostics_get_pointers('RSURF ',ipoint1,jpoint1,
113     . myThid)
114     call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,
115     . myThid)
116     C if fizhi is being used, may need to get physics grid pressures
117     #ifdef ALLOW_FIZHI
118     if(gdiag(ndId)(10:10) .EQ. 'L')then
119     call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,
120     . myThid)
121     endif
122     #endif
123 molod 1.2 if( jpoint1.ne.0 .and. jpoint2.ne.0) then
124     foundp = .true.
125     else
126     foundp = .false.
127     endif
128 jmc 1.1
129     if(.not. foundp) then
130     WRITE(msgBuf,'(3A)') 'DIAGNOSTICS_INTERP_VERT: ',
131     . ' Have asked for pressure interpolation but have not ',
132     . ' Activated surface and 3D pressure diagnostic, ',
133     . ' RSURF and PRESSURE'
134     CALL PRINT_ERROR( msgBuf , myThid )
135     STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'
136     endif
137    
138     DO bj = myByLo(myThid), myByHi(myThid)
139     DO bi = myBxLo(myThid), myBxHi(myThid)
140     call getdiag(1.,undef,qtmpsrf(1-OLx,1-OLy,bi,bj),
141     . jpoint1,0,ipoint1,0,bi,bj,myThid)
142     ENDDO
143     ENDDO
144     DO bj = myByLo(myThid), myByHi(myThid)
145     DO bi = myBxLo(myThid), myBxHi(myThid)
146     DO k = 1,nlevels(listId)
147     call getdiag(levs(k,listId),undef,
148     . qtmp2(1-OLx,1-OLy,k,bi,bj),jpoint2,0,ipoint2,0,
149     . bi,bj,myThid)
150     ENDDO
151     ENDDO
152     ENDDO
153     endif
154     #else
155     C If nonlinear free surf is off, get pressures from rC and rF arrays
156     DO bj = myByLo(myThid), myByHi(myThid)
157     DO bi = myBxLo(myThid), myBxHi(myThid)
158     DO j = 1-OLy,sNy+OLy
159     DO i = 1-OLx,sNx+OLx
160     qtmpsrf(i,j,bi,bj) = rF(1)
161     ENDDO
162     ENDDO
163     DO j = 1-OLy,sNy+OLy
164     DO i = 1-OLx,sNx+OLx
165     DO k = 1,nlevels(listId)
166     qtmp2(i,j,k,bi,bj) = rC(k)
167     ENDDO
168     ENDDO
169     ENDDO
170     ENDDO
171     ENDDO
172     #endif
173     C Load p to the kappa into a temporary array
174     DO bj = myByLo(myThid), myByHi(myThid)
175     DO bi = myBxLo(myThid), myBxHi(myThid)
176     DO j = 1,sNy
177     DO i = 1,sNx
178     pksrf(i,j) = qtmpsrf(i,j,bi,bj) ** kappa
179     DO k = 1,nlevels(listId)
180     if(gdiag(ndId)(10:10).eq.'R') then
181     if(hFacC(i,j,nlevels(listId)-k+1,bi,bj).ne.0.) then
182     qinp(i,j,k) = qtmp1(i,j,nlevels(listId)-k+1,bi,bj)
183     else
184     qinp(i,j,k) = undef
185     endif
186     pkz(i,j,k) = qtmp2(i,j,nlevels(listId)-k+1,bi,bj)**kappa
187     elseif(gdiag(ndId)(10:10).eq.'L') then
188     qinp(i,j,k) = qtmp1(i,j,k,bi,bj)
189     pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa
190     endif
191     ENDDO
192     ENDDO
193     ENDDO
194    
195 molod 1.2 if(fflags(listId)(3:3).eq.'1') then
196     nlevsout = nplevs1
197     DO k = 1,nplevs1
198 jmc 1.1 p = plevs1(k)
199 molod 1.2 call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,
200     . nlevels(listId) )
201     ENDDO
202     elseif(fflags(listId)(3:3).eq.'2')then
203     nlevsout = nplevs2
204     DO k = 1,nplevs2
205 jmc 1.1 p = plevs2(k)
206 molod 1.2 call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,
207 jmc 1.1 . nlevels(listId) )
208 molod 1.2 ENDDO
209     elseif(fflags(listId)(3:3).eq.'3')then
210     nlevsout = nplevs3
211     DO k = 1,nplevs3
212     p = plevs3(k)
213     call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,
214     . nlevels(listId) )
215     ENDDO
216     endif
217 jmc 1.1
218     DO j = 1,sNy
219     DO i = 1,sNx
220     DO k = 1,nlevsout
221     qtmp1(i,j,k,bi,bj) = qprs(i,j,k)
222     if(qtmp1(i,j,k,bi,bj).eq.undef) qtmp1(i,j,k,bi,bj) = 0.
223     ENDDO
224     ENDDO
225     ENDDO
226    
227     ENDDO
228     ENDDO
229    
230     ENDIF
231    
232     RETURN
233     END

  ViewVC Help
Powered by ViewVC 1.1.22