/[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.1 - (hide annotations) (download)
Thu Aug 25 21:54:54 2005 UTC (18 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint57s_post, checkpoint57v_post, checkpoint57r_post, checkpint57u_post
do vertical interpolation in dedicated S/R DIAGNOSTICS_INTERP_VERT

1 jmc 1.1 C $Header: $
2     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     integer nplevs
70     parameter(nplevs = 16)
71     _RL plevs1(nplevs)
72     data plevs1/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,
73     . 600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,
74     . 250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,
75     . 70.0 _d 2, 50.0 _d 2, 30.0 _d 2, 20.0 _d 2/
76     _RL plevs2(nplevs)
77     data plevs2/ 1000.0 _d 2, 950.0 _d 2, 900.0 _d 2, 850.0 _d 2,
78     . 800.0 _d 2, 750.0 _d 2, 700.0 _d 2, 600.0 _d 2,
79     . 500.0 _d 2, 400.0 _d 2, 300.0 _d 2, 250.0 _d 2,
80     . 200.0 _d 2, 150.0 _d 2, 100.0 _d 2, 50.0 _d 2/
81     _RL qprs(sNx,sNy,nplevs)
82     _RL qinp(sNx,sNy,Nr+Nrphys)
83     _RL pkz(sNx,sNy,Nr+Nrphys)
84     _RL pksrf(sNx,sNy)
85     _RL p
86     _RL kappa
87     #ifdef NONLIN_FRSURF
88     INTEGER jpoint1,ipoint1
89     INTEGER jpoint2,ipoint2
90     logical foundp
91     data foundp /.false./
92     #endif
93     CHARACTER*(MAX_LEN_MBUF) msgBuf
94    
95     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
96    
97     IF(fflags(listId)(2:2).eq.'P') then
98     kappa = getcon('KAPPA')
99    
100     c If nonlinear free surf is active, need averaged pressures
101     #ifdef NONLIN_FRSURF
102     if(select_rStar.GT.0)then
103     call diagnostics_get_pointers('RSURF ',ipoint1,jpoint1,
104     . myThid)
105     call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,
106     . myThid)
107     C if fizhi is being used, may need to get physics grid pressures
108     #ifdef ALLOW_FIZHI
109     if(gdiag(ndId)(10:10) .EQ. 'L')then
110     call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,
111     . myThid)
112     endif
113     #endif
114     if( jpoint1.ne.0 .and. jpoint2.ne.0) foundp = .true.
115    
116     if(.not. foundp) then
117     WRITE(msgBuf,'(3A)') 'DIAGNOSTICS_INTERP_VERT: ',
118     . ' Have asked for pressure interpolation but have not ',
119     . ' Activated surface and 3D pressure diagnostic, ',
120     . ' RSURF and PRESSURE'
121     CALL PRINT_ERROR( msgBuf , myThid )
122     STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'
123     endif
124    
125     DO bj = myByLo(myThid), myByHi(myThid)
126     DO bi = myBxLo(myThid), myBxHi(myThid)
127     call getdiag(1.,undef,qtmpsrf(1-OLx,1-OLy,bi,bj),
128     . jpoint1,0,ipoint1,0,bi,bj,myThid)
129     ENDDO
130     ENDDO
131     DO bj = myByLo(myThid), myByHi(myThid)
132     DO bi = myBxLo(myThid), myBxHi(myThid)
133     DO k = 1,nlevels(listId)
134     call getdiag(levs(k,listId),undef,
135     . qtmp2(1-OLx,1-OLy,k,bi,bj),jpoint2,0,ipoint2,0,
136     . bi,bj,myThid)
137     ENDDO
138     ENDDO
139     ENDDO
140     endif
141     #else
142     C If nonlinear free surf is off, get pressures from rC and rF arrays
143     DO bj = myByLo(myThid), myByHi(myThid)
144     DO bi = myBxLo(myThid), myBxHi(myThid)
145     DO j = 1-OLy,sNy+OLy
146     DO i = 1-OLx,sNx+OLx
147     qtmpsrf(i,j,bi,bj) = rF(1)
148     ENDDO
149     ENDDO
150     DO j = 1-OLy,sNy+OLy
151     DO i = 1-OLx,sNx+OLx
152     DO k = 1,nlevels(listId)
153     qtmp2(i,j,k,bi,bj) = rC(k)
154     ENDDO
155     ENDDO
156     ENDDO
157     ENDDO
158     ENDDO
159     #endif
160     C Load p to the kappa into a temporary array
161     nlevsout = nplevs
162     DO bj = myByLo(myThid), myByHi(myThid)
163     DO bi = myBxLo(myThid), myBxHi(myThid)
164     DO j = 1,sNy
165     DO i = 1,sNx
166     pksrf(i,j) = qtmpsrf(i,j,bi,bj) ** kappa
167     DO k = 1,nlevels(listId)
168     if(gdiag(ndId)(10:10).eq.'R') then
169     if(hFacC(i,j,nlevels(listId)-k+1,bi,bj).ne.0.) then
170     qinp(i,j,k) = qtmp1(i,j,nlevels(listId)-k+1,bi,bj)
171     else
172     qinp(i,j,k) = undef
173     endif
174     pkz(i,j,k) = qtmp2(i,j,nlevels(listId)-k+1,bi,bj)**kappa
175     elseif(gdiag(ndId)(10:10).eq.'L') then
176     qinp(i,j,k) = qtmp1(i,j,k,bi,bj)
177     pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa
178     endif
179     ENDDO
180     ENDDO
181     ENDDO
182    
183     DO k = 1,nplevs
184     if(fflags(listId)(3:3).eq.'1') then
185     p = plevs1(k)
186     elseif(fflags(listId)(3:3).eq.'2')then
187     p = plevs2(k)
188     endif
189     call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,
190     . nlevels(listId) )
191     ENDDO
192    
193     DO j = 1,sNy
194     DO i = 1,sNx
195     DO k = 1,nlevsout
196     qtmp1(i,j,k,bi,bj) = qprs(i,j,k)
197     if(qtmp1(i,j,k,bi,bj).eq.undef) qtmp1(i,j,k,bi,bj) = 0.
198     ENDDO
199     ENDDO
200     ENDDO
201    
202     ENDDO
203     ENDDO
204    
205     ENDIF
206    
207     RETURN
208     END

  ViewVC Help
Powered by ViewVC 1.1.22