Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 7493)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 7500)
@@ -3708,5 +3708,9 @@
              i = (index_g(ik)-(j-1)*iim_g)
 
-             basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,1)*resolution_g(ik,2))
+             IF ( resolution_g(ik,1) == 0 ) THEN
+                basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,2)*resolution_g(ik,2)*pi)
+             ELSE
+                basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,1)*resolution_g(ik,2))
+             ENDIF
              basinuparea(i,j) = MAX(basinuparea(i,j), lbasin_uparea(ib,ij))
              basincode(i,j) = lrivercode(ib,ij)
@@ -4708,5 +4712,11 @@
        basin_bx(iloc, jloc) = NINT(max_basins + 1)
        max_basins = max_basins + 1
-       area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib)
+       ! Check if we are at the poles : resolution(ib,1) = 0
+       IF ( resolution(ib,1) == 0 ) THEN
+          ! compute the pole cell area as the circle surface
+          area_bx(iloc, jloc) = pi*resolution(ib,2)*resolution(ib,2)*contfrac(ib)
+       ELSE
+          area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib)
+       ENDIF
        topoind_bx(iloc, jloc) = min_topoind
        hierarchy_bx(iloc, jloc) =  min_topoind
@@ -6519,5 +6529,11 @@
        !
        totbasins = SUM(basin_area(ib,1:basin_count(ib)))
-       contarea = resolution(ib,1)*resolution(ib,2)*contfrac(ib)
+       ! Check if we are at the poles (resolution(ib,1) = 0
+       IF ( resolution(ib,1) == 0 ) THEN
+          ! Hack to approximate the pole cell area by a circle
+          contarea = pi*resolution(ib,2)*resolution(ib,2)*contfrac(ib)
+       ELSE
+          contarea = resolution(ib,1)*resolution(ib,2)*contfrac(ib)
+       ENDIF
        !
        DO ij=1,basin_count(ib)
@@ -7034,5 +7050,11 @@
     !
     floflo(:,:) = zero
-    gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2)
+    ! if we are at the poles : resolution(:,1) = 0
+    WHERE (resolution(:,1) == 0)
+        ! compute grid area as the circle of radius resolution(:,2)
+        gridarea(:) = contfrac(:)*pi*resolution(:,2)*resolution(:,2)
+    ELSEWHERE
+        gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2)
+    END WHERE
     DO ib=1,nbpt
        gridbasinarea(ib) = SUM(routing_area(ib,:))
@@ -7761,4 +7783,5 @@
     REAL(r_std)                                    :: area_irrig            !! Irrigated surface in the grid box (m^2)
     REAL(r_std)                                    :: area_flood(ntype)     !! Flooded surface in the grid box (m^2)
+    REAL(r_std)                                    :: resolution_1          !! temporary variable
 !!$    REAL(r_std)                                :: irrigmap(nbpt)
 !!$    REAL(r_std)                                :: floodmap(nbpt)
@@ -7958,8 +7981,15 @@
        !
        IF ( init_irrig ) THEN
-          irrigated(ib) = MIN(area_irrig, resolution(ib,1)*resolution(ib,2)*contfrac(ib))
+          ! if we are at the poles resolution(ib,1) = 0
+          IF (resolution(ib,1) == 0) THEN
+             ! use pi*resolution(ib,2) to get the disc area
+             resolution_1 = pi*resolution(ib,2)
+          ELSE
+             resolution_1 = resolution(ib,1)
+          END IF
+          irrigated(ib) = MIN(area_irrig, resolution_1*resolution(ib,2)*contfrac(ib))
           IF ( irrigated(ib) < 0 ) THEN
              WRITE(numout,*) 'We have a problem here : ', irrigated(ib) 
-             WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
+             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
              WRITE(numout,*) area_irrig
              CALL ipslerr_p(3,'routing_irrigmap','Problem with irrigated...','','')
@@ -7975,9 +8005,16 @@
        !
        IF ( init_flood ) THEN
+          ! if we are at the poles resolution(ib,1) = 0
+          IF (resolution(ib,1) == 0) THEN 
+             ! use pi*resolution(ib,2) to get the disc area
+             resolution_1 = pi*resolution(ib,2)
+          ELSE
+             resolution_1 = resolution(ib,1)
+          END IF
           floodplains(ib) = MIN(area_flood(iflood)+area_flood(idam)+area_flood(isal), &
-               & resolution(ib,1)*resolution(ib,2)*contfrac(ib))
+               & resolution_1*resolution(ib,2)*contfrac(ib))
           IF ( floodplains(ib) < 0 ) THEN
              WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 
-             WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
+             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
              WRITE(numout,*) area_flood
              CALL ipslerr_p(3,'routing_irrigmap','Problem with floodplains..','','')
@@ -7992,8 +8029,15 @@
        !
        IF ( init_swamp ) THEN
-          swamp(ib) = MIN(area_flood(iswamp), resolution(ib,1)*resolution(ib,2)*contfrac(ib))
+          ! if we are at the poles resolution(ib,1) = 0
+          IF (resolution(ib,1) == 0) THEN 
+             ! use pi*resolution(ib,2) to get the disc area 
+             resolution_1 = pi*resolution(ib,2)
+          ELSE
+             resolution_1 = resolution(ib,1)
+          END IF
+          swamp(ib) = MIN(area_flood(iswamp), resolution_1*resolution(ib,2)*contfrac(ib))
           IF ( swamp(ib) < 0 ) THEN
              WRITE(numout,*) 'We have a problem here : ', swamp(ib) 
-             WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
+             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
              WRITE(numout,*) area_flood
              CALL ipslerr_p(3,'routing_irrigmap','Problem with swamp...','','')
