Changeset 277

Timestamp:

2011-06-23T11:25:25+02:00 (14 years ago)

Author:

didier.solyga

Message:

Update the externalized version with the last commit of the trunk (revision 275)

Location:

branches/ORCHIDEE_EXT/ORCHIDEE

Files:

• src_sechiba/intersurf.f90 (modified) (19 diffs)
• src_sechiba/routing.f90 (modified) (61 diffs)
• src_sechiba/sechiba.f90 (modified) (6 diffs)
• src_sechiba/slowproc.f90 (modified) (8 diffs)
• src_stomate/stomate.f90 (modified) (4 diffs)

branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/intersurf.f90

@@ -7 +7 @@
 !!
 !! @call sechiba_main
-!! @Version : $Revision: 221 $, $Date: 2011-05-16 17:26:17 +0200 (Mon, 16 May 2011) $
+!! @Version : $Revision: 275 $, $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
 !!
 !! @author Marie-Alice Foujols and Jan Polcher
 !! 
 !< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/intersurf.f90 $
-!< $Date: 2011-05-16 17:26:17 +0200 (Mon, 16 May 2011) $
+!< $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
 !< $Author: martial.mancip $
-!< $Revision: 221 $
+!< $Revision: 275 $
 !! IPSL (2006)
 !!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
@@ -162 +162 @@
     REAL(r_std),DIMENSION (kjpindex)                      :: dcoastal      !! Work array to keep coastalflow
     REAL(r_std),DIMENSION (kjpindex)                      :: driver        !! Work array to keep riverflow
+    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux
+    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use
     REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
     REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
@@ -357 +359 @@
        & zprecip_rain ,zprecip_snow,  zlwdown, zswnet, zswdown, zpb, &
 ! Output : Fluxes
-       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
+       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
 ! Surface temperatures and surface properties
        & ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0, &
@@ -701 +703 @@
     REAL(r_std),DIMENSION (kjpindex)                      :: dcoastal      !! Work array to keep coastal flow
     REAL(r_std),DIMENSION (kjpindex)                      :: driver        !! Work array to keep river out flow
+    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux
+    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use
     REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
     REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
@@ -874 +878 @@
        & zprecip_rain ,zprecip_snow,  zlwdown, zswnet, zswdown, zpb, &
 ! Output : Fluxes
-       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
+       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
 ! Surface temperatures and surface properties
        & ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0, &
@@ -1211 +1215 @@
     REAL(r_std),DIMENSION (kjpindex)                      :: dcoastal      !! Work array to keep coastal flow
     REAL(r_std),DIMENSION (kjpindex)                      :: driver        !! Work array to keep river out flow
+    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux
+    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use
     REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
     REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
@@ -1571 +1577 @@
        & zprecip_rain ,zprecip_snow,  lwdown, swnet, swdown, pb, &
 ! Output : Fluxes
-       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
+       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
 ! Surface temperatures and surface properties
        & ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0, &
@@ -1850 +1856 @@
      & tsol_rad, temp_sol_new, qsurf, albedo, emis, z0, lon_scat_g, lat_scat_g, &
 ! Ajout Nathalie - passage q2m/t2m pour calcul Rveget
-     & q2m, t2m)  
+     & q2m, t2m, &
+! Add emission/deposit fields
+     & field_out_names, fields_out, field_in_names, fields_in)  
 #else
   SUBROUTINE intersurf_gathered_2m (kjit, iim_glo, jjm_glo, kjpindex, kindex, xrdt, &
@@ -1866 +1874 @@
      & tsol_rad, temp_sol_new, qsurf, albedo, emis, z0, lon_scat_g, lat_scat_g, &
 ! Ajout Nathalie - passage q2m/t2m pour calcul Rveget
-     & q2m, t2m)  
+     & q2m, t2m, &
+! Add emission/deposit fields
+     & field_out_names, fields_out, field_in_names, fields_in)
 #endif
     ! routines called : sechiba_main
@@ -1925 +1935 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: fluxlat       !! Latent chaleur flux
     REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: emis          !! Emissivity
+    !
+    ! Optional arguments
+    !
+    ! Names and fields for emission variables : to be transport by GCM to chemistry model.
+    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_out_names
+    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: fields_out
+    !
+    ! Names and fields for deposit variables : to be transport from chemistry model by GCM to ORCHIDEE.
+    CHARACTER(LEN=*),DIMENSION(:), OPTIONAL, INTENT(IN) :: field_in_names
+    REAL(r_std),DIMENSION(:,:), OPTIONAL, INTENT(IN) :: fields_in
+    !
     ! LOCAL declaration
     ! work arrays to scatter and/or gather information just before/after sechiba_main call's
@@ -1937 +1958 @@
     REAL(r_std),DIMENSION (kjpindex)                      :: dcoastal      !! Work array to keep coastal flow
     REAL(r_std),DIMENSION (kjpindex)                      :: driver        !! Work array to keep river out flow
+    REAL(r_std),DIMENSION (kjpindex)                      :: znetco2       !! Work array to keep netco2flux
+    REAL(r_std),DIMENSION (kjpindex)                      :: zcarblu       !! Work array to keep fco2_land_use
     REAL(r_std),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
     REAL(r_std),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
@@ -1948 +1971 @@
     ! Optional arguments
     !
-    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN), OPTIONAL :: lon_scat_g, lat_scat_g !! The scattered values for longitude 
+    REAL(r_std),DIMENSION (iim_glo,jjm_glo), INTENT(IN) :: lon_scat_g, lat_scat_g !! The scattered values for longitude 
     !
     INTEGER(i_std)                          :: iim,jjm                                  !! local sizes
@@ -1980 +2003 @@
     LOGICAL, SAVE                                         :: fatmco2       !! Flag to force the value of atmospheric CO2 for vegetation.
     REAL(r_std), SAVE                                     :: atmco2        !! atmospheric CO2 
+    !
+    ! Number of fields to give (nb_fields_out) or get from (nb_fields_in) GCM :
+    INTEGER(i_std), SAVE                                  :: nb_fields_out, nb_fields_in
+    ! Id of fields to give (nb_fields_out) or get from (nb_fields_in) GCM :
+    INTEGER(i_std)                                        :: i_fields_out, i_fields_in
     !
     CALL ipslnlf(old_number=old_fileout)
@@ -2063 +2091 @@
        !  we have to do the work here.
        !
-       IF ( PRESENT(lon_scat_g) .AND. PRESENT(lat_scat_g)) THEN
+       IF ( .TRUE. ) THEN
           
           lon_scat(:,:)=zero
@@ -2081 +2109 @@
              lat_g(:,:) = lat_scat_g(:,:)
           ENDIF
-
-       ELSE IF ( PRESENT(lon_scat_g) .OR. PRESENT(lat_scat_g)) THEN
-
-          WRITE(numout,*) 'You need to provide the longitude AND latitude on the'
-          WRITE(numout,*) 'gathered grid in order to start ORCHIDEE.'
-          STOP 'intersurf_gathered'
 
        ELSE
@@ -2198 +2220 @@
        ENDIF
        !
+
+       ! Prepare fieds out/in for interface with GCM.
+       IF (PRESENT(field_out_names)) THEN
+          nb_fields_out=SIZE(field_out_names)
+       ELSE
+          nb_fields_out=0
+       ENDIF
+       IF (PRESENT(field_in_names)) THEN
+          nb_fields_in=SIZE(field_in_names)
+       ELSE
+          nb_fields_in=0
+       ENDIF
+
        IF ( check ) WRITE(numout,*) 'End of Initialisation of intersurf'
        !
@@ -2251 +2286 @@
        WRITE(numout,*) "Fraction of continent in the grid = ",zcontfrac
     ENDIF
+
+
+    ! Fields for deposit variables : to be transport from chemistry model by GCM to ORCHIDEE.
+    WRITE(numout,*) "Get fields from atmosphere."
+
+    DO i_fields_in=1,nb_fields_in
+       WRITE(numout,*) i_fields_in," Champ = ",TRIM(field_in_names(i_fields_in)) 
+       SELECT CASE(TRIM(field_in_names(i_fields_in)))
+       CASE DEFAULT 
+          CALL ipslerr (3,'intsurf_gathered_2m', &
+            &          'You ask in GCM an unknown field '//TRIM(field_in_names(i_fields_in))//&
+            &          ' to give to ORCHIDEE for this specific version.',&
+            &          'This model won''t be able to continue.', &
+            &          '(check your tracer parameters in GCM)')
+       END SELECT
+    ENDDO
+
     !
     ! 2. modification of co2
@@ -2301 +2353 @@
        & zprecip_rain ,zprecip_snow,  lwdown, swnet, swdown, pb, &
 ! Output : Fluxes
-       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, &
+       & zvevapp, zfluxsens, zfluxlat, zcoastal, zriver, znetco2, zcarblu, &
 ! Surface temperatures and surface properties
        & ztsol_rad, ztemp_sol_new, zqsurf, zalbedo, zemis, zz0, &
@@ -2553 +2605 @@
        riverflow(ik) = riverflow(ik)/xrdt
 
+    ENDDO
+    !
+    WRITE(numout,*) "Give fields to atmosphere."
+    
+    ! Fields for emission variables : to be transport by GCM to chemistry model.
+    DO i_fields_out=1,nb_fields_out
+       SELECT CASE(TRIM(field_out_names(i_fields_out)))
+       CASE("fCO2_land") 
+          fields_out(:,i_fields_out)=znetco2(:)
+       CASE("fCO2_land_use")
+          fields_out(:,i_fields_out)=zcarblu(:)
+       CASE DEFAULT 
+          CALL ipslerr (3,'intsurf_gathered_2m', &
+            &          'You ask from GCM an unknown field '//TRIM(field_out_names(i_fields_out))//&
+            &          ' to ORCHIDEE for this specific version.',&
+            &          'This model won''t be able to continue.', &
+            &          '(check your tracer parameters in GCM)')
+       END SELECT
     ENDDO
     !

branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/routing.f90

@@ -19 +19 @@
 !!
 !! @author Jan Polcher
-!! @Version : $Revision: 1.41 $, $Date: 2009/01/07 13:39:45 $
+!! @Version : $Revision: 274 $, $Date: 2011-06-21 15:18:18 +0200 (Tue, 21 Jun 2011) $
 !!
-!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/routing.f90,v 1.41 2009/01/07 13:39:45 ssipsl Exp $
+!< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/routing.f90 $
+!< $Date: 2011-06-21 15:18:18 +0200 (Tue, 21 Jun 2011) $
+!< $Author: martial.mancip $
+!< $Revision: 274 $
 !! IPSL (2006)
 !!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
@@ -285 +288 @@
        DO ig=1,nbpt
           IF ( lalo(ig,1) > 49.0 ) THEN
-             floodtemp(ig) = tp_00 -1.
+             floodtemp(ig) = tp_00 - un
           ENDIF
        ENDDO
@@ -499 +502 @@
     !Config If   = RIVER_ROUTING
     !Config Desc = Time step of th routing scheme
-    !Config Def  = 86400
+    !Config Def  = one_day
     !Config Help = This values gives the time step in seconds of the routing scheme. 
     !Config        It should be multiple of the main time step of ORCHIDEE. One day
@@ -569 +572 @@
        CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day)
        time_counter = tmp_day(1) 
-       CALL setvar (time_counter, val_exp, 'NO_KEYWORD', 0.0_r_std)
+       CALL setvar (time_counter, val_exp, 'NO_KEYWORD', zero)
     ENDIF
     CALL bcast(time_counter)
-!!$    CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', 0.0_r_std)
+!!$    CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', zero)
 
     !
@@ -650 +653 @@
     CALL ioconf_setatt('LONG_NAME','Water in the fast reservoir')
     CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g)
-    CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE (slow_reservoir(nbpt,nbasmax))
@@ -657 +660 @@
     CALL ioconf_setatt('LONG_NAME','Water in the slow reservoir')
     CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g)
-    CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE (stream_reservoir(nbpt,nbasmax))
@@ -664 +667 @@
     CALL ioconf_setatt('LONG_NAME','Water in the stream reservoir')
     CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g)
-    CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE (lake_reservoir(nbpt))
@@ -671 +674 @@
     CALL ioconf_setatt('LONG_NAME','Water in the lake reservoir')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g)
-    CALL setvar (lake_reservoir, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero)
     !
     ! Map of irrigated areas
@@ -689 +692 @@
     CALL ioconf_setatt('LONG_NAME','Previous outflow from this basin')
     CALL restget_p (rest_id, var_name, nbp_glo, nbasmax+3, 1, kjit, .TRUE., previous_outflow, "gather", nbp_glo, index_g)
-    CALL setvar_p (previous_outflow, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (previous_outflow, val_exp, 'NO_KEYWORD', zero)
     !
     IF ( dofloodplains ) THEN
@@ -707 +710 @@
     CALL ioconf_setatt('LONG_NAME','Lake inflow')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE (returnflow_mean(nbpt))
@@ -714 +717 @@
     CALL ioconf_setatt('LONG_NAME','Deep return flux')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero)
     returnflow(:) = returnflow_mean(:)
     !
@@ -726 +729 @@
        CALL ioconf_setatt('LONG_NAME','Artificial irrigation flux')
        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
-       CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
-       irrigation(:) = irrigation_mean(:) 
-   ELSE
+       CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
+    ELSE
        irrigation_mean(:) = zero
     ENDIF
+    irrigation(:) = irrigation_mean(:) 
     !
     ALLOCATE (riverflow_mean(nbpt))
@@ -737 +740 @@
     CALL ioconf_setatt('LONG_NAME','River flux into the sea')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
     riverflow(:) = riverflow_mean(:)
     !
@@ -745 +748 @@
     CALL ioconf_setatt('LONG_NAME','Diffuse flux into the sea')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
     coastalflow(:) = coastalflow_mean(:)
     !
@@ -759 +762 @@
     CALL ioconf_setatt('LONG_NAME','Hydrograph at outlow of grid')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g)
-    CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero)
     !
     ! The diagnostic variables, they are initialized from the above restart variables.
@@ -799 +802 @@
     CALL ioconf_setatt('LONG_NAME','Accumulated runoff for routing')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE(drainage_mean(nbpt))
@@ -806 +809 @@
     CALL ioconf_setatt('LONG_NAME','Accumulated drainage for routing')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE(evapot_mean(nbpt))
@@ -813 +816 @@
     CALL ioconf_setatt('LONG_NAME','Accumulated potential evaporation for routing')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (evapot_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (evapot_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE(precip_mean(nbpt))
@@ -820 +823 @@
     CALL ioconf_setatt('LONG_NAME','Accumulated rain precipitation for irrigation')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE(humrel_mean(nbpt))
@@ -827 +830 @@
     CALL ioconf_setatt('LONG_NAME','Mean humrel for irrigation')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', 1.0_r_std)
+    CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un)
     !
     ALLOCATE(totnobio_mean(nbpt))
@@ -834 +837 @@
     CALL ioconf_setatt('LONG_NAME','Last Total fraction of no bio for irrigation')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', 0.0_r_std)
+    CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero)
     !
     ALLOCATE(vegtot_mean(nbpt))
@@ -841 +844 @@
     CALL ioconf_setatt('LONG_NAME','Last Total fraction of vegetation')
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', 1.0_r_std)
+    CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
     !
     !
@@ -978 +981 @@
              !
           ELSE
-             fast_flow(ig,ib) = 0.0
-             slow_flow(ig,ib) = 0.0
-             stream_flow(ig,ib) = 0.0
+             fast_flow(ig,ib) = zero
+             slow_flow(ig,ib) = zero
+             stream_flow(ig,ib) = zero
           ENDIF
           inflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
@@ -994 +997 @@
              potflood(ig,ib) = inflow(ig,ib) - previous_outflow(ig,ib)
              !
-             IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN
+             IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > zero .AND. floodtemp(ig) > tp_00 ) THEN
                 !
                 IF (routing_area(ig,ib) > tobeflooded(ig)) THEN
                    floodindex = tobeflooded(ig) / routing_area(ig,ib)
                    ELSE
-                      floodindex = 1.0
+                      floodindex = un
                 ENDIF
                 !
@@ -1033 +1036 @@
 !ym mais n'est pas la plus efficace
 
-    IF (is_root_prc) &
-         ALLOCATE( fast_flow_g(nbp_glo, nbasmax), slow_flow_g(nbp_glo, nbasmax), &
-          stream_flow_g(nbp_glo, nbasmax), floods_g(nbp_glo, nbasmax), wdelay_g(nbp_glo, nbasmax) )
+    IF (is_root_prc)  THEN
+       ALLOCATE( fast_flow_g(nbp_glo, nbasmax), slow_flow_g(nbp_glo, nbasmax), &
+            stream_flow_g(nbp_glo, nbasmax), floods_g(nbp_glo, nbasmax),  &
+            wdelay_g(nbp_glo, nbasmax) )
+    ELSE
+       ALLOCATE( fast_flow_g(1,1), slow_flow_g(1,1), &
+            stream_flow_g(1, 1), floods_g(1,1),  &
+            wdelay_g(1,1) )
+    ENDIF
     
        
@@ -1057 +1066 @@
     ENDIF
 
-    IF (is_root_prc) &
-         DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g, floods_g, wdelay_g )
-   
+    DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g, floods_g, wdelay_g )
+      
     CALL scatter(transport_glo,transport)
 
@@ -1115 +1123 @@
        DO ig=1,nbpt
           
-          IF ((vegtot(ig) .GT. 0.0) .AND. (humrel(ig) .LT. 0.99)) THEN
-             irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(0.0, &
+          IF ((vegtot(ig) .GT. zero) .AND. (humrel(ig) .LT. 0.99)) THEN
+             irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, &
                   & crop_coef * evapot(ig) - &
                   & MAX(precip(ig)+returnflow(ig)-runoff(ig)-drainage(ig), zero) )
              irrig_netereq(ig) = 1 * irrig_netereq(ig)
              
-             IF(irrig_netereq(ig).LT.0.0) THEN
+             IF(irrig_netereq(ig).LT.zero) THEN
                 WRITE(numout,*) 'there is a probleme for irrig_netereq',ig,irrig_netereq(ig)
              ENDIF
@@ -1133 +1141 @@
                      &   stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) )
                 
-                slow_reservoir(ig,ib) = MAX(0.0, slow_reservoir(ig,ib) + &
-                     & MIN(0.0, fast_reservoir(ig,ib) + MIN(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib))))
+                slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + &
+                     & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib))))
 
-                fast_reservoir(ig,ib) = MAX( 0.0, &
-                     &  fast_reservoir(ig,ib) + MIN(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib)))
-                stream_reservoir(ig,ib) = MAX(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib) )
+                fast_reservoir(ig,ib) = MAX( zero, &
+                     &  fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib)))
+                stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib) )
 
-                IF(baseirrig(ig,ib) .LT. 0.0 .OR. slow_reservoir(ig,ib) .LT. 0.0 .OR. &
-                     & fast_reservoir(ig,ib) .LT. 0.0 .OR. stream_reservoir(ig,ib) .LT. 0.0) THEN
+                IF(baseirrig(ig,ib) .LT. zero .OR. slow_reservoir(ig,ib) .LT. zero .OR. &
+                     & fast_reservoir(ig,ib) .LT. zero .OR. stream_reservoir(ig,ib) .LT. zero) THEN
                    WRITE(numout,*) 'There is negative values related to irrigation', ig,ib,baseirrig(ig,ib), &
                         & slow_reservoir(ig,ib),fast_reservoir(ig,ib),stream_reservoir(ig,ib)
@@ -1378 +1386 @@
     !
     nb_pts(:) = 0
-    totarea(:) = 0.0
+    totarea(:) = zero
     hydrodiag(:,:) = 0
     DO ig=1,nbpt
@@ -1427 +1435 @@
     !
     !
-    basinmap(:) = 0.0
+    basinmap(:) = zero
     DO icc = 1, num_largest
        ff = MAXLOC(totarea)   
@@ -1445 +1453 @@
                & topids(ff(1)), name_str(1:15), totarea(ff(1))/1.e6,  nb_pts(ff(1))
        ENDIF
-       totarea(ff(1)) = 0.0
+       totarea(ff(1)) = zero
     ENDDO
     !
@@ -1471 +1479 @@
     WRITE(numout,*) 'Maximum topographic index :', MAXVAL(topo_resid)
     ic = COUNT(topo_resid .GT. 0.)
-    WRITE(numout,*) 'Mean topographic index :', SUM(topo_resid, MASK=topo_resid .GT. 0.)/ic
-    WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. 0.)
+    WRITE(numout,*) 'Mean topographic index :', SUM(topo_resid, MASK=topo_resid .GT. zero)/ic
+    WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. zero)
     !
     DEALLOCATE(pts)
@@ -1561 +1569 @@
     !  0.3 LOCAL
     !
-    !
     CHARACTER(LEN=80) :: filename
     INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, lastjp, nbexp
@@ -1674 +1681 @@
     nbexp = 0
     !
-    min_topoind = MINVAL(topoindex, MASK=topoindex .LT. undef_sechiba-1.)
+    min_topoind = MINVAL(topoindex, MASK=topoindex .LT. undef_sechiba-un)
     !
     DO ip=1,iml
@@ -1890 +1897 @@
        !
        !
-       ! Set everything to undef to locate easily empty points
-       !
-       trip_bx(:,:) = undef_int
-       basin_bx(:,:) = undef_int
-       topoind_bx(:,:) = undef_sechiba
-       area_bx(:,:) = undef_sechiba
-       hierarchy_bx(:,:) = undef_sechiba
-       !
        !  extract the information for this grid box
        !
@@ -2074 +2073 @@
     INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc
     REAL(r_std)    :: lonstr(nbvmax*nbvmax), latstr(nbvmax*nbvmax)
+    !
+    !
+    ! Set everything to undef to locate easily empty points
+    !
+    trip_bx(:,:) = undef_int
+    basin_bx(:,:) = undef_int
+    topoind_bx(:,:) = undef_sechiba
+    area_bx(:,:) = undef_sechiba
+    hierarchy_bx(:,:) = undef_sechiba
     !
     IF ( sub_pts(ib) > 0 ) THEN
@@ -3177 +3185 @@
        ! Compute the area of the basin
        !
-       basin_area(ib,ij) = 0.0
-       basin_hierarchy(ib,ij) = 0.0
+       basin_area(ib,ij) = zero
+       basin_hierarchy(ib,ij) = zero
        !
        SELECT CASE (hierar_method)
@@ -3186 +3194 @@
           !
        END SELECT
-       basin_topoind(ib,ij) = 0.0
+       basin_topoind(ib,ij) = zero
        !
        DO iz=1,basin_sz(ij)
@@ -3299 +3307 @@
     INTEGER(i_std) :: ff(1)
     !
+    ! WARNING
+    LOGICAL, PARAMETER :: check = .FALSE.
+    ! ERRORS 
+    LOGICAL :: error1, error2, error3, error4, error5
+    
+    error1=.FALSE.
+    error2=.FALSE.
+    error3=.FALSE.
+    error4=.FALSE.
+    error5=.FALSE.
+
     outflow_basin(:,:) = undef_int
     inflow_number(:,:) = 0
@@ -3373 +3392 @@
                       inflow_basin(inp, bop, inflow_number(inp,bop)) = sb
                    ELSE
-                      WRITE(numout,*) 'Increase nbvmax'
-                      STOP 'routing_linkup'
+                      error1=.TRUE.
+                      EXIT
                    ENDIF
                 ENDIF
@@ -3528 +3547 @@
                                   dop = sp
                                   bop = sbl
-                                  IF (basin_hierarchy(sp,sb) .LT. basin_hierarchy(sp,sbl)) THEN
-                                     WRITE(numout,*) '>>>>>>> POINT CORRECTED against hierarchy :',&
-                                          & sp, sb, 'into', sbl
+                                  IF (check) THEN
+                                     IF (basin_hierarchy(sp,sb) .LT. basin_hierarchy(sp,sbl)) THEN
+                                        WRITE(numout,*) '>>>>>>> POINT CORRECTED against hierarchy :',&
+                                             & sp, sb, 'into', sbl
+                                     ENDIF
                                   ENDIF
                                ENDIF
@@ -3542 +3563 @@
                 !
                 IF ( dop .EQ. undef_int .AND. bop .EQ. undef_int ) THEN
-                   WRITE(numout,*) 'Why are we here with point ', sp, sb
-                   WRITE(numout,*) 'Coodinates : (lon,lat) = ', lalo(sp,2), lalo(sp,1)
-                   WRITE(numout,*) 'Contfrac : = ', contfrac(sp)
-                   WRITE(numout,*) 'Local Basin ID :', basin_id(sp,1:basin_count(sp))
-                   WRITE(numout,*) 'Local hierarchies :', basin_hierarchy(sp,1:basin_count(sp))
-                   WRITE(numout,*) 'Local basin_flowdir :', basin_flowdir(sp,1:basin_count(sp))
-                   WRITE(numout,*) 'Local outflowgrid :', outflow_grid(sp,1:basin_count(sp))
-                   WRITE(numout,*) 'outflow_grid :', inp
-                   WRITE(numout,*) 'Coodinates outflow : (lon,lat) = ', lalo(inp,2), lalo(inp,1)
-                   WRITE(numout,*) 'Contfrac : = ', contfrac(inp)
-                   WRITE(numout,*) 'Outflow Basin ID :', basin_id(inp,1:basin_count(inp))
-                   WRITE(numout,*) 'Outflow hierarchies :', basin_hierarchy(inp,1:basin_count(inp))
-                   WRITE(numout,*) 'Outflow basin_flowdir :', basin_flowdir(inp,1:basin_count(inp))
-                   WRITE(numout,*) 'Explored options +1 :', dp1, bp1, outdp1
-                   WRITE(numout,*) 'Explored +1 Basin ID :', basin_id(dp1,1:basin_count(dp1))
-                   WRITE(numout,*) 'Explored +1 hierarchies :', basin_hierarchy(dp1,1:basin_count(dp1))
-                   WRITE(numout,*) 'Explored +1 basin_flowdir :', basin_flowdir(dp1,1:basin_count(dp1))
-                   WRITE(numout,*) 'Explored options -1 :', dm1, bm1, outdm1
-                   WRITE(numout,*) 'Explored -1 Basin ID :', basin_id(dm1,1:basin_count(dm1))
-                   WRITE(numout,*) 'Explored -1 hierarchies :', basin_hierarchy(dm1,1:basin_count(dm1))
-                   WRITE(numout,*) 'Explored -1 basin_flowdir :', basin_flowdir(dm1,1:basin_count(dm1))
-                   WRITE(numout,*) '****************************'
+                   IF (check) THEN
+                      WRITE(numout,*) 'Why are we here with point ', sp, sb
+                      WRITE(numout,*) 'Coodinates : (lon,lat) = ', lalo(sp,2), lalo(sp,1)
+                      WRITE(numout,*) 'Contfrac : = ', contfrac(sp)
+                      WRITE(numout,*) 'Local Basin ID :', basin_id(sp,1:basin_count(sp))
+                      WRITE(numout,*) 'Local hierarchies :', basin_hierarchy(sp,1:basin_count(sp))
+                      WRITE(numout,*) 'Local basin_flowdir :', basin_flowdir(sp,1:basin_count(sp))
+                      WRITE(numout,*) 'Local outflowgrid :', outflow_grid(sp,1:basin_count(sp))
+                      WRITE(numout,*) 'outflow_grid :', inp
+                      WRITE(numout,*) 'Coodinates outflow : (lon,lat) = ', lalo(inp,2), lalo(inp,1)
+                      WRITE(numout,*) 'Contfrac : = ', contfrac(inp)
+                      WRITE(numout,*) 'Outflow Basin ID :', basin_id(inp,1:basin_count(inp))
+                      WRITE(numout,*) 'Outflow hierarchies :', basin_hierarchy(inp,1:basin_count(inp))
+                      WRITE(numout,*) 'Outflow basin_flowdir :', basin_flowdir(inp,1:basin_count(inp))
+                      WRITE(numout,*) 'Explored options +1 :', dp1, bp1, outdp1
+                      WRITE(numout,*) 'Explored +1 Basin ID :', basin_id(dp1,1:basin_count(dp1))
+                      WRITE(numout,*) 'Explored +1 hierarchies :', basin_hierarchy(dp1,1:basin_count(dp1))
+                      WRITE(numout,*) 'Explored +1 basin_flowdir :', basin_flowdir(dp1,1:basin_count(dp1))
+                      WRITE(numout,*) 'Explored options -1 :', dm1, bm1, outdm1
+                      WRITE(numout,*) 'Explored -1 Basin ID :', basin_id(dm1,1:basin_count(dm1))
+                      WRITE(numout,*) 'Explored -1 hierarchies :', basin_hierarchy(dm1,1:basin_count(dm1))
+                      WRITE(numout,*) 'Explored -1 basin_flowdir :', basin_flowdir(dm1,1:basin_count(dm1))
+                      WRITE(numout,*) '****************************'
+                   ENDIF
                    IF ( contfrac(sp) > 0.01 ) THEN
-                      CALL ipslerr(3,'routing_linkup', &
-                        &      'In the routine which make connections between the basins and ensure global coherence,', & 
-                        &      'there is a problem with outflow linkup without any valid direction.', &
-                        &      '(Perhaps there is a problem with the grid.)')
+                      error2=.TRUE.
+                      EXIT
                    ENDIF
                 ENDIF
@@ -3585 +3606 @@
                    inflow_basin(dop, bop, inflow_number(dop,bop)) = sb
                 ELSE
-                   WRITE(numout,*) 'Increase nbvmax'
-                   STOP 'routing_linkup'
+                   error3=.TRUE.
+                   EXIT
                 ENDIF
                 !
@@ -3605 +3626 @@
              !
              
-             WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb
-             
+             IF (check) &
+                  WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb
              !
              DO sbl=1,basin_count(sp)
@@ -3622 +3643 @@
                    IF ( inflow_number(sp,sbl) .LE. nbvmax ) THEN
                       IF ( sp .EQ. 42 .AND. sbl .EQ. 1) THEN
-                         WRITE(numout,*) 'ADD INFLOW (3):', sp, sb
+                         IF (check) &
+                              WRITE(numout,*) 'ADD INFLOW (3):', sp, sb
                       ENDIF
                       inflow_grid(sp, sbl, inflow_number(sp,sbl)) = sp
                       inflow_basin(sp, sbl, inflow_number(sp,sbl)) = sb
                    ELSE
-                      WRITE(numout,*) 'Increase nbvmax'
-                      STOP 'routing_linkup'
+                      error4=.TRUE.
+                      EXIT
                    ENDIF
                 ENDIF
@@ -3639 +3661 @@
                & .AND. basin_flowdir(sp,sb) .GT. 0) THEN
              !
-             WRITE(numout,*) 'We could not find the basin into which we need to flow'
-             WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb
-             WRITE(numout,*) 'Explored neighbours :', dm1, dp1 
-             WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb)
-             WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb)
-             WRITE(numout,*) 'basin ID:',basin_id(sp,sb)
-             WRITE(numout,*) 'Hierarchy :', basin_hierarchy(sp,sb)
-             STOP 'routing_linkup'
+             error5=.TRUE.
+             EXIT
           ENDIF
        ENDDO
        !
     ENDDO
+    IF (error1) THEN
+       WRITE(numout,*) " routing_linkup : bop .LT. undef_int",bop
+       CALL ipslerr(3,'routing_linkup', &
+            "bop .LT. undef_int",'Increase nbvmax','stop routing_linkup')
+    ENDIF
+    IF (error2) THEN
+       CALL ipslerr(3,'routing_linkup', &
+            &      'In the routine which make connections between the basins and ensure global coherence,', & 
+            &      'there is a problem with outflow linkup without any valid direction. Try with check=.TRUE.', &
+            &      '(Perhaps there is a problem with the grid.)')
+    ENDIF
+    IF (error3) THEN
+       WRITE(numout,*) " routing_linkup : dop .GT. 0 .AND. dop .NE. undef_int",dop
+       CALL ipslerr(3,'routing_linkup', &
+            "dop .GT. 0 .AND. dop .NE. undef_int",'Increase nbvmax. Try with check=.TRUE.','stop routing_linkup')
+    ENDIF
+    IF (error4) THEN
+       WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", & 
+            &  " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, & 
+            &  basin_hierarchy(sp,sb),basin_hierarchy(sp,sbl)
+       CALL ipslerr(3,'routing_linkup', &
+            "(sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))" &
+            ,'Increase nbvmax. Try with check=.TRUE.','stop routing_linkup')
+    ENDIF
+    IF (error5) THEN
+       WRITE(numout,*) 'We could not find the basin into which we need to flow'
+       WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb
+       WRITE(numout,*) 'Explored neighbours :', dm1, dp1 
+       WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb)
+       WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb)
+       WRITE(numout,*) 'Outlfow basin :',outflow_basin(sp,sb)
+       WRITE(numout,*) 'basin ID:',basin_id(sp,sb)
+       WRITE(numout,*) 'Hierarchy :', basin_hierarchy(sp,sb)
+       CALL ipslerr(3,'routing_linkup', &
+            "We could not find the basin into which we need to flow",'Try with check=.TRUE.','stop routing_linkup')
+    ENDIF
     !
     ! Check for each outflow basin that it exists
@@ -3724 +3776 @@
     ! Compute the area upstream of each basin
     !
-    fetch_basin(:,:) = 0.0
+    fetch_basin(:,:) = zero
     !
     !
@@ -3793 +3845 @@
        ff = MAXLOC(tmp_area(1:nboutflow))
        outflow_grid(tmpindex(ff(1),1), tmpindex(ff(1),2)) = -1
-       tmp_area(ff(1)) = 0.0
+       tmp_area(ff(1)) = zero
     ENDDO
     !
@@ -3944 +3996 @@
              ! Now the take the smalest to be transfered to the largest
              !
-             iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. 0.)
+             iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
              sbas = multbas_sz(iml(1))
-             iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. 0.)
+             iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
              kbas = multbas_sz(iml(1))
              !
@@ -3991 +4043 @@
                    tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
                 ENDDO
-                iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.)
+                iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
                 sbas = multbas_list(ik,iml(1))
-                iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.)
+                iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
                 kbas = multbas_list(ik,iml(1))
                 !
@@ -4037 +4089 @@
                 ! If one of the basins goes to the ocean then it is going to have the priority
                 !
-                tmp_area(:) = 0.
+                tmp_area(:) = zero
                 IF ( COUNT(outflow_grid(ib,multbas_list(ik,1:multbas_sz(ik))) .LT. 0) .GT. 0) THEN
                    DO ii=1,multbas_sz(ik)
@@ -4047 +4099 @@
                    ENDDO
                    ! take the smalest of the subbasins
-                   iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.)
+                   iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
                    kbas = multbas_list(ik,iml(1))
                 ELSE
@@ -4056 +4108 @@
                       tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
                    ENDDO
-                   iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.)
+                   iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
                    sbas = multbas_list(ik,iml(1))
-                   iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.)
+                   iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
                    kbas = multbas_list(ik,iml(1))
                    !
@@ -4147 +4199 @@
        route_togrid(ib,:) = ib
        route_tobasin(ib,:) = 0
-       routing_area(ib,:) = 0.0
+       routing_area(ib,:) = zero
        !
     ENDDO
@@ -4210 +4262 @@
     ! Verify areas of the contienents
     !
-    floflo(:,:) = 0.0
+    floflo(:,:) = zero
     gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2)
     DO ib=1,nbpt
@@ -4249 +4301 @@
     !
     DO ib=1,nbpt
-       IF ( gridbasinarea(ib) > 0. ) THEN
+       IF ( gridbasinarea(ib) > zero ) THEN
           ratio = gridarea(ib)/gridbasinarea(ib)
           routing_area(ib,:) = routing_area(ib,:)*ratio
@@ -4268 +4320 @@
           largest_basins(ibf,:) = ff(:)
        ENDIF
-       floflo(ff(1), ff(2)) = 0.0
+       floflo(ff(1), ff(2)) = zero
     ENDDO
     !
@@ -4405 +4457 @@
     basin_flowdir(ib, tokill:basin_count(ib)-1) = basin_flowdir(ib, tokill+1:basin_count(ib))
     basin_area(ib, tokill:basin_count(ib)-1) = basin_area(ib, tokill+1:basin_count(ib))
-    basin_area(ib, basin_count(ib):nwbas) = 0.0
+    basin_area(ib, basin_count(ib):nwbas) = zero
     basin_topoind(ib, tokill:basin_count(ib)-1) = basin_topoind(ib, tokill+1:basin_count(ib))
-    basin_topoind(ib, basin_count(ib):nwbas) = 0.0
+    basin_topoind(ib, basin_count(ib):nwbas) = zero
     fetch_basin(ib, tokill:basin_count(ib)-1) = fetch_basin(ib, tokill+1:basin_count(ib))
-    fetch_basin(ib, basin_count(ib):nwbas) = 0.0
+    fetch_basin(ib, basin_count(ib):nwbas) = zero
     !
     ! Before we remove the information from the outflow fields we have to correct the corresponding inflow fields
@@ -4919 +4971 @@
        DO jp=1,jml
           !
-          IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN
+          IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN
              irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100.
-             IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = 0.0
+             IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero
           ENDIF
           !
-          IF ( flood_frac(ip,jp) .LT. undef_sechiba-1.) THEN
+          IF ( flood_frac(ip,jp) .LT. undef_sechiba-un) THEN
              flood_frac(ip,jp) = flood_frac(ip,jp)/100
-             IF ( flood_frac(ip,jp) < 0.005 )  flood_frac(ip,jp) = 0.0
+             IF ( flood_frac(ip,jp) < 0.005 )  flood_frac(ip,jp) = zero
           ENDIF
           !
@@ -5082 +5134 @@
                    ENDIF
                    !
-                   IF (flood_frac(ip,jp) .LT. undef_sechiba-1.) THEN
+                   IF (flood_frac(ip,jp) .LT. undef_sechiba-un) THEN
                       area_flood = area_flood + ax*ay*flood_frac(ip,jp)
                    ENDIF                      
@@ -5105 +5157 @@
           ENDIF
           ! Compute a diagnostic of the map.
-          IF(contfrac(ib).GT.0.0) THEN
+          IF(contfrac(ib).GT.zero) THEN
              irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) )
           ELSE

branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/sechiba.f90

@@ -4 +4 @@
 !!
 !! @author Marie-Alice Foujols and Jan Polcher
-!! @Version : $Revision: 45 $, $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
+!! @Version : $Revision: 275 $, $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
 !! 
 !< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/sechiba.f90 $
-!< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
-!< $Author: mmaipsl $
-!< $Revision: 45 $
+!< $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
+!< $Author: martial.mancip $
+!< $Revision: 275 $
 !! IPSL (2006)
 !!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
@@ -186 +186 @@
     & precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
          ! Output : Fluxes
-    & vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
+    & vevapp, fluxsens, fluxlat, coastalflow, riverflow, netco2flux, fco2_lu, &
          ! Surface temperatures and surface properties
     & tsol_rad, temp_sol_new, qsurf_out, albedo_out, emis_out, z0_out, &
@@ -249 +249 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat          !! Latent chaleur flux
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out         !! Emissivity
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux       !! Sum CO2 flux over PFTs (gC/m**2 of average ground/s)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu          !! Land Cover Change CO2 flux (gC/m**2 of average ground/s)
 
     REAL(r_std), ALLOCATABLE, DIMENSION (:)                  :: runoff1,drainage1, soilcap1,soilflx1
@@ -322 +324 @@
             lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
             rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-            co2_flux)
+            co2_flux, fco2_lu)
+       netco2flux(:) = zero
+       DO jv = 2,nvm
+          netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
+       ENDDO
        ! 
        ! computes initialisation of diffusion coeff
@@ -570 +576 @@
          lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
          rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-         co2_flux)
-
+         co2_flux, fco2_lu)
+    !
+    ! Compute global CO2 flux
+    !
+    netco2flux(:) = zero
+    DO jv = 2,nvm
+       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
+    ENDDO
     !
     ! call swap from new computed variables  
@@ -831 +843 @@
             lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
             rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-            co2_flux)
-
+            co2_flux, fco2_lu)
+       netco2flux(:) = zero
+       DO jv = 2,nvm
+          netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
+       ENDDO
 
        var_name= 'shumdiag'  

branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/slowproc.f90

@@ -3 +3 @@
 !
 !< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/slowproc.f90 $
-!< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
-!< $Author: mmaipsl $
-!< $Revision: 45 $
+!< $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
+!< $Author: martial.mancip $
+!< $Revision: 275 $
 !! IPSL (2006)
 !!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
@@ -78 +78 @@
        lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
        rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-       co2_flux)
+       co2_flux, fco2_lu)
 
 
@@ -117 +117 @@
     ! output fields
     REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out)      :: co2_flux         !! CO2 flux in gC/m**2 of average ground/second
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)          :: fco2_lu          !! Land Cover Change CO2 flux (gC/m**2 of average ground/s)
     ! modified scalar
     ! modified fields
@@ -190 +191 @@
                veget_nextyear, totfrac_nobio_nextyear, &
                hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-               co2_flux,resp_maint,resp_hetero,resp_growth)
+               co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
           !
        ENDIF
@@ -286 +287 @@
                veget_nextyear, totfrac_nobio_nextyear, &
                hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-               co2_flux,resp_maint,resp_hetero,resp_growth)
+               co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
        ENDIF
 
@@ -384 +385 @@
             veget_nextyear, totfrac_nobio_nextyear, &
             hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-            co2_flux,resp_maint,resp_hetero,resp_growth)
+            co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
        IF ( control%ok_stomate .AND. control%ok_sechiba ) THEN
           CALL histwrite(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
@@ -756 +757 @@
     CALL getin_p ("HYDROL_SOIL_DEPTH", dpu_cste)
     dpu(:)=dpu_cste
-    !
-    !Config Key  = HYDROL_HUMCSTE
-    !Config Desc = Root profile
-    !Config Def  = 5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4.
-    !Config Help = Default values were defined for 2 meters soil depth.
-    !Config        For 4 meters soil depth, you may use those ones :
-    !Config        5., .4, .4, 1., .8, .8, 1., 1., .8, 4., 1., 4., 1.
-    !
-!    humcste(:)= &
-!         & (/5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4./)
-!!$    CALL getin_p ("HYDROL_HUMCSTE", humcste)
-
 !MM, T. d'O. : before in constantes_soil :
 !          diaglev = &
@@ -2917 +2906 @@
     !
     IF (MAXVAL(vegmap) .LT. nolson) THEN
-      WRITE(*,*) 'WARNING -- WARNING'
-      WRITE(*,*) 'The vegetation map has to few vegetation types.'
-      WRITE(*,*) 'If you are lucky it will work but please check'
+       WRITE(*,*) 'WARNING -- WARNING'
+       WRITE(*,*) 'The vegetation map has to few vegetation types.'
+       WRITE(*,*) 'If you are lucky it will work but please check'
     ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
-      WRITE(*,*) 'More vegetation types in file than the code can'
-      WRITE(*,*) 'deal with.: ',  MAXVAL(vegmap),  nolson
-      STOP 'slowproc_interpol'
+       WRITE(*,*) 'More vegetation types in file than the code can'
+       WRITE(*,*) 'deal with.: ',  MAXVAL(vegmap),  nolson
+       STOP 'slowproc_interpol'
     ENDIF
     !

branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate.f90

@@ -328 +328 @@
        &  veget_max_new, totfrac_nobio_new, &
        &  hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
-       &  co2_flux,resp_maint,resp_hetero,resp_growth)
+       &  co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
     !---------------------------------------------------------------------
     !
@@ -418 +418 @@
     !NV champs 2D 
     REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out)      :: co2_flux
+    REAL(r_std),DIMENSION(kjpindex),INTENT(out)      :: fco2_lu
     ! autotrophic respiration in gC/m**2 of surface/dt
     REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out)  :: resp_maint
@@ -598 +599 @@
        !
        ! 1.3.1 read STOMATE's start file
+       !
+       co2_flux(:,:) = zero
+       fco2_lu(:) = zero
        !
        CALL readstart &
@@ -1337 +1341 @@
                &             fpc_max)
 
+          !
+          ! fco2_lu --> luccarb
+          fco2_lu(:) = convflux(:) &
+               &             + cflux_prod10(:)  &
+               &             + cflux_prod100(:) &
+               &             + harvest_above(:)
           !
           ! 6.4 output: transform from dimension nvm to dimension nvm