Merge-Hydro Branch

This version will merge the version of ORCHIDEE developed at LMD with the 11-layer hydrological scheme into the latest version of the trunk (1.9.5.2). Based on the differences between the LMD version and those from which it originates (either 1_9 or 1_9_4_1) (see the attachments), several notes have been written in order to prepare the merge. Here below are the links to these notes.

The hydrology group is : isabelle.gouttevin, jan.polcher, aducharne, aurelien.campoy, matthieu.guimberteau, gerhard.krinner, nathalie.de-noblet, catherine.ottle, pierre.brender, tao.wang, frederique.cheruy, marie-alice.foujols, patricia.cadule.
Everybody with a svn login can(since 12 april 2012) download this branch but only people in the group can commit changes.

How to download ORCHIDEE Hydrology branch for use in offline configuration

1. Get modipsl

   > svn co http://forge.ipsl.jussieu.fr/igcmg/svn/modipsl/trunk modipsl
   

2. Get the ORCHIDEE TOOLS-box. You will be prompted for username and password. Use personal svn login (firstname.lastname) or old sechiba cvs login.

   > cd modipsl/util
   > svn co svn://forge.ipsl.jussieu.fr/orchidee/trunk/TOOLS
   

3. Get the hydrology version. svn_login can be your personal svn login or sechiba

   > TOOLS/recup_my_ORCHIDEE svn_login branches/Hydrology
   

How to download ORCHIDEE Hydrology branch for use in coupled LMDZOR configuration

1. Get modipsl

   > svn co http://forge.ipsl.jussieu.fr/igcmg/svn/modipsl/trunk modipsl
   

2. Modifiy in modipsl/util/mod.def to get ORCHIDEE Hydrology sources instead of default.

   In modipsl/util/mod.def change following line :
   #-C- LMDZOR_v4  ORCHIDEE                    orchidee_1_9_5   2  .   modeles
   into
   #-C- LMDZOR_v4  branches/Hydrology/ORCHIDEE HEAD             15 .  modeles
   

3. Extract configuration LMDZOR_v4 and compile, for more information see ​wiki LMDZOR_v4

   cd modipsl/util
   ./model LMDZOR_v4       => you will be asked for passwd for sechiba
   ./ins_make
   cd ../config/LMDZOR_v4
   gmake
   

4. Some changes are needed to run this version using Choisnel 2 layers. No default set up is done yet for using CWRR.
   • Add in PARAM/orchidee.def :

     SOILCLASS_FILE = soils_param.nc
     TOPOGRAPHY_FILE =cartepente2d_15min.nc
     

   • Add in orchidee.card, in section [InitialStateFiles] following line (and add ", \" at line above):

     (${R_INIT}/SRF/${config_UserChoices_TagName}/cartepente2d_15min.nc, .)
     

   • NB! For the moment this version bug with default set up using RIVER_ROUTING=y. To have it runnig, change in orchidee.def to RIVER_ROUTING=n.

Notes to prepare the merge

Note on CWRR

Note on evapnu

Note on dpu

Note on flag couple

Note on logz0

Note on LAI

Note on Soil Map

Note on Interpolations

Martial notes for the merge

Testing the Merge-Hydro version

Matthieu Guimberteau: Off-line tests in the Amazon with bug corrections

Aurélien Campoy (Feb 2012) - test peformed off-line using 50 years of SAFRAN met forcing for the SIRTA site, without routing, without STOMATE : compares reasonably well with the CWRR version of ORCHIDEE before the merge

Jan Polcher (March 2012) - Off-line simulation over Europe using the WATCH forcing data set : simulation performed, to be analyzed

Nicolas Vuichard (March 2012) - off-line tests at FLUXNET sites using STOMATE : ongoing. See details

Fabienne Maignan (Feb-Jul 2012) - global off-line simulation

• forcing: ERA-Interim regular grid 0.72° 3-hourly
• computer: LSCE/obelix*
• 16 procs
• CWRR + STOMATE + ROUTING
  • segmentation fault SIGSEGV here:

     Projection arrays for Slope map                      : 
     nbvmax =           16
     We will work with      1000105  points of the fine grid
     Aggregate_2d : Slope map                     
     
     aggregate_2D nbvmax =           16 max used =           16
    

    The problem is in slowproc_slope after the call to aggregate_vec:

        DO ib = 1, nbpt
          idi=1
          !-
          !- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro 
          !-
          slopecoef = zero
    !@BUG SOURCE      DO WHILE ( sub_area(ib,idi) > zero ) !@END BUG SOURCE 
    !@PROPOSED CORRECTION
           DO WHILE ( sub_area(ib,idi) > zero .AND. idi .LE. nbvmax) 
    !@END PROPOSED CORRECTION
             ip = sub_index(ib,idi,1)
             jp = sub_index(ib,idi,2)
             !
             slopecoef = slopecoef + MIN(slopemap(ip,jp)/slope_noreinf, un) * sub_area(ib,idi)
             idi = idi +1
          ENDDO
    

    idi may become greater than nbvmax, which causes the segmentation fault. The problem is resolved by the correction proposed in the above code.
  • stop in hydrol_soil_infilt:

    Error in the calculation of infilt tot -4.176114610371373E-008
     k, ji, jst, mc   1.17784270216726        2.16374183793683             2394
               1  0.277241865204730
    
    FATAL ERROR FROM ROUTINE hydrol_soil_infilt
     --> We will STOP after hydrol_soil_infilt.
     -->
     -->
    
    Fatal error from IOIPSL. STOP in ipslerr with code
    

    This is caused by slightly negative values in the ERA-Interim precipitations, as already mentionned by Nicolas Vuichard. Proposed (and validated) modification in intersurf_main_2d and intersurf_main_1d:

    !@ PROPOSED MODIFICATION
        WHERE(zprecip_rain(:) .LT. 0.)
    	zprecip_rain(:)=0.		
        ENDWHERE				
    !@ END PROPOSED MODIFICATION
        !
        IF (check_INPUTS) THEN
    ...
    

• test option CHECK_CWRR: OK (no stop, much longer ~*3)
• Comparison Choisnel against 196: on-going

Frédérique Cheruy - global simulation coupled to LMDz : to be done

Merge into the trunk (rev 941) (Didier Solyga)

This section lists the modifications needed to merge the DOC/Hydro branch into the trunk :

• hydrol : Move flag ok_throughfall_by_pft to pft_parameters :

       IF ( active_flags%hydrol_cwrr ) THEN
           
           !! 2.1 Read the flag ok_throughfall_by_pft to know if 
           !!      we have to use the parameter throughfall_by_pft
  
           !Config Key   = OK_THROUGHFALL_PFT
           !Config Desc  = Activate use of PERCENT_THROUGHFALL_PFT
           !Config If    = HYDROL_CWRR
           !Config Def   = FALSE
           !Config Help  = If NOT OFF_LINE_MODE it is always TRUE (coupled with a GCM)
           !Config Units = [FLAG]
           IF ( .NOT. OFF_LINE_MODE ) ok_throughfall_by_pft = .TRUE.
           CALL getin_p('OK_THROUGHFALL_PFT',ok_throughfall_by_pft)   
  
        END IF
  
  

The pft parameter throughfall_by_pft is initiliazed and read in pft_parameters.f90. Correct memory allocation and initialization for parameter throughfall_by_pft :

      IF ( .NOT.(active_flags%hydrol_cwrr) .OR. (active_flags%hydrol_cwrr .AND. ok_throughfall_by_pft) ) THEN
         ALLOCATE(throughfall_by_pft(nvm),stat=ier)
         l_error = l_error .OR. (ier /= 0)
         IF (l_error) THEN
            WRITE(numout,*) ' Memory allocation error for throughfall_by_pft. We stop. We need nvm words = ',nvm
            STOP 'pft_parameters_alloc'
         END IF
      END IF

      IF ( .NOT.(active_flags%hydrol_cwrr) .OR.  (active_flags%hydrol_cwrr .AND. ok_throughfall_by_pft) ) THEN
         throughfall_by_pft(:) = throughfall_by_mtc(pft_to_mtc(:))
      ENDIF

because throughfall_by_pft is still used with Choisnel hydrology but not automatically with CWRR. ===> Q: Inconsistency ?

• pft_parameters : humcste has different default values according the value of dpu_max. We know that if we use the 11-layers hydrology, dpu_max should be set to 2 and humcste the corresponding values

     If (active_flags%hydrol_cwrr ) THEN
        humcste(:) = humcste_cwrr(pft_to_mtc(:)) ! values for 2m soil depth
     ELSE
        humcste(:) = humcste_mtc(pft_to_mtc(:))  ! values for 4m soil depth 
     END IF
  

• constantes_soil : reintegrate the module constantes_soil. Add the corresponding "USE" in the code. The old hydrological parameters externalized which are obsolete now have been commented in constantes.f90. dpu_max is set at 4 meters by default (value used for AR5 + Choisnel)
• intersurf : add a consistency test for dpu_max value. dpu_max can't be set to 4 if hydrol_cwrr is activated.

   IF (control_flags%hydrol_cwrr .AND. (dpu_max /= 2.)) THEN
       WRITE (numout,*) "You can not use the 11-layers hydrology with dpu_max /= 2. We set it to 2."
       dpu_max = 2.
    END IF

• routing : Add documentation. the following lines have been adapted for the externalization :

    DO ig = 1, nbpt
       IF (MAXVAL(veget_max(ig,(nvm-3):nvm)) .GT. min_sechiba) THEN
          DO jv = nvm-3, nvm
             transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,(nvm-3):nvm))
          ENDDO
       ELSE

The corresponding code now is :

   tot_vegfrac_nowoody(:) = zero
   DO jv  = 1, nvm
      IF (is_c3(jv) .OR. is_c4(jv)) THEN
         tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv) 
      END IF
   END DO

   DO ig = 1, nbpt
      IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN
         DO jv = 1,nvm
            IF ( is_c3(jv) .OR. is_c4(jv) ) THEN
               transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig)  
            END IF
         END DO

• MERGE DONE : Merge Hydrology(+Doc) branch revision 932 into trunk revision 945. Done on 19/07/2012 look at [947].

Differences explaining the differences between the trunk and the Hydrology branch :

• Initialization of lai : in the branch (and tag 196), the lai is initialized twice if we have no restart files, once in slowproc_lai, the second time in stomate.f90. In slowproc_lai, the lai is non-zero :

    IF  ( .NOT. read_lai ) THEN
    
       lai(: ,1) = zero
       ! On boucle sur 2,nvm au lieu de 1,nvm
       DO jv = 2,nvm
          SELECT CASE (type_of_lai(jv))
             
          CASE ("mean ")
             !
             ! 1. do the interpolation between laimax and laimin
             !
             lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
             !
          CASE ("inter")
             !
             ! 2. do the interpolation between laimax and laimin
             !
             DO ji = 1,kjpindex
                lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
             ENDDO
             !
          CASE default
             !
             ! 3. Problem
             !
             WRITE (numout,*) 'This kind of lai choice is not possible. '// &
                  ' We stop with type_of_lai ',jv,' = ', type_of_lai(jv) 
             STOP 'slowproc_lai'
             
          END SELECT

but once in stomate_f90, it is recalculated :

       IF (control%ok_pheno) THEN
          !! 5.1.1 Update LAI 
          ! Set lai of bare soil to zero
          lai(:,ibare_sechiba) = zero
          ! lai for all PFTs
          DO j = 2, nvm
             lai(:,j) = biomass(:,j,ileaf)*sla(j)
          ENDDO
       ELSE 
          ! 5.1.2 Use a prescribed lai
          ! WARNING: code in setlai is identical to the lines above
          ! Update subroutine if LAI has to be forced 
          CALL  setlai(kjpindex,lai) 
       ENDIF

lai could be positive in sechiba but once in stomate, it could be equal to zero! It was inconsistent with the BVOC : we have isopren emissions without vegetation! I add the following line (see [890]) in slowproc (after calling slowproc_lai):

          IF ( .NOT. read_lai ) THEN
             lai(:,:) = zero
          ENDIF

This modification was not possible in slowproc_lai because the error occurs only when we don't use restart files.

• The second is the modification of solar module by N.Vuichard (see [861]) for using daily forcings. lhour is considered as a real and not an integer anymore.
• The sla parameter is calculated in Hydrology branch but prescribed in the trunk. The values are relatively the same but could explain some little differences.
• The CMOR outputs mrros, mrr, prveg, evspsblveg, evspsblsoi, tran have been corrected in the trunk (see ticket #9). Set SECHIBA_HIST_LEVEL to 10 to not activate them.

Differences explaining the differences between the trunk and the tag 196 :

• In tag 196, when IMPOSE_VEG is activated veget is not calculated by the subroutine slowproc_veget. So as we have veget=veget_max, no fraction of bare soil is calculated.
• In trunk (rev [947]), veget(:,1) is replaced by tot_frac_bare and slowproc_veget is called whatever the case. So a fraction is always calculated even in IMPOSE_VEG. This modification implies a modification of the value of z0 (and albedo). To have the same results, a call to slowproc_veget should be added in the tag version of slowproc :

         !
         !Config Key   = SLOWPROC_HEIGHT
         !Config Desc  = Height for all vegetation types
         !Config Def   = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0
         !Config If    = OK_SECHIBA
         !Config Help  = The height used in the 0dim mode. The values should be found
         !Config         in the restart file. The new values of height will be computed anyway
         !Config         at the end of the current day. The need for this variable is caused
         !Config         by the fact that the model may stop during a day and thus we have not
         !Config         yet been through the routines which compute the new surface conditions.
         !Config Units = [m]
         !
         CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc)
  
         CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget)
  
      ELSE
  

This line lets ORCHIDEE to calculate a fraction of bare soil even in IMPOSE_VEG.