Changeset 7339

Timestamp:

2021-11-06T13:15:01+01:00 (23 months ago)

Author:

agnes.ducharne

Message:

Contributes to ticket #416: field capacity and wilting point are now calculated in slowproc_soilt, instead of being defined in constantes_soil_var for each texture class. The latter values had been calculated with the same functions as in slowproc_soilt, but truncated, so the results are not perfectly preserved: on global scale average, FC is reduced by 0.005%, WP by 0.02%, and total ET by 0.02% (mean for 1 day = Jan 5th).

File:

• branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90 (modified) (4 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

@@ -2833 +2833 @@
     INTEGER(i_std), DIMENSION(:), ALLOCATABLE            :: vecpos
     REAL(r_std)                                          :: sgn              !! sum of fractions excluding glaciers and ocean
+
+    ! For the calculation of field capacity and wilting point
+    REAL(r_std),DIMENSION (nbpt)                         :: mvan             !! Van Genuchten parameter m
+    REAL(r_std),DIMENSION (nbpt)                         :: psi_w            !! Matrix potential characterizing the wilting point (mm)
+    REAL(r_std),DIMENSION (nbpt)                         :: psi_fc           !! Matrix potential characterizing the field capacity (mm)
+    
 !_ ================================================================================================================================
 
@@ -3466 +3472 @@
               WRITE(numout,*) 'thetas map read _______'
 
-              variablename = 'thetapwpvg'
+              variablename = 'thetapwpvg' ! mcw
               CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
                    contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
@@ -3473 +3479 @@
               WRITE(numout,*) 'thetapwpvg map read _______'
 
-              variablename = 'thetafcvg'
+              variablename = 'thetafcvg' !mcfc
               CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
                    contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
@@ -3495 +3501 @@
              mcs(:) = mcs_usda(njsc(:))
              ks(:) = ks_usda(njsc(:))
-             mcfc(:) = mcf_usda(njsc(:))
-             mcw(:) = mcw_usda(njsc(:))    
+!!$             mcfc(:) = mcf_usda(njsc(:))
+!!$             mcw(:) = mcw_usda(njsc(:))
+             
+             !! Calculation of FC and WP based on above 5 parameters
+             mvan(:) = un - (un / nvan(:))
+             ! Define matrix potential in mm for wilting point and field capacity (with sand vs clay-silt variation)
+             psi_w(:) = 150000.
+             DO ib=1, nbpt
+                IF ( ks(ib) .GE. 560 ) THEN ! Sandy soils (560 is equivalent of 2.75 at log scale of Ks, mm/d)
+                   psi_fc(ib) = 1000.
+                ELSE ! Finer soils
+                   psi_fc(ib) = 3300. 
+                ENDIF
+             ENDDO
+             mcfc(:) = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) * psi_fc(:))** nvan(:))** mvan(:))
+             mcw(:)  = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) *  psi_w(:))** nvan(:))** mvan(:))
+             
          ELSE ! if spmipexp is not among texture or maps or unif
             WRITE(numout,*) "Unsupported spmipexp=",spmipexp