Changeset 7338 for branches/ORCHIDEE_2_2


Ignore:
Timestamp:
2021-11-06T12:00:01+01:00 (3 years ago)
Author:
agnes.ducharne
Message:

Improvement of r7337: (a) correct default value over Greenland with Zobler maps, involves the suppression of soilclass_default, replaced by scalar usda_default=6; (b) bug correction on max index for xios interpolation of Zobler map. Compiles and runs OK (though not tested with xios_interpolation).

Location:
branches/ORCHIDEE_2_2/ORCHIDEE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90

    r7337 r7338  
    113113 
    114114  INTEGER(i_std), PARAMETER,DIMENSION(nscm_fao) :: fao2usda = (/ 3,6,9 /) !! To find the values of Coarse, Medium, Fine in Zobler map 
    115   !! from the USDA lookup tables 
    116    
    117 !!$  REAL(r_std),DIMENSION(nscm_fao),SAVE :: soilclass_default_fao = &   !! Default soil texture distribution for fao : 
    118 !!$ & (/ 0.28, 0.52, 0.20 /)                                             !! in the following order : COARSE, MEDIUM, FINE (unitless) 
    119 !!$!$OMP THREADPRIVATE(soilclass_default_fao) 
    120  
     115                                                                          !! from the USDA lookup tables 
     116   
    121117  !!  2. Parameters for USDA Classification 
    122118 
    … …  
    127123  ! Greenland, where the Zobler map has no data, were changed, even when having 0.28 at the 3rd place below 
    128124  ! I kept the original soilclass_default_usda to keep the Reynolds map unchanged 
    129    
    130   REAL(r_std),DIMENSION(nscm_usda),SAVE :: soilclass_default_usda = &    !! Default soil texture distribution in the above order : 
    131  & (/ 0.28, 0.52, 0.20, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)   !! Thus different from "FAO"'s COARSE, MEDIUM, FINE 
    132   !! which have indices 3,6,9 in the 12-texture vector 
    133   !$OMP THREADPRIVATE(soilclass_default_usda) 
    134  
    135 !!$  REAL(r_std),DIMENSION(nscm_usda),SAVE :: soilclass_default_usda = &       !! Default soil texture distribution, to be coherent with the  
    136 !!$ & (/ 0.0, 0.0, 0.28, 0.0, 0.0, 0.52, 0.0, 0.0, 0.20, 0.0, 0.0, 0.0, 0.0 /) !! original FAO/Zobler values 
    137 !!$  !$OMP THREADPRIVATE(soilclass_default_usda)   
    138  
     125 
     126   INTEGER(i_std), SAVE      :: usda_default = 6            !! Default USDA texture class if no value found from map 
     127!$OMP THREADPRIVATE(usda_default) 
     128   
    139129  REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: nvan_usda = &            !! Van Genuchten coefficient n (unitless) 
    140130 & (/ 2.68_r_std, 2.28_r_std, 1.89_r_std, 1.41_r_std, &                   !  RK: 1/n=1-m 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

    r7337 r7338  
    7171!$OMP THREADPRIVATE(siltfraction)   
    7272  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: laimap              !! LAI map when the LAI is prescribed and not calculated by STOMATE 
    73 !$OMP THREADPRIVATE(laimap) 
    74   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: soilclass_default 
    75 !$OMP THREADPRIVATE(soilclass_default) 
    7673  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: veget_max_new       !! New year fraction of vegetation type (0-1, unitless) 
    7774!$OMP THREADPRIVATE(veget_max_new) 
    … …  
    948945    siltfraction(:)=undef_sechiba 
    949946 
    950     ! Initialisation of the fraction of the different vegetation: Start with 100% of bare soil 
    951     ALLOCATE (soilclass_default(nscm),stat=ier) 
    952     IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilclass_default','','') 
    953     soilclass_default(:)=undef_sechiba 
    954  
    955947    ! Allocation of last year vegetation fraction in case of land use change 
    956948    ALLOCATE(veget_max_new(kjpindex, nvm), STAT=ier) 
    … …  
    13001292       IF (impsoilt) THEN 
    13011293 
    1302           ! If njsc is not in restart file, then initialize soilclass from values  
    1303           ! from run.def file and recalculate njsc 
     1294          ! If njsc is not in restart file, then initialize to default value 
     1295          ! tbd - reading from run.def file  
    13041296          IF ( ALL(njsc(:) .EQ. undef_int )) THEN 
    1305              !Config Key   = SOIL_FRACTIONS 
    1306              !Config Desc  = Fraction of the 3 soil types (0-dim mode) 
    1307              !Config Def   = undef_sechiba 
    1308              !Config If    = IMPOSE_VEG and IMPOSE_SOILT 
    1309              !Config Help  = Determines the fraction for the 3 soil types 
    1310              !Config         in the mesh in the following order : sand loam and clay. 
    1311              !Config Units = [-] 
    1312            
    1313              soilclass(1,:) = soilclass_default(:) 
    1314              CALL getin_p('SOIL_FRACTIONS',soilclass(1,:)) 
    1315              ! Assign for each grid-cell the % of the different textural classes (up to 12 if 'usda') 
    1316              DO ji=2,kjpindex 
    1317                 ! here we read, for the prescribed grid-cell, the % occupied by each of the soil texture classes  
    1318                 soilclass(ji,:) = soilclass(1,:) 
    1319              ENDDO 
    1320  
    1321              ! Simplify an heterogeneous grid-cell into an homogeneous one with the dominant texture 
    1322              njsc(:) = 0 
    1323              DO ji = 1, kjpindex 
    1324                 ! here we reduce to the dominant texture class 
    1325                 njsc(ji) = MAXLOC(soilclass(ji,:),1) 
    1326              ENDDO 
     1297             njsc(:) = usda_default ! 6 = Loam 
    13271298          END IF 
    13281299 
    … …  
    13571328             CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, lalo, neighbours, resolution, contfrac, soilclass, & 
    13581329                  clayfraction, sandfraction, siltfraction) 
    1359              njsc(:) = 0 
    1360              DO ji = 1, kjpindex 
    1361                 njsc(ji) = MAXLOC(soilclass(ji,:),1) 
    1362              ENDDO 
    13631330             call_slowproc_soilt=.FALSE. 
    13641331          ENDIF 
    … …  
    14941461            clayfraction, sandfraction, siltfraction) 
    14951462       IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_soilt' 
    1496        njsc(:) = 0 
    1497        DO ji = 1, kjpindex 
    1498           njsc(ji) = MAXLOC(soilclass(ji,:),1) 
    1499        ENDDO 
    15001463       call_slowproc_soilt=.FALSE. 
    15011464    ENDIF 
    … …  
    16041567    IF (ALLOCATED (woodharvest)) DEALLOCATE (woodharvest) 
    16051568    IF (ALLOCATED (frac_nobio_new)) DEALLOCATE (frac_nobio_new) 
    1606     IF ( ALLOCATED (soilclass_default)) DEALLOCATE (soilclass_default) 
    16071569 
    16081570 ! 2. Clear all the variables in stomate  
    … …  
    29202882          sandfraction=0.81 
    29212883          siltfraction=0.13 
    2922          ! njsc=2 
    29232884          DO ib=1 , nbpt 
    29242885             njsc(ib) = 2 
    … …  
    29312892             mcw(ib) = 0.05221 
    29322893          ENDDO 
    2933           ! definir clayfraction à partir du barycentre du domaine loamy sand 
    2934           ! definir njsc pour cette classe => pcent dans hydrol 
    29352894 
    29362895       CASE ('b') !loam 
    … …  
    29382897          sandfraction=0.4 
    29392898          siltfraction=0.4 
    2940           !njsc=6 
    29412899          DO ib=1, nbpt 
    29422900             njsc(ib) = 6 
    … …  
    29492907             mcw(ib) = 0.09115 
    29502908          ENDDO 
    2951           ! definir clayfraction à partir du barycentre du domaine 
    2952           ! definir njsc pour cette classe 
    29532909 
    29542910       CASE ('c') !silt 
    … …  
    29562912          sandfraction=0.06 
    29572913          siltfraction=0.84 
    2958           !njsc=5 
    29592914          DO ib=1, nbpt 
    29602915             njsc(ib)=5 
    … …  
    29682923          ENDDO 
    29692924 
    2970           ! definir clayfraction à partir du barycentre du domaine 
    2971           ! definir njsc pour cette classe 
    29722925       CASE ('d')!clay 
    29732926          clayfraction=0.55 
    29742927          sandfraction=0.15 
    29752928          siltfraction=0.3 
    2976           !njsc=12 
    29772929          DO ib=1, nbpt 
    29782930             njsc(ib)=12 
    … …  
    29852937             mcw(ib) = 0.1897 
    29862938          ENDDO 
    2987           ! definir clayfraction à partir du barycentre du domaine 
    2988           ! definir njsc pour cette classe 
    29892939 
    29902940       CASE DEFAULT 
    … …  
    30322982             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    30332983             DO ib=1, nbpt 
    3034                 soilclass(ib,:) = soilclass_default_usda  
     2984                njsc(ib) = usda_default ! 6 = Loam 
    30352985                clayfraction(ib) = clayfraction_default 
    30362986             ENDDO 
    30372987 
    30382988          CASE('zobler') 
    3039              ! 
    3040              soilclass_default=soilclass_default_usda ! USDA means here 13 final texture classes, owing to fao2usda 
    3041              ! 
     2989             !             ! 
    30422990             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification, to be read using XIOS" 
    30432991             ! 
    … …  
    30763024                     textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7) 
    30773025 
    3078                 sgn=SUM(soilclass(ib,1:3)) ! grid-cell fraction with texture info 
    3079  
    3080                 IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume 28%/52%, 20% of texture classes 3/6/9 
    3081                    soilclass(ib,:) = soilclass_default(:) 
     3026                sgn=SUM(soilclass(ib,:)) ! grid-cell fraction with texture info 
     3027 
     3028                IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume that texture = Loam 
     3029                   njsc(ib) = usda_default ! 6 = Loam 
    30823030                   clayfraction(ib) = clayfraction_default 
    30833031                   sandfraction(ib) = sandfraction_default 
    … …  
    30893037                   sandfraction(ib) = sandfraction(ib) / sgn 
    30903038                   siltfraction(ib) = siltfraction(ib) / sgn 
    3091                    soilclass(ib,1:3) = soilclass(ib,1:3) / sgn 
    3092                 ENDIF 
     3039                   soilclass(ib,:) = soilclass(ib,:) / sgn 
     3040                   njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class 
     3041                ENDIF   
    30933042 
    30943043             ENDDO 
    … …  
    30973046 
    30983047             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda' 
    3099  
    3100              soilclass_default=soilclass_default_usda 
    31013048             ! 
    31023049             WRITE(numout,*) "Using a soilclass map with usda classification, to be read using XIOS" 
    … …  
    31343081                   sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf) 
    31353082                   siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf) 
     3083                   ! textfrac_table holds the %silt,%sand,%clay 
    31363084                ENDDO 
    31373085 
    31383086                sgn=SUM(soilclass(ib,:)) ! grid-cell fraction with texture info 
    31393087 
    3140                 IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume 28%/52%, 20% of texture classes 3/6/9 
    3141                    soilclass(ib,:) = soilclass_default(:) 
     3088                IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume that texture = Loam 
     3089                   njsc(ib) = usda_default ! 6 = Loam 
    31423090                   clayfraction(ib) = clayfraction_default 
    31433091                   sandfraction(ib) = sandfraction_default 
    … …  
    31503098                   siltfraction(ib) = siltfraction(ib) / sgn 
    31513099                   atext(ib)=sgn 
    3152                 ENDIF 
     3100                   njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class 
     3101                ENDIF                
     3102                     
    31533103             ENDDO 
    31543104 
    … …  
    32353185             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    32363186             DO ib=1, nbpt 
    3237                 soilclass(ib,:) = soilclass_default_usda 
     3187                njsc(ib) = usda_default ! 6 = Loam 
    32383188                clayfraction(ib) = clayfraction_default 
    32393189                sandfraction(ib) = sandfraction_default 
    … …  
    32413191             ENDDO 
    32423192          CASE('zobler') 
    3243              ! 
    3244              soilclass(ib,:) = soilclass_default_usda ! USDA means here 13 final texture classes, owing to fao2usda 
    3245              ! 
     3193             !             ! 
    32463194             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification" 
    32473195             ! 
    … …  
    32653213                   ! No points were found for current grid box, use default values 
    32663214                   nbexp = nbexp + 1 
    3267                    soilclass(ib,:) = soilclass_default(:) 
     3215                   njsc(ib) = usda_default ! 6=Loam 
    32683216                   clayfraction(ib) = clayfraction_default 
    32693217                   sandfraction(ib) = sandfraction_default 
    … …  
    33323280                      ! or if now information on the source grid was found. 
    33333281                      nbexp = nbexp + 1 
    3334                       soilclass(ib,:) = soilclass_default(:) 
     3282                      njsc(ib) = usda_default ! 6 = Loam 
    33353283                      clayfraction(ib) = clayfraction_default 
    33363284                      sandfraction(ib) = sandfraction_default 
    … …  
    33413289                      clayfraction(ib) = clayfraction(ib)/sgn 
    33423290                      sandfraction(ib) = sandfraction(ib)/sgn 
    3343                       siltfraction(ib) = siltfraction(ib)/sgn 
     3291                      siltfraction(ib) = siltfraction(ib)/sgn               
     3292                      njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class 
    33443293                   ENDIF 
    33453294                ENDIF 
    … …  
    33503299          CASE("usda") 
    33513300             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification" 
    3352  
    3353              soilclass_default=soilclass_default_usda 
    3354  
    33553301 
    33563302             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
    … …  
    33813327                   IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib 
    33823328                   nbexp = nbexp + 1 
    3383                    soilclass(ib,:) = soilclass_default 
     3329                   njsc(ib) = usda_default ! 6 = Loam 
    33843330                   clayfraction(ib) = clayfraction_default 
    33853331                   sandfraction(ib) = sandfraction_default 
    … …  
    34133359                      ! 
    34143360                   ENDDO 
     3361                   njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class 
    34153362 
    34163363                   ! Set default values if the surface in source file is too small 
     3364                   ! Warning - This test is donne differently for Zobler (based on sgn, related to class 6=ice) 
    34173365                   IF ( atext(ib) .LT. min_sechiba) THEN 
    34183366                      nbexp = nbexp + 1 
    3419                       soilclass(ib,:) = soilclass_default(:) 
     3367                      njsc(ib) = usda_default ! 6 = Loam 
    34203368                      clayfraction(ib) = clayfraction_default 
    34213369                      sandfraction(ib) = sandfraction_default 
    … …  
    34263374             ENDDO 
    34273375 
    3428              IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda' 
     3376          IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda' 
    34293377 
    34303378          CASE DEFAULT 
    … …  
    34483396 
    34493397       ENDIF        !      xios_interpolation 
    3450  
    3451        !salma: calculate njsc 
    3452        njsc(:) = 0 
    3453        DO ib = 1, nbpt 
    3454           njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Here we get 3/6/9 for the Zobler classes Coarse/Medium/Fine 
    3455        ENDDO 
    34563398 
    34573399       ! End of soil texture reading, for 'maps' and classical behavior 
Note: See TracChangeset for help on using the changeset viewer.