/[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.4 by molod, Fri Feb 18 19:44:11 2005 UTC revision 1.12 by jmc, Sun Jan 3 00:42:45 2010 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:
57  C ===============  C ===============
58        INTEGER m, n        INTEGER m, n
59        INTEGER ndiagnum, ipointer        INTEGER ndiagnum, ipointer
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        _RL undef,getcon        INTEGER ivt, ij, i
65          _RL undef
66          INTEGER iSp, ndId, j,l
67          INTEGER region2fill(0:nRegions)
68          _RL     scaleFact
69          _RL     gridField(sNx*sNy,nlevs), gridFrac(sNx*sNy)
70    #ifndef REAL4_IS_SLOW
71          _RS     dummyRS(1)
72    #endif
73    
74    #ifdef ALLOW_FIZHI
75          _RL   getcon
76          EXTERNAL getcon
77    #endif
78    
79  C Run through list of active diagnostics to make sure  C Run through list of active diagnostics to make sure
80  C we are trying to fill a valid diagnostic  C we are trying to fill a valid diagnostic
81    
82          undef = UNSET_RL
83    #ifdef ALLOW_FIZHI
84          IF ( check ) undef = getcon('UNDEF')
85    #endif
86        ndiagnum = 0        ndiagnum = 0
87        ipointer = 0        ipointer = 0
88        DO n=1,nlists        DO n=1,nlists
89         DO m=1,nActive(n)         DO m=1,nActive(n)
90          IF ( chardiag.EQ.flds(m,n) ) THEN          IF ( chardiag.EQ.flds(m,n) .AND. idiag(m,n).GT.0 ) THEN
91           ndiagnum = jdiag(m,n)           ndiagnum = jdiag(m,n)
92           IF (ndiag(ndiagnum).GE.0) ipointer = idiag(ndiagnum)           ipointer = idiag(m,n)
93          ENDIF           IF ( ndiagnum.NE.0 .AND. ndiag(ipointer,1,1).GE.0 ) THEN
94         ENDDO  C--   do the filling: start here:
95        ENDDO  
96              IF ( (ABS(kLev).LE.1) .AND. (npeice.EQ.1) ) THEN
97        print *,' in diag_vegtile_fill for ',chardiag,  C Increment the counter for the diagnostic
98       .     ' ndiagnum ',ndiagnum,' ipointer ',ipointer              ndiag(ipointer,bi,bj) = ndiag(ipointer,bi,bj) + 1
99              ENDIF
100    
101  C If-sequence to see if we are a valid and an active diagnostic            offset = ib*(npeice-1)
102              Lena   = MIN(ib,numpts-offset)
       IF ( ndiagnum.NE.0 .AND. ipointer.NE.0 ) THEN  
   
 C Increment the counter for the diagnostic (if we are at bi=bj=myThid=1)  
        _BEGIN_MASTER(myThid)  
        IF((bi.EQ.1).AND.(bj.EQ.1).AND.(ABS(kLev).LE.1)  
      .         .AND.(npeice.eq.1))  
      .                     ndiag(ndiagnum) = ndiag(ndiagnum) + 1  
        _END_MASTER(myThid)  
   
        offset = ib*(npeice-1)  
        Lena    = min(ib,numpts-offset)  
        offset = offset+1  
103    
104  C-      Which part of field to add : k = 3rd index,  C-      Which part of field to add : k = 3rd index,
105  C         and do the loop >> do k=kFirst,kLast <<  C         and do the loop >> do k=kFirst,kLast <<
106         IF (kLev.LE.0) THEN            IF (kLev.LE.0) THEN
107          kFirst = 1             kFirst = 1
108          kLast  = nLevs             kLast  = nLevs
109         ELSEIF ( nLevs.EQ.1 ) THEN            ELSEIF ( nLevs.EQ.1 ) THEN
110          kFirst = 1             kFirst = 1
111          kLast  = 1             kLast  = 1
112         ELSEIF ( kLev.LE.nLevs ) THEN            ELSEIF ( kLev.LE.nLevs ) THEN
113          kFirst = kLev             kFirst = kLev
114          kLast  = kLev             kLast  = kLev
115         ELSE            ELSE
116          STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL kLev > nLevs > 0'             STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL kLev > nLevs > 0'
117         ENDIF            ENDIF
118  C-      Which part of qdiag to update: kd = 3rd index,  C-      Which part of qdiag to update: kd = 3rd index,
119  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 <<
120         IF ( kLev.EQ.-1 ) THEN            IF ( kLev.EQ.-1 ) THEN
121          ksgn = -1             ksgn = -1
122          kd0 = ipointer + nLevs             kd0 = ipointer + nLevs
123         ELSEIF ( kLev.EQ.0 ) THEN            ELSEIF ( kLev.EQ.0 ) THEN
124          ksgn = 1             ksgn = 1
125          kd0 = ipointer - 1             kd0 = ipointer - 1
126         ELSE            ELSE
127          ksgn = 0             ksgn = 0
128          kd0 = ipointer + kLev - 1             kd0 = ipointer + kLev - 1
129         ENDIF            ENDIF
130    
131  C-      Check for consistency with Nb of levels reserved in storage array  C-      Check for consistency with Nb of levels reserved in storage array
132          kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - ipointer + 1            kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - ipointer + 1
133         IF ( kStore.GT.kdiag(ndiagnum) ) THEN            IF ( kStore.GT.kdiag(ndiagnum) ) THEN
134          _BEGIN_MASTER(myThid)             _BEGIN_MASTER(myThid)
135          WRITE(msgBuf,'(2A,I3,A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I4,A)') 'DIAGNOSTICS_FILL: ',
136       .     'exceed Nb of levels(=',kdiag(ndiagnum),' ) reserved '       &      'exceed Nb of levels(=',kdiag(ndiagnum),' ) reserved '
137          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
138          WRITE(msgBuf,'(2A,I4,2A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_FILL: ',
139       .     'for Diagnostics #', ndiagnum, ' : ', chardiag       &      'for Diagnostics #', ndiagnum, ' : ', chardiag
140          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
141          WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGNOSTICS_FILL ',             WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGNOSTICS_FILL ',
142       .     'with kLev,nLevs=', kLev,nLevs       &      'with kLev,nLevs=', kLev,nLevs
143          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
144          WRITE(msgBuf,'(2A,I6,A)') 'DIAGNOSTICS_FILL: ',             WRITE(msgBuf,'(2A,I6,A)') 'DIAGNOSTICS_FILL: ',
145       .     '==> trying to store up to ', kStore, ' levels'       &      '==> trying to store up to ', kStore, ' levels'
146          CALL PRINT_ERROR( msgBuf , myThid )             CALL PRINT_ERROR( msgBuf , myThid )
147          STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL'             STOP 'ABNORMAL END: S/R DIAGNOSTICS_FILL'
148          _END_MASTER(myThid)             _END_MASTER(myThid)
149         ENDIF            ENDIF
150    
151         DO k = kFirst,kLast            DO k = kFirst,kLast
152          kd = kd0 + ksgn*k             kd = kd0 + ksgn*k
153          print *,' level slot in qdiag= ',kd, ' level from input= ',k             IF ( check ) THEN
154          if( check ) then              DO ivt = 1,Lena
155           undef = getcon('UNDEF')               ij = indx(ivt+offset) - 1
156           do i= 1,Lena               j = 1 + INT(ij/sNx)
157            jindx = 1 + int((indx(i+offset-1)-1)/sNx)               i = 1 + MOD(ij,sNx)
158            newindx = indx(i+offset-1)+(jindx-1)*2*Olx               IF ( field(ivt,k).EQ.undef ) THEN
159            if(qdiag(newindx,1,kd,bi,bj).eq.undef                qdiag(i,j,kd,bi,bj) = undef
160       .                                  .or.field(i,k).eq.undef)then               ELSEIF ( qdiag(i,j,kd,bi,bj).NE.undef ) THEN
161             qdiag(newindx,1,kd,bi,bj) = undef                qdiag(i,j,kd,bi,bj) = qdiag(i,j,kd,bi,bj)
162            else       &                            + field(ivt,k)*chfr(ivt)
163             qdiag(newindx,1,kd,bi,bj)=qdiag(newindx,1,kd,bi,bj)+               ENDIF
164       .                                                field(i,k)*chfr(i)              ENDDO
165            endif             ELSE
166           enddo              DO ivt = 1,Lena
167          else                ij = indx(ivt+offset) - 1
168           do i= 1,Lena                j = 1 + INT(ij/sNx)
169            jindx = 1 + int((indx(i+offset-1)-1)/sNx)                i = 1 + MOD(ij,sNx)
170            newindx = indx(i+offset-1)+(jindx-1)*2*Olx                qdiag(i,j,kd,bi,bj) = qdiag(i,j,kd,bi,bj)
171            qdiag(newindx,1,kd,bi,bj)=qdiag(newindx,1,kd,bi,bj)+       &                            + field(ivt,k)*chfr(ivt)
172       .                                                field(i,k)*chfr(i)              ENDDO
173           enddo             ENDIF
174          endif            ENDDO
175    
176    C--   do the filling: ends here.
177             ENDIF
178            ENDIF
179         ENDDO         ENDDO
180          ENDDO
181    
182    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
183    C--   Global/Regional Statistics :
184          scaleFact = 1. _d 0
185    
186    C Run through list of active statistics-diagnostics to make sure
187    C we are trying to compute & fill a valid diagnostic
188    
189          DO n=1,diagSt_nbLists
190           DO m=1,diagSt_nbActv(n)
191            IF ( chardiag.EQ.diagSt_Flds(m,n) .AND. iSdiag(m,n).GT.0 ) THEN
192             iSp = iSdiag(m,n)
193             IF ( qSdiag(0,0,iSp,bi,bj).GE.0. ) THEN
194               ndId = jSdiag(m,n)
195    C-    Find list of regions to fill:
196               DO j=0,nRegions
197                region2fill(j) = diagSt_region(j,n)
198               ENDDO
199    C-    if this diagnostics appears in several lists (with same freq)
200    C     then add regions from other lists
201               DO l=1,diagSt_nbLists
202                DO k=1,diagSt_nbActv(l)
203                 IF ( iSdiag(k,l).EQ.-iSp ) THEN
204                  DO j=0,nRegions
205                   region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))
206                  ENDDO
207                 ENDIF
208                ENDDO
209               ENDDO
210    
211    C-      Which part of field to add : k = 3rd index,
212    C         and do the loop >> do k=kFirst,kLast <<
213               IF (kLev.LE.0) THEN
214                kFirst = 1
215                kLast  = nLevs
216               ELSE
217                kFirst = 1
218                kLast  = 1
219               ENDIF
220    
221    C-    Fill local array with grid-space field after conversion.
222               offset = ib*(npeice-1)
223               Lena    = MIN(ib,numpts-offset)
224    
225               DO ij = 1,sNx*sNy
226                 gridFrac(ij)= 0.
227               ENDDO
228               DO ivt = 1,Lena
229                 ij = indx(ivt+offset)
230                 gridFrac(ij)=gridFrac(ij)+chfr(ivt)
231               ENDDO
232    
233               DO k = kFirst,kLast
234                DO ij = 1,sNx*sNy
235                 gridField(ij,k)= 0.
236                ENDDO
237                IF ( check ) THEN
238                 DO ivt = 1,Lena
239                  ij = indx(ivt+offset)
240                  IF ( field(ivt,k).EQ.undef ) THEN
241                   gridField(ij,k) = undef
242                  ELSEIF ( gridFrac(ij).GT.0. _d 0 ) THEN
243                   gridField(ij,k) = gridField(ij,k)
244         &                         + field(ivt,k)*chfr(ivt)/gridFrac(ij)
245                  ENDIF
246                 ENDDO
247                ELSE
248                 DO ivt = 1,Lena
249                  ij = indx(ivt+offset)
250                  IF ( gridFrac(ij).GT.0. _d 0 ) THEN
251                   gridField(ij,k) = gridField(ij,k)
252         &                         + field(ivt,k)*chfr(ivt)/gridFrac(ij)
253                  ENDIF
254                 ENDDO
255                ENDIF
256               ENDDO
257    
258    C-    diagnostics is valid and Active: Now do the filling
259               CALL DIAGSTATS_FILL(
260         I              gridField, gridFrac,
261    #ifndef REAL4_IS_SLOW
262         I               dummyRS, dummyRS,
263    #endif
264         I              scaleFact, 1, 0, 1,
265         I              ndId, iSp, region2fill, kLev, nLevs,
266         I              3, bi, bj, myThid )
267             ENDIF
268            ENDIF
269           ENDDO
270          ENDDO
271    
272        ELSE        RETURN
 C     IF (myThid.EQ.1) WRITE(6,1000) chardiag  
       ENDIF  
   
  1000 format(' ',' Warning: Trying to write to diagnostic ',a8,  
      .        ' But it is not a valid (or active) name ')  
       RETURN  
273        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22