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

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

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

revision 1.9 by jmc, Tue Nov 18 21:41:07 2008 UTC revision 1.11 by jmc, Thu Sep 3 20:46:12 2009 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: DIAG_VEGTILE_FILL  C     !ROUTINE: DIAG_VEGTILE_FILL
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE DIAG_VEGTILE_FILL (field,indx,chfr,ib,numpts,npeice,        SUBROUTINE DIAG_VEGTILE_FILL(
10       . check, chardiag, kLev, nLevs, bi, bj, myThid)       &                field,indx,chfr,ib,numpts,npeice,
11         &                check, chardiag, kLev, nLevs, bi, bj, myThid )
12  C     !DESCRIPTION:  C     !DESCRIPTION:
13  C***********************************************************************  C***********************************************************************
14  C Increment the diagnostics array with a vegetation tile space field  C Increment the diagnostics array with a vegetation tile space field
# Line 23  C     == Global variables === Line 24  C     == Global variables ===
24    
25  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
26  C***********************************************************************  C***********************************************************************
27  C   field  - array to be mapped to grid space [ib,levs] and added to qdiag  C   field    :: array to be mapped to grid space [ib,levs] and added to qdiag
28  C   indx   - array of horizontal indeces of grid points to convert to  C   indx     :: array of horizontal indices of grid points to convert to
29  C            tile space[numpts]  C               tile space[numpts]
30  C   chfr   - fractional area covered by the tile [ib]  C   chfr     :: fractional area covered by the tile [ib]
31  C   ib     - inner dimension of source array AND number of points in  C   ib       :: inner dimension of source array and number of points in
32  C            array a that need to be pasted  C               array a that need to be pasted
33  C   numpts - total number of points which were stripped  C   numpts   :: total number of points which were stripped
34  C   npeice - the current strip number to be filled  C   npeice   :: the current strip number to be filled
35  C   check  - logical to check for undefined values  C   check    :: logical to check for undefined values
36  C   chardiag ... Character expression for diag to fill  C   chardiag :: Character expression for diag to fill
37  C   kLev   ..... Integer flag for vertical levels:  C   kLev     :: Integer flag for vertical levels:
38  C                > 0 (any integer): WHICH single level to increment  C                > 0 (any integer): which single level to increment
39  C                0,-1 to increment "nLevs" levels in qdiag:  C                0,-1 to increment "nLevs" levels in qdiag:
40  C                 0 : fill-in in the same order as the input array  C                  0 : fill-in in the same order as the input array
41  C                -1 : fill-in in reverse order.  C                 -1 : fill-in in reverse order.
42  C   nLevs ...... indicates Number of levels of the input field array  C   nLevs    :: indicates Number of levels of the input field array
43  C   bi    ...... X-direction tile number  C   bi       :: X-direction tile number
44  C   bj    ...... Y-direction tile number  C   bj       :: Y-direction tile number
45  C   myThid     ::  my thread Id number  C   myThid   :: my thread Id number
 C  
 c IMPORTANT NOTE:  
 c  
 c This routine will result in roundoff differences if called from  
 c within a parallel region.  
46  C***********************************************************************  C***********************************************************************
47        CHARACTER*8 chardiag        CHARACTER*8 chardiag
48        INTEGER kLev, nLevs, bi, bj        INTEGER kLev, nLevs, bi, bj
49        INTEGER myThid        INTEGER myThid
50        integer ib,numpts,npeice        INTEGER ib,numpts,npeice
51        integer indx(numpts)        INTEGER indx(numpts)
52        _RL field(ib,nlevs), chfr(ib)        _RL field(ib,nlevs), chfr(ib)
53        logical check        LOGICAL check
54  CEOP  CEOP
55    
56  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
# Line 64  C =============== Line 60  C ===============
60        INTEGER k, kFirst, kLast        INTEGER k, kFirst, kLast
61        INTEGER kd, kd0, ksgn, kStore        INTEGER kd, kd0, ksgn, kStore
62        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
63        integer i,offset,Lena,newindx,jindx        INTEGER offset, Lena
64          INTEGER ivt, ij, i
65        _RL undef        _RL undef
66        INTEGER iSp, ndId, j,l        INTEGER iSp, ndId, j,l
67        INTEGER region2fill(0:nRegions)        INTEGER region2fill(0:nRegions)
68        _RL     scaleFact        _RL     scaleFact
69        _RL     gridField(sNx*sNy,nlevs), gridFrac(sNx*sNy)        _RL     gridField(sNx*sNy,nlevs), gridFrac(sNx*sNy)
70          _RS     dummyRS(1)
71    
72  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
73        _RL   getcon        _RL   getcon
# Line 80  C Run through list of active diagnostics Line 78  C Run through list of active diagnostics
78  C we are trying to fill a valid diagnostic  C we are trying to fill a valid diagnostic
79    
80        undef = UNSET_RL        undef = UNSET_RL
81    #ifdef ALLOW_FIZHI
82          IF ( check ) undef = getcon('UNDEF')
83    #endif
84        ndiagnum = 0        ndiagnum = 0
85        ipointer = 0        ipointer = 0
86        DO n=1,nlists        DO n=1,nlists
# Line 90  C we are trying to fill a valid diagnost Line 91  C we are trying to fill a valid diagnost
91           IF ( ndiagnum.NE.0 .AND. ndiag(ipointer,1,1).GE.0 ) THEN           IF ( ndiagnum.NE.0 .AND. ndiag(ipointer,1,1).GE.0 ) THEN
92  C--   do the filling: start here:  C--   do the filling: start here:
93    
94         IF ( (ABS(kLev).LE.1) .AND. (npeice.eq.1) ) THEN            IF ( (ABS(kLev).LE.1) .AND. (npeice.EQ.1) ) THEN
95  C Increment the counter for the diagnostic  C Increment the counter for the diagnostic
96           ndiag(ipointer,bi,bj) = ndiag(ipointer,bi,bj) + 1              ndiag(ipointer,bi,bj) = ndiag(ipointer,bi,bj) + 1
97         ENDIF            ENDIF
98    
99         offset = ib*(npeice-1)            offset = ib*(npeice-1)
100         Lena    = min(ib,numpts-offset)            Lena   = MIN(ib,numpts-offset)
        offset = offset+1  
101    
102  C-      Which part of field to add : k = 3rd index,  C-      Which part of field to add : k = 3rd index,
103  C         and do the loop >> do k=kFirst,kLast <<  C         and do the loop >> do k=kFirst,kLast <<
104         IF (kLev.LE.0) THEN            IF (kLev.LE.0) THEN
105          kFirst = 1             kFirst = 1
106          kLast  = nLevs             kLast  = nLevs
107         ELSEIF ( nLevs.EQ.1 ) THEN            ELSEIF ( nLevs.EQ.1 ) THEN
108          kFirst = 1             kFirst = 1
109          kLast  = 1             kLast  = 1
110         ELSEIF ( kLev.LE.nLevs ) THEN            ELSEIF ( kLev.LE.nLevs ) THEN
111          kFirst = kLev             kFirst = kLev
112          kLast  = kLev             kLast  = kLev
113         ELSE            ELSE
114          STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL kLev > nLevs > 0'             STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL kLev > nLevs > 0'
115         ENDIF            ENDIF
116  C-      Which part of qdiag to update: kd = 3rd index,  C-      Which part of qdiag to update: kd = 3rd index,
117  C         and do the loop >> do k=kFirst,kLast ; kd = kd0 + k*ksgn <<  C         and do the loop >> do k=kFirst,kLast ; kd = kd0 + k*ksgn <<
118         IF ( kLev.EQ.-1 ) THEN            IF ( kLev.EQ.-1 ) THEN
119          ksgn = -1             ksgn = -1
120          kd0 = ipointer + nLevs             kd0 = ipointer + nLevs
121         ELSEIF ( kLev.EQ.0 ) THEN            ELSEIF ( kLev.EQ.0 ) THEN
122          ksgn = 1             ksgn = 1
123          kd0 = ipointer - 1             kd0 = ipointer - 1
124         ELSE            ELSE
125          ksgn = 0             ksgn = 0
126          kd0 = ipointer + kLev - 1             kd0 = ipointer + kLev - 1
127         ENDIF            ENDIF
128    
129  C-      Check for consistency with Nb of levels reserved in storage array  C-      Check for consistency with Nb of levels reserved in storage array
130          kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - ipointer + 1            kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - ipointer + 1
131         IF ( kStore.GT.kdiag(ndiagnum) ) THEN            IF ( kStore.GT.kdiag(ndiagnum) ) THEN
132          _BEGIN_MASTER(myThid)             _BEGIN_MASTER(myThid)
133          WRITE(msgBuf,'(2A,I4,A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I4,A)') 'DIAGNOSTICS_FILL: ',
134       &     'exceed Nb of levels(=',kdiag(ndiagnum),' ) reserved '       &      'exceed Nb of levels(=',kdiag(ndiagnum),' ) reserved '
135          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
136          WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_FILL: ',
137       &     'for Diagnostics #', ndiagnum, ' : ', chardiag       &      'for Diagnostics #', ndiagnum, ' : ', chardiag
138          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
139          WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGNOSTICS_FILL ',             WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGNOSTICS_FILL ',
140       &     'with kLev,nLevs=', kLev,nLevs       &      'with kLev,nLevs=', kLev,nLevs
141          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
142          WRITE(msgBuf,'(2A,I6,A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I6,A)') 'DIAGNOSTICS_FILL: ',
143       &     '==> trying to store up to ', kStore, ' levels'       &      '==> trying to store up to ', kStore, ' levels'
144          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
145          STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL'             STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL'
146          _END_MASTER(myThid)             _END_MASTER(myThid)
147         ENDIF            ENDIF
148    
149         DO k = kFirst,kLast            DO k = kFirst,kLast
150          kd = kd0 + ksgn*k             kd = kd0 + ksgn*k
151          if( check ) then             IF ( check ) THEN
152  #ifdef ALLOW_FIZHI              DO ivt = 1,Lena
153           undef = getcon('UNDEF')               ij = indx(ivt+offset) - 1
154  #endif               j = 1 + INT(ij/sNx)
155           do i= 1,Lena               i = 1 + MOD(ij,sNx)
156            jindx = 1 + int((indx(i+offset-1)-1)/sNx)               IF ( field(ivt,k).EQ.undef ) THEN
157            newindx = indx(i+offset-1)+(jindx-1)*2*Olx                qdiag(i,j,kd,bi,bj) = undef
158            if(qdiag(newindx,1,kd,bi,bj).eq.undef               ELSEIF ( qdiag(i,j,kd,bi,bj).NE.undef ) THEN
159       .                                  .or.field(i,k).eq.undef)then                qdiag(i,j,kd,bi,bj) = qdiag(i,j,kd,bi,bj)
160             qdiag(newindx,1,kd,bi,bj) = undef       &                            + field(ivt,k)*chfr(ivt)
161            else               ENDIF
162             qdiag(newindx,1,kd,bi,bj)=qdiag(newindx,1,kd,bi,bj)+              ENDDO
163       .                                                field(i,k)*chfr(i)             ELSE
164            endif              DO ivt = 1,Lena
165           enddo                ij = indx(ivt+offset) - 1
166          else                j = 1 + INT(ij/sNx)
167           do i= 1,Lena                i = 1 + MOD(ij,sNx)
168            jindx = 1 + int((indx(i+offset-1)-1)/sNx)                qdiag(i,j,kd,bi,bj) = qdiag(i,j,kd,bi,bj)
169            newindx = indx(i+offset-1)+(jindx-1)*2*Olx       &                            + field(ivt,k)*chfr(ivt)
170            qdiag(newindx,1,kd,bi,bj)=qdiag(newindx,1,kd,bi,bj)+              ENDDO
171       .                                                field(i,k)*chfr(i)             ENDIF
172           enddo            ENDDO
         endif  
        ENDDO  
173    
174  C--   do the filling: ends here.  C--   do the filling: ends here.
175           ENDIF           ENDIF
# Line 210  C     then add regions from other lists Line 208  C     then add regions from other lists
208    
209  C-      Which part of field to add : k = 3rd index,  C-      Which part of field to add : k = 3rd index,
210  C         and do the loop >> do k=kFirst,kLast <<  C         and do the loop >> do k=kFirst,kLast <<
211         IF (kLev.LE.0) THEN             IF (kLev.LE.0) THEN
212          kFirst = 1              kFirst = 1
213          kLast  = nLevs              kLast  = nLevs
214         ELSE             ELSE
215          kFirst = 1              kFirst = 1
216          kLast  = 1              kLast  = 1
217         ENDIF             ENDIF
218    
219  C-    Fill local array with grid-space field after conversion.  C-    Fill local array with grid-space field after conversion.
220         offset = ib*(npeice-1)             offset = ib*(npeice-1)
221         Lena    = min(ib,numpts-offset)             Lena    = MIN(ib,numpts-offset)
        offset = offset+1  
222    
223         DO i=1,sNx*sNy             DO ij = 1,sNx*sNy
224           gridFrac(i)= 0.               gridFrac(ij)= 0.
225         ENDDO             ENDDO
226         DO i=1,Lena             DO ivt = 1,Lena
227           newindx = indx(i+offset-1)               ij = indx(ivt+offset)
228           gridFrac(newindx)=gridFrac(newindx)+chfr(i)               gridFrac(ij)=gridFrac(ij)+chfr(ivt)
229         ENDDO             ENDDO
230    
231         DO k = kFirst,kLast             DO k = kFirst,kLast
232          DO i=1,sNx*sNy              DO ij = 1,sNx*sNy
233           gridField(i,k)= 0.               gridField(ij,k)= 0.
234          ENDDO              ENDDO
235          if( check ) then              IF ( check ) THEN
236  #ifdef ALLOW_FIZHI               DO ivt = 1,Lena
237           undef = getcon('UNDEF')                ij = indx(ivt+offset)
238  #endif                IF ( field(ivt,k).EQ.undef ) THEN
239           do i= 1,Lena                 gridField(ij,k) = undef
240            newindx = indx(i+offset-1)                ELSEIF ( gridFrac(ij).GT.0. _d 0 ) THEN
241            if(gridField(newindx,k).eq.undef                 gridField(ij,k) = gridField(ij,k)
242       .                           .or.field(i,k).eq.undef)then       &                         + field(ivt,k)*chfr(ivt)/gridFrac(ij)
243             gridField(newindx,k)= undef                ENDIF
244            else               ENDDO
245             gridField(newindx,k)=gridField(newindx,k)              ELSE
246       &                         +field(i,k)*chfr(i)/gridFrac(newindx)               DO ivt = 1,Lena
247            endif                ij = indx(ivt+offset)
248           enddo                IF ( gridFrac(ij).GT.0. _d 0 ) THEN
249          else                 gridField(ij,k) = gridField(ij,k)
250           do i= 1,Lena       &                         + field(ivt,k)*chfr(ivt)/gridFrac(ij)
251            newindx = indx(i+offset-1)                ENDIF
252            gridField(newindx,k)=gridField(newindx,k)               ENDDO
253       &                        +field(i,k)*chfr(i)/gridFrac(newindx)              ENDIF
254           enddo             ENDDO
         endif  
        ENDDO  
255    
256  C-    diagnostics is valid and Active: Now do the filling  C-    diagnostics is valid and Active: Now do the filling
257             CALL DIAGSTATS_FILL(             CALL DIAGSTATS_FILL(
258       I              gridField, gridFrac, scaleFact, 1, 1,       I              gridField, gridFrac,
259    #ifndef REAL4_IS_SLOW
260         I               dummyRS, dummyRS,
261    #endif
262         I              scaleFact, 1, 0, 1,
263       I              ndId, iSp, region2fill, kLev, nLevs,       I              ndId, iSp, region2fill, kLev, nLevs,
264       I              3, bi, bj, myThid )       I              3, bi, bj, myThid )
265           ENDIF           ENDIF
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22