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

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

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

revision 1.3 by molod, Tue Oct 25 16:27:29 2005 UTC revision 1.8 by jmc, Wed Jan 31 21:47:55 2007 UTC
# Line 9  C     !ROUTINE: DIAGNOSTICS_INTERP_VERT Line 9  C     !ROUTINE: DIAGNOSTICS_INTERP_VERT
9    
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE  DIAGNOSTICS_INTERP_VERT(        SUBROUTINE  DIAGNOSTICS_INTERP_VERT(
12       I     listId, md, ndId, ip, im,       I                        listId, md, ndId, ip, im, lm,
13       U     nlevsout,       U                        qtmp1,
14       U     qtmp1,       I                        undef,
15       I     undef,       I                        myTime, myIter, myThid )
      I     myTime, myIter, myThid )  
16    
17  C     !DESCRIPTION:  C     !DESCRIPTION:
18  C     Interpolate vertically a diagnostics field before writing to file.  C     Interpolate vertically a diagnostics field before writing to file.
19  C       presently implemented (for Atmospheric fields only):  C       presently implemented (for Atmospheric fields only):
20  C         Interpolation (linear in p^kappa) to standard pressure levels  C         Interpolation (linear in p^kappa) to standard pressure levels
21  C                                C
22    
23  C     !USES:  C     !USES:
24        IMPLICIT NONE        IMPLICIT NONE
# Line 30  C     !USES: Line 29  C     !USES:
29  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
30  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
31    
32  #ifdef ALLOW_FIZHI        INTEGER NrMax
33  #include "fizhi_SIZE.h"        PARAMETER( NrMax = numLevels )
 #else  
       INTEGER Nrphys  
       PARAMETER (Nrphys=0)  
 #endif  
34    
35    
36  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
# Line 44  C     md      :: field number in the lis Line 39  C     md      :: field number in the lis
39  C     ndId    :: diagnostics  Id number (in available diagnostics list)  C     ndId    :: diagnostics  Id number (in available diagnostics list)
40  C     ip      :: diagnostics  pointer to storage array  C     ip      :: diagnostics  pointer to storage array
41  C     im      :: counter-mate pointer to storage array  C     im      :: counter-mate pointer to storage array
42  C     nlevsout:: number of levels  C     lm      :: index in the averageCycle
43  C     qtmp1   :: diagnostics field output array  C     qtmp1   :: diagnostics field output array
44  C     undef   ::  C     undef   ::
45  C     myTime  :: current time of simulation (s)  C     myTime  :: current time of simulation (s)
46  C     myIter  :: current iteration number  C     myIter  :: current iteration number
47  C     myThid  :: my Thread Id number  C     myThid  :: my Thread Id number
48        INTEGER listId, md, ndId, ip, im        INTEGER listId, md, ndId, ip, im, lm
49        INTEGER nlevsout        _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
       _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)  
50        _RL     undef        _RL     undef
51        _RL     myTime        _RL     myTime
52        INTEGER myIter, myThid        INTEGER myIter, myThid
# Line 62  C     !LOCAL VARIABLES: Line 56  C     !LOCAL VARIABLES:
56  C     i,j,k :: loop indices  C     i,j,k :: loop indices
57        INTEGER i, j, k        INTEGER i, j, k
58        INTEGER bi, bj        INTEGER bi, bj
59        _RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL   qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
60        _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _RL   qtmp2  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
61        _RL getcon        _RL   getcon
62        EXTERNAL getcon        EXTERNAL getcon
63        integer nplevs1, nplevs2, nplevs3        INTEGER kLev
64        parameter(nplevs1 = 16)        _RL   qprs (sNx,sNy)
65        parameter(nplevs2 = 16)        _RL   qinp (sNx,sNy,NrMax)
66        parameter(nplevs3 = 17)        _RL   pkz  (sNx,sNy,NrMax)
67        _RL plevs1(nplevs1)        _RL   pksrf(sNx,sNy)
68        data plevs1/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,        _RL   pk, pkTop, tmpLev
69       .              600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,        _RL   kappa
      .              250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  
      .               70.0 _d 2,  50.0 _d 2,  30.0 _d 2,  20.0 _d 2/  
       _RL plevs2(nplevs2)  
       data plevs2/ 1000.0 _d 2, 950.0 _d 2, 900.0 _d 2, 850.0 _d 2,  
      .              800.0 _d 2, 750.0 _d 2, 700.0 _d 2, 600.0 _d 2,  
      .              500.0 _d 2, 400.0 _d 2, 300.0 _d 2, 250.0 _d 2,  
      .              200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  50.0 _d 2/  
       _RL plevs3(nplevs3)  
       data plevs3/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,  
      .              600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,  
      .              250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  
      .               70.0 _d 2,  50.0 _d 2,  30.0 _d 2,  20.0 _d 2,  
      .               10.0 _d 2/  
 C Use the biggest of nplevs 1-3 (or any others) for the size of qprs  
       _RL qprs(sNx,sNy,nplevs3)  
       _RL qinp(sNx,sNy,Nr+Nrphys)  
       _RL pkz(sNx,sNy,Nr+Nrphys)  
       _RL pksrf(sNx,sNy)  
       _RL p  
       _RL kappa  
 #ifdef NONLIN_FRSURF  
70        INTEGER jpoint1,ipoint1        INTEGER jpoint1,ipoint1
71        INTEGER jpoint2,ipoint2        INTEGER jpoint2,ipoint2
72        logical foundp        LOGICAL pInc
       data foundp /.false./  
 #endif  
73        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
74    
75  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
76    
77          IF(fflags(listId)(2:2).eq.'P') then        IF (fflags(listId)(2:2).EQ.'P') THEN
78            kappa = getcon('KAPPA')          kappa = getcon('KAPPA')
79            pkTop = 0. _d 0
80  c If nonlinear free surf is active, need averaged pressures  
81  #ifdef NONLIN_FRSURF  C--   If nonlinear free surf is active, need averaged pressures
82            if(select_rStar.GT.0)then          IF (select_rStar.GT.0) THEN
83             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,            CALL DIAGNOSTICS_GET_POINTERS( 'RSURF   ', listId,
84       .                                                           myThid)       &                                   jpoint1, ipoint1, myThid )
85             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,  C-    IF fizhi is being  used, may need to get physics grid pressures
86       .                                                           myThid)            IF ( useFIZHI .AND.
87  C if fizhi is being  used, may need to get physics grid pressures       &          gdiag(ndId)(10:10) .EQ. 'L') THEN
88  #ifdef ALLOW_FIZHI             CALL DIAGNOSTICS_GET_POINTERS('FIZPRES ', listId,
89             if(gdiag(ndId)(10:10) .EQ. 'L')then       &                                   jpoint2, ipoint2, myThid )
90             call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,            ELSE
91       .                                                           myThid)             CALL DIAGNOSTICS_GET_POINTERS('RCENTER ', listId,
92             endif       &                                   jpoint2, ipoint2, myThid )
93  #endif            ENDIF
94             if( jpoint1.ne.0 .and. jpoint2.ne.0) then            IF ( ipoint1.EQ.0 .OR. ipoint2.EQ.0 ) THEN
95              foundp = .true.              WRITE(msgBuf,'(2A,I4,2A)') 'DIAGNOSTICS_INTERP_VERT: ',
96             else       &      'fails to interpolate diag.(#', ndId,'): ',flds(md,listId)
             foundp = .false.  
            endif  
   
            if(.not. foundp) then  
             WRITE(msgBuf,'(3A)') 'DIAGNOSTICS_INTERP_VERT: ',  
      .    ' Have asked for pressure interpolation but have not ',  
      .    ' Activated surface and 3D pressure diagnostic, ',  
      .    ' RSURF and PRESSURE'  
97              CALL PRINT_ERROR( msgBuf , myThid )              CALL PRINT_ERROR( msgBuf , myThid )
98              STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'              STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'
99             endif            ENDIF
100    C-    averageCycle: move pointer
101              ipoint1 = ipoint1 + kdiag(jpoint1)*(lm-1)
102              ipoint2 = ipoint2 + kdiag(jpoint2)*(lm-1)
103    
104             DO bj = myByLo(myThid), myByHi(myThid)            DO bj = myByLo(myThid), myByHi(myThid)
105              DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
106               call getdiag(1.,undef,qtmpsrf(1-OLx,1-OLy,bi,bj),               tmpLev = 1. _d 0
107       .                       jpoint1,0,ipoint1,0,bi,bj,myThid)               CALL GETDIAG( tmpLev,undef,
108              ENDDO       O                     qtmpsrf(1-OLx,1-OLy,bi,bj),
109             ENDDO       I                     jpoint1,0,ipoint1,0, bi,bj,myThid )
110             DO bj = myByLo(myThid), myByHi(myThid)  c            WRITE(0,*) 'rSurf:',bi,bj,qtmpsrf(15,15,bi,bj)
111              DO bi = myBxLo(myThid), myBxHi(myThid)               DO k = 1,kdiag(jpoint2)
112               DO k = 1,nlevels(listId)                tmpLev = k
113                call getdiag(levs(k,listId),undef,                CALL GETDIAG(tmpLev,undef,
114       .          qtmp2(1-OLx,1-OLy,k,bi,bj),jpoint2,0,ipoint2,0,       O                     qtmp2(1-OLx,1-OLy,k,bi,bj),
115       .          bi,bj,myThid)       I                     jpoint2,0,ipoint2,0, bi,bj,myThid )
116               ENDDO               ENDDO
             ENDDO  
117             ENDDO             ENDDO
118            endif            ENDDO
119  #else  
120  C If nonlinear free surf is off, get pressures from rC and rF arrays          ELSE
121    C-    If nonlinear free surf is off, get pressures from rC and rF arrays
122    
123            DO bj = myByLo(myThid), myByHi(myThid)            DO bj = myByLo(myThid), myByHi(myThid)
124             DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
125              DO j = 1-OLy,sNy+OLy              DO j = 1-OLy,sNy+OLy
126               DO i = 1-OLx,sNx+OLx               DO i = 1-OLx,sNx+OLx
127                qtmpsrf(i,j,bi,bj) = rF(1)                 qtmpsrf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
128               ENDDO               ENDDO
129              ENDDO              ENDDO
130              DO j = 1-OLy,sNy+OLy              DO k = 1,kdiag(ndId)
131               DO i = 1-OLx,sNx+OLx               DO j = 1-OLy,sNy+OLy
132                DO k = 1,nlevels(listId)                DO i = 1-OLx,sNx+OLx
133                 qtmp2(i,j,k,bi,bj) = rC(k)                 qtmp2(i,j,k,bi,bj) = rC(k)
134                ENDDO                ENDDO
135               ENDDO               ENDDO
136              ENDDO              ENDDO
137             ENDDO             ENDDO
138            ENDDO            ENDDO
139  #endif  
140  C Load p to the kappa into a temporary array  C-    end if nonlinear/linear free-surf
141            DO bj = myByLo(myThid), myByHi(myThid)          ENDIF
142             DO bi = myBxLo(myThid), myBxHi(myThid)  
143              DO j = 1,sNy  C--   start loops on tile indices bi,bj:
144               DO i = 1,sNx          DO bj = myByLo(myThid), myByHi(myThid)
145                pksrf(i,j) = qtmpsrf(i,j,bi,bj) ** kappa           DO bi = myBxLo(myThid), myBxHi(myThid)
146                DO k = 1,nlevels(listId)  C-    Load p to the kappa into a temporary array
147                 if(gdiag(ndId)(10:10).eq.'R') then             DO j = 1,sNy
148                  if(hFacC(i,j,nlevels(listId)-k+1,bi,bj).ne.0.) then              DO i = 1,sNx
149                   qinp(i,j,k) =  qtmp1(i,j,nlevels(listId)-k+1,bi,bj)                pksrf(i,j)    = qtmpsrf(i,j,bi,bj)**kappa
150                  else              ENDDO
151                   qinp(i,j,k) =  undef             ENDDO
152                  endif             IF ( useFIZHI.AND.gdiag(ndId)(10:10).EQ.'L') THEN
153                  pkz(i,j,k) = qtmp2(i,j,nlevels(listId)-k+1,bi,bj)**kappa              pInc = .TRUE.
154                 elseif(gdiag(ndId)(10:10).eq.'L') then              DO k = 1,kdiag(ndId)
155                  qinp(i,j,k) =  qtmp1(i,j,k,bi,bj)               DO j = 1,sNy
156                  pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa                DO i = 1,sNx
157                 endif                  qinp(i,j,k) = qtmp1(i,j,k,bi,bj)
158                    pkz(i,j,k)  = qtmp2(i,j,k,bi,bj)**kappa
159                ENDDO                ENDDO
160               ENDDO               ENDDO
161              ENDDO              ENDDO
162               ELSE
163              if(fflags(listId)(3:3).eq.'1') then              DO k = 1,kdiag(ndId)
164               nlevsout = nplevs1               pInc = .TRUE.
165               DO k = 1,nplevs1               kLev = kdiag(ndId)-k+1
166                p = plevs1(k)  c            pInc = .FALSE.
167                call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,  c            kLev = k
168       .                                     nlevels(listId),myThid )               DO j = 1,sNy
169               ENDDO                DO i = 1,sNx
170              elseif(fflags(listId)(3:3).eq.'2')then                  IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
171               nlevsout = nplevs2                   qinp(i,j,k)= qtmp1(i,j,kLev,bi,bj)
172               DO k = 1,nplevs2                  ELSE
173                p = plevs2(k)                   qinp(i,j,k)= undef
174                call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,                  ENDIF
175       .                                     nlevels(listId),myThid )                  pkz(i,j,k)  = qtmp2(i,j,kLev,bi,bj)**kappa
              ENDDO  
             elseif(fflags(listId)(3:3).eq.'3')then  
              nlevsout = nplevs3  
              DO k = 1,nplevs3  
               p = plevs3(k)  
               call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,  
      .                                     nlevels(listId),myThid )  
              ENDDO  
             endif  
   
             DO j = 1,sNy  
              DO i = 1,sNx  
               DO k = 1,nlevsout  
                qtmp1(i,j,k,bi,bj) =  qprs(i,j,k)  
                if(qtmp1(i,j,k,bi,bj).eq.undef) qtmp1(i,j,k,bi,bj) = 0.  
176                ENDDO                ENDDO
177               ENDDO               ENDDO
178              ENDDO              ENDDO
179               ENDIF
180    
181    C-    Interpolate, level per level, and put interpolated field in qprs:
182               DO k = 1,nlevels(listId)
183                 pk = levs(k,listId)**kappa
184                 CALL DIAGNOSTICS_INTERP_P2P(
185         O                        qprs,
186         I                        qinp,pkz,pksrf,pkTop,pk,
187         I                        undef,pInc,sNx*sNy,kdiag(ndId),myThid )
188    C-    Transfert qprs to qtmp1:
189                 DO j = 1,sNy
190                  DO i = 1,sNx
191                   IF (qprs(i,j).EQ.undef) THEN
192                     qtmp1(i,j,k,bi,bj) = 0.
193                   ELSE
194                     qtmp1(i,j,k,bi,bj) =  qprs(i,j)
195                   ENDIF
196                  ENDDO
197                 ENDDO
198             ENDDO             ENDDO
           ENDDO  
199    
200          ENDIF  C-   end bi,bj loops
201             ENDDO
202            ENDDO
203    
204          ENDIF
205    
206        RETURN        RETURN
207        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22