Changeset 7337 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba

Timestamp:

2021-11-06T08:24:04+01:00 (3 years ago)

Author:

agnes.ducharne

Message:

Simplification of soil texture processing, cf ticket #416: when using the Zobler map, the soil parameters are no more taken from 3-value "FAO" vectors in constantes_soil_var.f90, but from 13-value USDA vectors, owing to a pointer fao2usda. A 5-day running test with the Zobler map shows no changes, but over Greenland (class 6=ice in Zobler map).

Location:

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba

Files:

• hydrol.f90 (modified) (1 diff)
• slowproc.f90 (modified) (10 diffs)
• thermosoil.f90 (modified) (3 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90

@@ -1449 +1449 @@
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','')
        
-    !! 2.2 Soil texture choose
-
-    SELECTCASE (nscm)
-    CASE (3)              
-       pcent(:) = pcent_fao(:) 
-       mc_awet(:) = mc_awet_fao(:)
-       mc_adry(:) = mc_adry_fao(:)
-    CASE (13)           
-       pcent(:) = pcent_usda(:) 
-       mc_awet(:) = mc_awet_usda(:)
-       mc_adry(:) = mc_adry_usda(:)       
-    CASE DEFAULT
-       WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler and usda according to the map'
-       CALL ipslerr_p(3,'hydrol_init','Unsupported soil type classification. ',&
-            'Choose between zobler and usda according to the map','')
-    ENDSELECT
-
+    !! 2.2 Soil texture parameters
+         
+    pcent(:) = pcent_usda(:) 
+    mc_awet(:) = mc_awet_usda(:)
+    mc_adry(:) = mc_adry_usda(:) 
 
     !! 2.3 Read in the run.def the parameters values defined by the user

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

@@ -3032 +3032 @@
              IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
              DO ib=1, nbpt
-                soilclass(ib,:) = soilclass_default_fao
+                soilclass(ib,:) = soilclass_default_usda 
                 clayfraction(ib) = clayfraction_default
              ENDDO
@@ -3038 +3038 @@
           CASE('zobler')
              !
-             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
+             soilclass_default=soilclass_default_usda ! USDA means here 13 final texture classes, owing to fao2usda
              !
              IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification, to be read using XIOS"
@@ -3057 +3057 @@
              !             !
              DO ib =1, nbpt
-                soilclass(ib,1)=textrefrac(ib,1)
-                soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
-                soilclass(ib,3)=textrefrac(ib,5)
+                soilclass(ib,:)=0.
+                soilclass(ib,fao2usda(1))=textrefrac(ib,1)
+                soilclass(ib,fao2usda(2))=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
+                soilclass(ib,fao2usda(3))=textrefrac(ib,5)
 
                 ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
@@ -3075 +3076 @@
                      textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
 
-                sgn=SUM(soilclass(ib,1:3))
-
-                IF (sgn < min_sechiba) THEN
+                sgn=SUM(soilclass(ib,1:3)) ! grid-cell fraction with texture info
+
+                IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume 28%/52%, 20% of texture classes 3/6/9
                    soilclass(ib,:) = soilclass_default(:)
                    clayfraction(ib) = clayfraction_default
@@ -3135 +3136 @@
                 ENDDO
 
-                sgn=SUM(soilclass(ib,:))
-
-                IF (sgn < min_sechiba) THEN
+                sgn=SUM(soilclass(ib,:)) ! grid-cell fraction with texture info
+
+                IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume 28%/52%, 20% of texture classes 3/6/9
                    soilclass(ib,:) = soilclass_default(:)
                    clayfraction(ib) = clayfraction_default
@@ -3234 +3235 @@
              IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
              DO ib=1, nbpt
-                soilclass(ib,:) = soilclass_default_fao
+                soilclass(ib,:) = soilclass_default_usda
                 clayfraction(ib) = clayfraction_default
                 sandfraction(ib) = sandfraction_default
@@ -3241 +3242 @@
           CASE('zobler')
              !
-             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
+             soilclass(ib,:) = soilclass_default_usda ! USDA means here 13 final texture classes, owing to fao2usda
              !
              IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
@@ -3294 +3295 @@
                          SELECT CASE(solt(ilf))
                          CASE(1)
-                            soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(1)) = soilclass(ib,fao2usda(1)) + textrefrac(ib,solt(ilf))
                          CASE(2)
-                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
                          CASE(3)
-                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
                          CASE(4)
-                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
                          CASE(5)
-                            soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(3)) = soilclass(ib,fao2usda(3)) + textrefrac(ib,solt(ilf))
                          CASE(7)
-                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                            soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
                          CASE DEFAULT
                             WRITE(numout,*) 'We should not be here, an impossible case appeared'
@@ -3451 +3452 @@
        njsc(:) = 0
        DO ib = 1, nbpt
-          njsc(ib) = MAXLOC(soilclass(ib,:),1)
+          njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Here we get 3/6/9 for the Zobler classes Coarse/Medium/Fine
        ENDDO
 
@@ -3545 +3546 @@
 
        ELSE ! spmipexp is not maps nor unif, then it must be texture
-         IF (spmipexp == 'texture') THEN
-            IF (soil_classif == 'usda') THEN
-               ! Texture map from SP-MIP, thus Soilgrids modified
-               nvan(:) = nvan_usda(njsc(:))
-               avan(:) = avan_usda(njsc(:))
-               mcr(:) = mcr_usda(njsc(:))
-               mcs(:) = mcs_usda(njsc(:))
-               ks(:) = ks_usda(njsc(:))
-               mcfc(:) = mcf_usda(njsc(:))
-               mcw(:) = mcw_usda(njsc(:))
-               ! on aura pcent(:) = pcent(njsc(:)) dans hydrol
-            ELSE ! soil_classif == 'zobler' or 'none'
-               ! salma: here we are in exp3 -- Zobler map
-               nvan(:) = nvan_fao(njsc(:))
-               avan(:) = avan_fao(njsc(:))
-               mcr(:) = mcr_fao(njsc(:))
-               mcs(:) = mcs_fao(njsc(:))
-               ks(:) = ks_fao(njsc(:))
-               mcfc(:) = mcf_fao(njsc(:))
-               mcw(:) = mcw_fao(njsc(:))
-            ENDIF !if spmipexp is texture       
+          IF (spmipexp == 'texture') THEN
+             ! Whichever the soil texture map, we can use the USDA parameter vectors with 13 values 
+             nvan(:) = nvan_usda(njsc(:))
+             avan(:) = avan_usda(njsc(:))
+             mcr(:) = mcr_usda(njsc(:))
+             mcs(:) = mcs_usda(njsc(:))
+             ks(:) = ks_usda(njsc(:))
+             mcfc(:) = mcf_usda(njsc(:))
+             mcw(:) = mcw_usda(njsc(:))    
          ELSE ! if spmipexp is not among texture or maps or unif
             WRITE(numout,*) "Unsupported spmipexp=",spmipexp

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90

@@ -379 +379 @@
 
 
-!! Soil texture choose
-    SELECTCASE (nscm)
-    CASE (3)
-       QZ(:) = QZ_fao(:)
-       so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
-    CASE (13) !Salma changed from 12 to 13 for the new class Oxisols
-       QZ(:) = QZ_usda(:)
-       so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
-    CASE DEFAULT
-       WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
-       STOP 'thermosoil_initialize'
-    ENDSELECT
-
+    !! Soil texture choose : Now useless since njsc defines the dominant texture within 13 classes whichever the soil map
+    QZ(:) = QZ_usda(:)
+    so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
     
     !! 2. Initialize variable from restart file or with default values 
@@ -1330 +1320 @@
           ! Eq 11 in Peters-Lidard et al., 1998
           IF ( satratio(ji,jg) >  0.1 ) THEN
-            IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') )  THEN
+            IF (jst < 4 )  THEN
                 ! Coarse 
                 ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0
@@ -1338 +1328 @@
             ENDIF
           ELSEIF ( satratio(ji,jg) >  0.05 .AND. satratio(ji,jg) <=  0.1 ) THEN
-            IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') )  THEN
+            IF (jst < 4 )  THEN
                 ! Coarse 
                 ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0