source: CONFIG/trunk/SCRIPT/REGRID_forcage/ETAT0/grid_correspondance.f @ 1356

      SUBROUTINE grid_correspondance
!--------------------------------------------------------------
!     Determine the grid correspondance between initial grid
!     and final grid and calculates the weighting function based
!     on the corresponding grid-cell area.
!
!     It is assumed that the initial grid is defined as:
!
!         ---------- lat_edge(j+1)
!         | (i,j) |
!         |   +   |
!         |       |
!         ---------- lat_edge(j)
!         |       |
!    lon_edge(i) lon_edge(i+1)
!
!--------------------------------------------------------------
      USE final_grid_lmdz
      USE grid
      USE correspondance
      IMPLICIT NONE

!--------------------------------------------------------------
!       ... Local variables
!--------------------------------------------------------------
      INTEGER :: error
      INTEGER :: i, j
      INTEGER :: oldi, oldj
      INTEGER :: from_lon, to_lon 
      INTEGER :: from_lat, to_lat
      INTEGER, SAVE :: countmax = 200000

      REAL    :: lonw, lone, lats, latn
      REAL    :: d2r, pi
      REAL    :: area

      pi      = 4.*ATAN(1.)
      d2r     = pi/180.

      PRINT *, 'Grid correspondance ...'

      ! ... Loop over countmax if needed
      outer : DO

      ! ... Free space if needed and allocate correspondance array
        IF ( ALLOCATED(corr) )  DEALLOCATE(corr)
        ALLOCATE(corr(countmax), STAT=error)
        IF (error /= 0) STOP 'Space requested not possible'
        
        count = 0

      ! ... loop over new longitudes 
!      print*,'nlon_len',nlon_len
!      print*,'lon_len',lon_len
!      print*,'nlon_edge',nlon_edge
!      print*,'lon_edge',lon_edge
!      print*,'nlat_len',nlat_len
!      print*,'lat_len',lat_len
!      print*,'nlat_edge',nlat_edge
!      print*,'lat_edge',lat_edge
        DO i = 1, nlon_len             
      
          DO oldi = 1, lon_len
            from_lon = oldi
            IF (nlon_edge(i) <= lon_edge(oldi+1)) EXIT   
          END DO
          DO oldi = from_lon, lon_len
            to_lon = oldi
            IF (nlon_edge(i+1) <= lon_edge(oldi+1)) EXIT   
          END DO
          
      ! ... loop over new latitudes
          DO j = 1, nlat_len
      
            DO oldj = 1, lat_len
              from_lat = oldj
              IF (nlat_edge(j+1) > lat_edge(oldj) .AND. nlat_edge(j) < lat_edge(oldj+1) ) EXIT
              from_lat = 0
            END DO
            IF ( from_lat == 0 ) CYCLE
            DO oldj = from_lat, lat_len
              IF (nlat_edge(j+1) > lat_edge(oldj) .AND. nlat_edge(j) < lat_edge(oldj+1) ) THEN
                  to_lat = oldj
                  CYCLE
              END IF
              EXIT
            END DO
     
      ! ... loop over old longitudes and latitudes
            DO oldi = from_lon, to_lon
              DO oldj = from_lat, to_lat

                lonw = MAX(nlon_edge(i),lon_edge(oldi)) 
                lone = MIN(nlon_edge(i+1),lon_edge(oldi+1))
                lats = MAX(nlat_edge(j),lat_edge(oldj))
                latn = MIN(nlat_edge(j+1),lat_edge(oldj+1))
                
      ! Validity checks
                IF ( (sin(latn*d2r)-sin(lats*d2r)) <0) THEN
                    print *, 'WARNING: latn < lats ', latn, lats
                END IF
                IF ( lone < lonw ) THEN
                    print *, 'WARNING: lone < lonw ', lone, lonw
                END IF
                
                count = count + 1
      ! ... If needed increase countmax, reset, and restart
                IF (count > countmax) THEN
                    PRINT *, 'count (=', count, ') larger than countmax (=', countmax, ')'
                    PRINT *, 'Increase countmax to ', countmax*2, ' and restart'
                    countmax = countmax * 2
                    CYCLE outer
                END IF

                area = (sin(latn*d2r)-sin(lats*d2r))*(lone-lonw)*d2r
                
                corr(count)%olon = oldi
                corr(count)%olat = oldj
                corr(count)%nlon = i
                corr(count)%nlat = j
                corr(count)%wght = area / area_out(i,j)
                
              END DO
            END DO
            
          END DO 
          
        END DO

        EXIT outer
      END DO outer

!      PRINT *, '     OK.'

      END SUBROUTINE grid_correspondance