Changeset 7500 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba

Timestamp:

2022-02-24T11:04:47+01:00 (2 years ago)

Author:

josefine.ghattas

Message:

Correction to allow routing scheme to work on the south pole: Check if resoultion(ib,1)=0 and if so, calculate the area in a different way.
Solution done by Sebastien Nguyen. See ticket #811

File:

• branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90 (modified) (8 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90

@@ -3708 +3708 @@
              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 +4712 @@
        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 +6529 @@
        !
        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 +7050 @@
     !
     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 +7783 @@
     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 +7981 @@
        !
        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 +8005 @@
        !
        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 +8029 @@
        !
        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...','','')