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