pkg/diagnostics/diagnostics_interp_vert.F

C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_interp_vert.F,v 1.6 2006/11/19 22:03:13 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C     !ROUTINE: DIAGNOSTICS_INTERP_VERT

C     !INTERFACE:
      SUBROUTINE  DIAGNOSTICS_INTERP_VERT(
     I                        listId, md, ndId, ip, im, lm,
     U                        qtmp1,
     I                        undef,
     I                        myTime, myIter, myThid )

C     !DESCRIPTION:
C     Interpolate vertically a diagnostics field before writing to file.
C       presently implemented (for Atmospheric fields only):
C         Interpolation (linear in p^kappa) to standard pressure levels
C

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"

      INTEGER NrMax
#ifdef ALLOW_FIZHI
#include "fizhi_SIZE.h"
      PARAMETER( NrMax = Nr+Nrphys )
#else
      PARAMETER( NrMax = Nr )
#endif


C     !INPUT PARAMETERS:
C     listId  :: Diagnostics list number being written
C     md      :: field number in the list "listId".
C     ndId    :: diagnostics  Id number (in available diagnostics list)
C     ip      :: diagnostics  pointer to storage array
C     im      :: counter-mate pointer to storage array
C     lm      :: index in the averageCycle
C     qtmp1   :: diagnostics field output array
C     undef   ::
C     myTime  :: current time of simulation (s)
C     myIter  :: current iteration number
C     myThid  :: my Thread Id number
      INTEGER listId, md, ndId, ip, im, lm
      _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
      _RL     undef
      _RL     myTime
      INTEGER myIter, myThid
CEOP

C     !LOCAL VARIABLES:
C     i,j,k :: loop indices
      INTEGER i, j, k
      INTEGER bi, bj
      _RL   qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL   qtmp2  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
      _RL   getcon
      EXTERNAL getcon
      INTEGER kLev
      _RL   qprs (sNx,sNy)
      _RL   qinp (sNx,sNy,NrMax)
      _RL   pkz  (sNx,sNy,NrMax)
      _RL   pksrf(sNx,sNy)
      _RL   pk, pkTop, tmpLev
      _RL   kappa
      INTEGER jpoint1,ipoint1
      INTEGER jpoint2,ipoint2
      LOGICAL pInc
      CHARACTER*(MAX_LEN_MBUF) msgBuf

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF (fflags(listId)(2:2).EQ.'P') THEN
        kappa = getcon('KAPPA')
        pkTop = 0. _d 0

C--   If nonlinear free surf is active, need averaged pressures
        IF (select_rStar.GT.0) THEN
          CALL DIAGNOSTICS_GET_POINTERS( 'RSURF   ', listId,
     &                                   jpoint1, ipoint1, myThid )
C-    IF fizhi is being  used, may need to get physics grid pressures
          IF ( useFIZHI .AND.
     &          gdiag(ndId)(10:10) .EQ. 'L') THEN
           CALL DIAGNOSTICS_GET_POINTERS('FIZPRES ', listId,
     &                                   jpoint2, ipoint2, myThid )
          ELSE
           CALL DIAGNOSTICS_GET_POINTERS('RCENTER ', listId,
     &                                   jpoint2, ipoint2, myThid )
          ENDIF
          IF ( ipoint1.EQ.0 .OR. ipoint2.EQ.0 ) THEN
            WRITE(msgBuf,'(2A,I4,2A)') 'DIAGNOSTICS_INTERP_VERT: ',
     &      'fails to interpolate diag.(#', ndId,'): ',flds(md,listId)
            CALL PRINT_ERROR( msgBuf , myThid )
            STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'
          ENDIF
C-    averageCycle: move pointer
          ipoint1 = ipoint1 + kdiag(jpoint1)*(lm-1)
          ipoint2 = ipoint2 + kdiag(jpoint2)*(lm-1)

          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
             tmpLev = 1. _d 0
             CALL GETDIAG( tmpLev,undef,
     O                     qtmpsrf(1-OLx,1-OLy,bi,bj),
     I                     jpoint1,0,ipoint1,0, bi,bj,myThid )
c            WRITE(0,*) 'rSurf:',bi,bj,qtmpsrf(15,15,bi,bj)
             DO k = 1,kdiag(jpoint2)
              tmpLev = k
              CALL GETDIAG(tmpLev,undef,
     O                     qtmp2(1-OLx,1-OLy,k,bi,bj),
     I                     jpoint2,0,ipoint2,0, bi,bj,myThid )
             ENDDO
           ENDDO
          ENDDO

        ELSE
C-    If nonlinear free surf is off, get pressures from rC and rF arrays

          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j = 1-OLy,sNy+OLy
             DO i = 1-OLx,sNx+OLx
               qtmpsrf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
             ENDDO
            ENDDO
            DO k = 1,kdiag(ndId)
             DO j = 1-OLy,sNy+OLy
              DO i = 1-OLx,sNx+OLx
               qtmp2(i,j,k,bi,bj) = rC(k)
              ENDDO
             ENDDO
            ENDDO
           ENDDO
          ENDDO

C-    end if nonlinear/linear free-surf
        ENDIF

C--   start loops on tile indices bi,bj:
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
C-    Load p to the kappa into a temporary array
           DO j = 1,sNy
            DO i = 1,sNx
              pksrf(i,j)    = qtmpsrf(i,j,bi,bj)**kappa
            ENDDO
           ENDDO
           IF ( useFIZHI.AND.gdiag(ndId)(10:10).EQ.'L') THEN
            pInc = .TRUE.
            DO k = 1,kdiag(ndId)
             DO j = 1,sNy
              DO i = 1,sNx
                qinp(i,j,k) = qtmp1(i,j,k,bi,bj)
                pkz(i,j,k)  = qtmp2(i,j,k,bi,bj)**kappa
              ENDDO
             ENDDO
            ENDDO
           ELSE
            DO k = 1,kdiag(ndId)
             pInc = .TRUE.
             kLev = kdiag(ndId)-k+1
c            pInc = .FALSE.
c            kLev = k
             DO j = 1,sNy
              DO i = 1,sNx
                IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
                 qinp(i,j,k)= qtmp1(i,j,kLev,bi,bj)
                ELSE
                 qinp(i,j,k)= undef
                ENDIF
                pkz(i,j,k)  = qtmp2(i,j,kLev,bi,bj)**kappa
              ENDDO
             ENDDO
            ENDDO
           ENDIF

C-    Interpolate, level per level, and put interpolated field in qprs:
           DO k = 1,nlevels(listId)
             pk = levs(k,listId)**kappa
             CALL DIAGNOSTICS_INTERP_P2P(
     O                        qprs,
     I                        qinp,pkz,pksrf,pkTop,pk,
     I                        undef,pInc,sNx*sNy,kdiag(ndId),myThid )
C-    Transfert qprs to qtmp1:
             DO j = 1,sNy
              DO i = 1,sNx
               IF (qprs(i,j).EQ.undef) THEN
                 qtmp1(i,j,k,bi,bj) = 0.
               ELSE
                 qtmp1(i,j,k,bi,bj) =  qprs(i,j)
               ENDIF
              ENDDO
             ENDDO
           ENDDO

C-   end bi,bj loops
         ENDDO
        ENDDO

      ENDIF

      RETURN
      END
1	jmc	1.7	C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_interp_vert.F,v 1.6 2006/11/19 22:03:13 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	jmc	1.7	I listId, md, ndId, ip, im, lm,
13			U qtmp1,
14			I undef,
15			I myTime, myIter, myThid )
16	jmc	1.1
17			C !DESCRIPTION:
18			C Interpolate vertically a diagnostics field before writing to file.
19	jmc	1.7	C presently implemented (for Atmospheric fields only):
20	jmc	1.1	C Interpolation (linear in p^kappa) to standard pressure levels
21	jmc	1.7	C
22	jmc	1.1
23			C !USES:
24			IMPLICIT NONE
25			#include "SIZE.h"
26			#include "EEPARAMS.h"
27			#include "PARAMS.h"
28			#include "GRID.h"
29			#include "DIAGNOSTICS_SIZE.h"
30			#include "DIAGNOSTICS.h"
31
32	jmc	1.7	INTEGER NrMax
33	jmc	1.1	#ifdef ALLOW_FIZHI
34			#include "fizhi_SIZE.h"
35	jmc	1.7	PARAMETER( NrMax = Nr+Nrphys )
36	jmc	1.1	#else
37	jmc	1.7	PARAMETER( NrMax = Nr )
38	jmc	1.1	#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	jmc	1.7	C lm :: index in the averageCycle
48	jmc	1.1	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	jmc	1.7	INTEGER listId, md, ndId, ip, im, lm
54			_RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
55	jmc	1.1	_RL undef
56			_RL myTime
57			INTEGER myIter, myThid
58			CEOP
59
60			C !LOCAL VARIABLES:
61			C i,j,k :: loop indices
62			INTEGER i, j, k
63			INTEGER bi, bj
64	jmc	1.7	_RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65			_RL qtmp2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66			_RL getcon
67	jmc	1.1	EXTERNAL getcon
68	jmc	1.7	INTEGER kLev
69			_RL qprs (sNx,sNy)
70			_RL qinp (sNx,sNy,NrMax)
71			_RL pkz (sNx,sNy,NrMax)
72			_RL pksrf(sNx,sNy)
73			_RL pk, pkTop, tmpLev
74			_RL kappa
75	jmc	1.1	INTEGER jpoint1,ipoint1
76			INTEGER jpoint2,ipoint2
77	jmc	1.7	LOGICAL pInc
78	jmc	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
79
80			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
81
82	jmc	1.7	IF (fflags(listId)(2:2).EQ.'P') THEN
83			kappa = getcon('KAPPA')
84			pkTop = 0. _d 0
85
86			C-- If nonlinear free surf is active, need averaged pressures
87			IF (select_rStar.GT.0) THEN
88			CALL DIAGNOSTICS_GET_POINTERS( 'RSURF ', listId,
89			& jpoint1, ipoint1, myThid )
90			C- IF fizhi is being used, may need to get physics grid pressures
91			IF ( useFIZHI .AND.
92			& gdiag(ndId)(10:10) .EQ. 'L') THEN
93			CALL DIAGNOSTICS_GET_POINTERS('FIZPRES ', listId,
94			& jpoint2, ipoint2, myThid )
95			ELSE
96			CALL DIAGNOSTICS_GET_POINTERS('RCENTER ', listId,
97			& jpoint2, ipoint2, myThid )
98			ENDIF
99			IF ( ipoint1.EQ.0 .OR. ipoint2.EQ.0 ) THEN
100			WRITE(msgBuf,'(2A,I4,2A)') 'DIAGNOSTICS_INTERP_VERT: ',
101			& 'fails to interpolate diag.(#', ndId,'): ',flds(md,listId)
102	jmc	1.1	CALL PRINT_ERROR( msgBuf , myThid )
103			STOP 'ABNORMAL END: S/R DIAGNOSTICS_INTERP_VERT'
104	jmc	1.7	ENDIF
105			C- averageCycle: move pointer
106			ipoint1 = ipoint1 + kdiag(jpoint1)*(lm-1)
107			ipoint2 = ipoint2 + kdiag(jpoint2)*(lm-1)
108	jmc	1.1
109	jmc	1.7	DO bj = myByLo(myThid), myByHi(myThid)
110			DO bi = myBxLo(myThid), myBxHi(myThid)
111			tmpLev = 1. _d 0
112			CALL GETDIAG( tmpLev,undef,
113			O qtmpsrf(1-OLx,1-OLy,bi,bj),
114			I jpoint1,0,ipoint1,0, bi,bj,myThid )
115			c WRITE(0,*) 'rSurf:',bi,bj,qtmpsrf(15,15,bi,bj)
116			DO k = 1,kdiag(jpoint2)
117			tmpLev = k
118			CALL GETDIAG(tmpLev,undef,
119			O qtmp2(1-OLx,1-OLy,k,bi,bj),
120			I jpoint2,0,ipoint2,0, bi,bj,myThid )
121	jmc	1.1	ENDDO
122			ENDDO
123	jmc	1.7	ENDDO
124
125			ELSE
126			C- If nonlinear free surf is off, get pressures from rC and rF arrays
127
128	jmc	1.1	DO bj = myByLo(myThid), myByHi(myThid)
129			DO bi = myBxLo(myThid), myBxHi(myThid)
130			DO j = 1-OLy,sNy+OLy
131			DO i = 1-OLx,sNx+OLx
132	jmc	1.7	qtmpsrf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
133	jmc	1.1	ENDDO
134			ENDDO
135	jmc	1.7	DO k = 1,kdiag(ndId)
136			DO j = 1-OLy,sNy+OLy
137			DO i = 1-OLx,sNx+OLx
138	jmc	1.1	qtmp2(i,j,k,bi,bj) = rC(k)
139			ENDDO
140			ENDDO
141			ENDDO
142			ENDDO
143			ENDDO
144	jmc	1.7
145			C- end if nonlinear/linear free-surf
146			ENDIF
147
148			C-- start loops on tile indices bi,bj:
149			DO bj = myByLo(myThid), myByHi(myThid)
150			DO bi = myBxLo(myThid), myBxHi(myThid)
151			C- Load p to the kappa into a temporary array
152			DO j = 1,sNy
153			DO i = 1,sNx
154			pksrf(i,j) = qtmpsrf(i,j,bi,bj)**kappa
155			ENDDO
156			ENDDO
157			IF ( useFIZHI.AND.gdiag(ndId)(10:10).EQ.'L') THEN
158			pInc = .TRUE.
159			DO k = 1,kdiag(ndId)
160			DO j = 1,sNy
161			DO i = 1,sNx
162			qinp(i,j,k) = qtmp1(i,j,k,bi,bj)
163			pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa
164	jmc	1.1	ENDDO
165			ENDDO
166			ENDDO
167	jmc	1.7	ELSE
168			DO k = 1,kdiag(ndId)
169			pInc = .TRUE.
170			kLev = kdiag(ndId)-k+1
171			c pInc = .FALSE.
172			c kLev = k
173			DO j = 1,sNy
174			DO i = 1,sNx
175			IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
176			qinp(i,j,k)= qtmp1(i,j,kLev,bi,bj)
177			ELSE
178			qinp(i,j,k)= undef
179			ENDIF
180			pkz(i,j,k) = qtmp2(i,j,kLev,bi,bj)**kappa
181	jmc	1.1	ENDDO
182			ENDDO
183			ENDDO
184	jmc	1.7	ENDIF
185	jmc	1.1
186	jmc	1.7	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	jmc	1.1	ENDDO
204
205	jmc	1.7	C- end bi,bj loops
206			ENDDO
207			ENDDO
208
209			ENDIF
210	jmc	1.1
211			RETURN
212			END