- Timestamp:
- 2022-02-24T11:04:47+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
r6189 r7500 3708 3708 i = (index_g(ik)-(j-1)*iim_g) 3709 3709 3710 basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,1)*resolution_g(ik,2)) 3710 IF ( resolution_g(ik,1) == 0 ) THEN 3711 basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,2)*resolution_g(ik,2)*pi) 3712 ELSE 3713 basinfrac(i,j) = basinfrac(i,j) + lbasin_area(ib,ij)/(resolution_g(ik,1)*resolution_g(ik,2)) 3714 ENDIF 3711 3715 basinuparea(i,j) = MAX(basinuparea(i,j), lbasin_uparea(ib,ij)) 3712 3716 basincode(i,j) = lrivercode(ib,ij) … … 4708 4712 basin_bx(iloc, jloc) = NINT(max_basins + 1) 4709 4713 max_basins = max_basins + 1 4710 area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib) 4714 ! Check if we are at the poles : resolution(ib,1) = 0 4715 IF ( resolution(ib,1) == 0 ) THEN 4716 ! compute the pole cell area as the circle surface 4717 area_bx(iloc, jloc) = pi*resolution(ib,2)*resolution(ib,2)*contfrac(ib) 4718 ELSE 4719 area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib) 4720 ENDIF 4711 4721 topoind_bx(iloc, jloc) = min_topoind 4712 4722 hierarchy_bx(iloc, jloc) = min_topoind … … 6519 6529 ! 6520 6530 totbasins = SUM(basin_area(ib,1:basin_count(ib))) 6521 contarea = resolution(ib,1)*resolution(ib,2)*contfrac(ib) 6531 ! Check if we are at the poles (resolution(ib,1) = 0 6532 IF ( resolution(ib,1) == 0 ) THEN 6533 ! Hack to approximate the pole cell area by a circle 6534 contarea = pi*resolution(ib,2)*resolution(ib,2)*contfrac(ib) 6535 ELSE 6536 contarea = resolution(ib,1)*resolution(ib,2)*contfrac(ib) 6537 ENDIF 6522 6538 ! 6523 6539 DO ij=1,basin_count(ib) … … 7034 7050 ! 7035 7051 floflo(:,:) = zero 7036 gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2) 7052 ! if we are at the poles : resolution(:,1) = 0 7053 WHERE (resolution(:,1) == 0) 7054 ! compute grid area as the circle of radius resolution(:,2) 7055 gridarea(:) = contfrac(:)*pi*resolution(:,2)*resolution(:,2) 7056 ELSEWHERE 7057 gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2) 7058 END WHERE 7037 7059 DO ib=1,nbpt 7038 7060 gridbasinarea(ib) = SUM(routing_area(ib,:)) … … 7761 7783 REAL(r_std) :: area_irrig !! Irrigated surface in the grid box (m^2) 7762 7784 REAL(r_std) :: area_flood(ntype) !! Flooded surface in the grid box (m^2) 7785 REAL(r_std) :: resolution_1 !! temporary variable 7763 7786 !!$ REAL(r_std) :: irrigmap(nbpt) 7764 7787 !!$ REAL(r_std) :: floodmap(nbpt) … … 7958 7981 ! 7959 7982 IF ( init_irrig ) THEN 7960 irrigated(ib) = MIN(area_irrig, resolution(ib,1)*resolution(ib,2)*contfrac(ib)) 7983 ! if we are at the poles resolution(ib,1) = 0 7984 IF (resolution(ib,1) == 0) THEN 7985 ! use pi*resolution(ib,2) to get the disc area 7986 resolution_1 = pi*resolution(ib,2) 7987 ELSE 7988 resolution_1 = resolution(ib,1) 7989 END IF 7990 irrigated(ib) = MIN(area_irrig, resolution_1*resolution(ib,2)*contfrac(ib)) 7961 7991 IF ( irrigated(ib) < 0 ) THEN 7962 7992 WRITE(numout,*) 'We have a problem here : ', irrigated(ib) 7963 WRITE(numout,*) 'resolution :', resolution (ib,1), resolution(ib,2)7993 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 7964 7994 WRITE(numout,*) area_irrig 7965 7995 CALL ipslerr_p(3,'routing_irrigmap','Problem with irrigated...','','') … … 7975 8005 ! 7976 8006 IF ( init_flood ) THEN 8007 ! if we are at the poles resolution(ib,1) = 0 8008 IF (resolution(ib,1) == 0) THEN 8009 ! use pi*resolution(ib,2) to get the disc area 8010 resolution_1 = pi*resolution(ib,2) 8011 ELSE 8012 resolution_1 = resolution(ib,1) 8013 END IF 7977 8014 floodplains(ib) = MIN(area_flood(iflood)+area_flood(idam)+area_flood(isal), & 7978 & resolution (ib,1)*resolution(ib,2)*contfrac(ib))8015 & resolution_1*resolution(ib,2)*contfrac(ib)) 7979 8016 IF ( floodplains(ib) < 0 ) THEN 7980 8017 WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 7981 WRITE(numout,*) 'resolution :', resolution (ib,1), resolution(ib,2)8018 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 7982 8019 WRITE(numout,*) area_flood 7983 8020 CALL ipslerr_p(3,'routing_irrigmap','Problem with floodplains..','','') … … 7992 8029 ! 7993 8030 IF ( init_swamp ) THEN 7994 swamp(ib) = MIN(area_flood(iswamp), resolution(ib,1)*resolution(ib,2)*contfrac(ib)) 8031 ! if we are at the poles resolution(ib,1) = 0 8032 IF (resolution(ib,1) == 0) THEN 8033 ! use pi*resolution(ib,2) to get the disc area 8034 resolution_1 = pi*resolution(ib,2) 8035 ELSE 8036 resolution_1 = resolution(ib,1) 8037 END IF 8038 swamp(ib) = MIN(area_flood(iswamp), resolution_1*resolution(ib,2)*contfrac(ib)) 7995 8039 IF ( swamp(ib) < 0 ) THEN 7996 8040 WRITE(numout,*) 'We have a problem here : ', swamp(ib) 7997 WRITE(numout,*) 'resolution :', resolution (ib,1), resolution(ib,2)8041 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 7998 8042 WRITE(numout,*) area_flood 7999 8043 CALL ipslerr_p(3,'routing_irrigmap','Problem with swamp...','','')
Note: See TracChangeset
for help on using the changeset viewer.