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