|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE ICE_SOILDIF(PTSTEP, PTAUICE, PKSFC_IVEG, PLEGI, & 00003 PSOILHCAPZ, PWSATZ, PMPOTSATZ, PBCOEFZ, & 00004 PTG, PWGI, PWG, KWG_LAYER, & 00005 PDELTAT, PDZG, PWGI_EXCESS ) 00006 ! ########################################################################## 00007 ! 00008 !!**** *ICE_SOILDIF* 00009 ! 00010 !! PURPOSE 00011 !! ------- 00012 ! This subroutine calculates soil water phase changes using the 00013 ! available/excess energy approach. Soil temperature and volumetric 00014 ! ice content are adjusted due to phase changes. See the references 00015 ! Boone et al., 39, JAM, 2000 and Giard and Bazile, 128, MWR, 2000. 00016 ! NOTE that more recently a modification was made: freeze/thaw follows 00017 ! a relationship between liquid water and temperature derriving from 00018 ! the Clausius Clapeyron Eq. This results in little to no freezing for 00019 ! sufficiently dry but cold (below freezing) soils. Scatter about this 00020 ! curve results due to 'phase change efficiencies' and the surface insolation 00021 ! coefficient. 00022 ! 00023 !!** METHOD 00024 !! ------ 00025 ! 00026 ! Direct calculation 00027 ! 00028 !! EXTERNAL 00029 !! -------- 00030 ! 00031 ! None 00032 !! 00033 !! IMPLICIT ARGUMENTS 00034 !! ------------------ 00035 !! 00036 !! 00037 !! REFERENCE 00038 !! --------- 00039 !! 00040 !! Boone (2000) 00041 !! Boone et al., (2000) 00042 !! 00043 !! AUTHOR 00044 !! ------ 00045 !! A. Boone * Meteo-France * 00046 !! 00047 !! MODIFICATIONS 00048 !! ------------- 00049 !! Original 28/02/00 Boone 00050 !! Modified 24/11/09 Boone 00051 !! Limit energy available for phase change by 00052 ! local amount owing to diffusion. Has almost 00053 ! no impact except under rare circumstances 00054 ! (avoids rare but possible oscillatory behavior) 00055 ! Also, add minimum (numerical) melt/freeze efficieny to prevent 00056 ! prolonged periods of small ice amounts approaching zero. 00057 ! 00058 !! Modified 01/06/11 Boone 00059 ! Use apparent heat capacity linearization for freezing 00060 ! (when temperature depenence on ice change is direct) 00061 ! Do away with efficiency coefficients as they acted to provide 00062 ! numerical stability: not needed as apparent heat capacity increases 00063 ! the effective heat capacity thus increasing stability. 00064 ! NOTE: for now considers Brooks & Corey type water retention. 00065 ! 00066 ! Modified 08/2011 Decharme 00067 ! Optimization 00068 !------------------------------------------------------------------------------- 00069 ! 00070 !* 0. DECLARATIONS 00071 ! ------------ 00072 ! 00073 USE MODD_CSTS, ONLY : XLMTT, XTT, XG, XCI, XRHOLI, XRHOLW, XLSTT 00074 USE MODD_ISBA_PAR, ONLY : XWGMIN 00075 ! 00076 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00077 USE PARKIND1 ,ONLY : JPRB 00078 ! 00079 IMPLICIT NONE 00080 ! 00081 !* 0.1 declarations of arguments 00082 ! 00083 REAL, INTENT(IN) :: PTSTEP ! Model time step (s) 00084 ! 00085 REAL, DIMENSION(:), INTENT(IN) :: PTAUICE, PKSFC_IVEG, PLEGI 00086 ! PKSFC_IVEG = effect of surface layer insolation on phase changes 00087 ! Giard and Bazile (2000): non-dimensional 00088 ! PTAUICE = soil phase change characteristic time scale (s) 00089 ! PLEGI = ice sublimation 00090 ! 00091 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ, PWSATZ, PDELTAT 00092 ! PSOILHCAPZ = soil heat capacity [J/(m3 K)] 00093 ! PWSATZ = soil porosity (m3/m3) 00094 ! PDELTAT = change in temperature over the time 00095 ! step before adjustment owing to phase 00096 ! changes (K) 00097 ! 00098 REAL, DIMENSION(:,:), INTENT(IN) :: PDZG 00099 ! PDZG = Layer thickness (DIF option) 00100 ! 00101 REAL, DIMENSION(:,:), INTENT(IN) :: PMPOTSATZ, PBCOEFZ 00102 ! PMPOTSATZ = matric potential at saturation (m) 00103 ! PBCOEFZ = slope of the water retention curve (-) 00104 ! 00105 REAL, DIMENSION(:,:), INTENT(INOUT) :: PTG, PWGI, PWG 00106 ! PTG = soil temperature (K) 00107 ! PWGI = soil volumetric ice content (m3/m3) 00108 ! PWGI = soil volumetric liquid water content (m3/m3) 00109 ! 00110 INTEGER, DIMENSION(:), INTENT(IN) :: KWG_LAYER 00111 ! KWG_LAYER = Number of soil moisture layers (DIF option) 00112 ! 00113 REAL, DIMENSION(:), INTENT(OUT) :: PWGI_EXCESS 00114 ! PWGI_EXCESS = Soil ice excess water content 00115 ! 00116 !* 0.2 declarations of local variables 00117 ! 00118 INTEGER :: JJ, JL ! loop control 00119 ! 00120 INTEGER :: INI ! Number of point 00121 INTEGER :: INL ! Number of explicit soil layers 00122 INTEGER :: IDEPTH ! Total moisture soil depth 00123 ! 00124 REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2)) :: ZK, ZEXCESSFC 00125 ! 00126 REAL, DIMENSION(SIZE(PTG,1)) :: ZEXCESS 00127 ! 00128 REAL :: ZWGMAX, ZPSIMAX, ZPSI, ZDELTAT, 00129 ZPHASE, ZTGM, ZWGM, ZWGIM, ZLOG, 00130 ZEFFIC, ZPHASEM, ZPHASEF, ZWORK, 00131 ZAPPHEATCAP 00132 ! 00133 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00134 !------------------------------------------------------------------------------- 00135 ! 00136 ! Initialization: 00137 ! --------------- 00138 ! 00139 IF (LHOOK) CALL DR_HOOK('ICE_SOILDIF',0,ZHOOK_HANDLE) 00140 ! 00141 INI = SIZE(PTG(:,:),1) 00142 INL = MAXVAL(KWG_LAYER(:)) 00143 ! 00144 ZEXCESSFC (:,:)=0.0 00145 ZEXCESS (: )=0.0 00146 PWGI_EXCESS(: )=0.0 00147 ! 00148 !------------------------------------------------------------------------------- 00149 ! 00150 ! 1. Surface layer vegetation insulation coefficient (-) 00151 ! --------------------------------------------------- 00152 ! 00153 ZK(:,:) = 1.0 00154 ZK(:,1) = PKSFC_IVEG(:) 00155 ! 00156 ! 2. Soil ice evolution computation: 00157 ! ------------------------------- 00158 ! 00159 DO JL=1,INL 00160 DO JJ=1,INI 00161 IDEPTH=KWG_LAYER(JJ) 00162 IF(JL<=IDEPTH)THEN 00163 ! 00164 ZWGIM = PWGI(JJ,JL) 00165 ZWGM = PWG(JJ,JL) 00166 ZTGM = PTG(JJ,JL) 00167 00168 ! The maximum liquid water content as 00169 ! as function of temperature (sub-freezing) 00170 ! based on Gibbs free energy (Fuchs et al., 1978): 00171 ! 00172 ZPSIMAX = MIN(PMPOTSATZ(JJ,JL),XLMTT*(PTG(JJ,JL)-XTT)/(XG*PTG(JJ,JL))) 00173 ! 00174 ZWORK = ZPSIMAX/PMPOTSATZ(JJ,JL) 00175 ZLOG = LOG(ZWORK)/PBCOEFZ(JJ,JL) 00176 ZWGMAX = PWSATZ(JJ,JL)*EXP(-ZLOG) 00177 ! 00178 ! Calculate maximum temperature for ice based on Gibbs free energy: first 00179 ! compute soil water potential using Brook and Corey (1966) model: 00180 ! psi=mpotsat*(w/wsat)**(-bcoef) 00181 ! 00182 ZWORK = PWG(JJ,JL)/PWSATZ(JJ,JL) 00183 ZLOG = PBCOEFZ(JJ,JL)*LOG(ZWORK) 00184 ZPSI = PMPOTSATZ(JJ,JL)*EXP(-ZLOG) 00185 ! 00186 ZDELTAT = PTG(JJ,JL) - XLMTT*XTT/(XLMTT-XG*ZPSI) 00187 ! 00188 ! Compute apparent heat capacity. This is considered 00189 ! only when there is available energy (cold) and liquid water 00190 ! available...freezing front. 00191 ! This also has the secondary effect of increasing numerical stability 00192 ! during freezing (as there is a strong temperature dependence) by 00193 ! i) potentially significantly increasing the "apparent" heat capacity and 00194 ! ii) this part is also treated implicitly herein. 00195 ! 00196 ZWORK = (XCI*XRHOLI/(XLMTT*XRHOLW))*ZK(JJ,JL)*MAX(0.0,-ZDELTAT) 00197 ! 00198 ZAPPHEATCAP=0.0 00199 IF(ZDELTAT<0.0.AND.ZWGM>=ZWGMAX.AND.ZWORK>=MAX(0.0,ZWGM-ZWGMAX))THEN 00200 ZAPPHEATCAP = -(XTT*XRHOLW*XLMTT*XLMTT/XG)*ZWGMAX/(ZPSIMAX*PBCOEFZ(JJ,JL)*ZTGM*ZTGM) 00201 ENDIF 00202 ! 00203 ! *Melt* ice if energy and ice available: 00204 ZPHASEM = (PTSTEP/PTAUICE(JJ))*MIN(ZK(JJ,JL)*XCI*XRHOLI*MAX(0.0,ZDELTAT),ZWGIM*XLMTT*XRHOLW) 00205 ! 00206 ! *Freeze* liquid water if energy and water available: 00207 ZPHASEF = (PTSTEP/PTAUICE(JJ))*MIN(ZK(JJ,JL)*XCI*XRHOLI*MAX(0.0,-ZDELTAT),MAX(0.0,ZWGM-ZWGMAX)*XLMTT*XRHOLW) 00208 ! 00209 ! Update heat content if melting or freezing 00210 PTG(JJ,JL) = ZTGM + (ZPHASEF - ZPHASEM)/(PSOILHCAPZ(JJ,JL)+ZAPPHEATCAP) 00211 ! 00212 ! Get estimate of actual total phase change (J/m3) for equivalent soil water changes: 00213 ZPHASE = ZPHASEF - ZPHASEM - ZAPPHEATCAP*(PTG(JJ,JL)-ZTGM) 00214 ! 00215 ! Adjust ice and liquid water conents (m3/m3) accordingly : 00216 PWGI(JJ,JL) = ZWGIM + ZPHASE/(XLMTT*XRHOLW) 00217 PWG (JJ,JL) = ZWGM - ZPHASE/(XLMTT*XRHOLW) 00218 ! 00219 ENDIF 00220 ENDDO 00221 ENDDO 00222 ! 00223 ! 3. Adjust surface soil ice content for sublimation 00224 ! ----------------------------------------------- 00225 ! 00226 PWGI(:,1) = PWGI(:,1) - PLEGI(:)*PTSTEP/(XLSTT*XRHOLW*PDZG(:,1)) 00227 ! 00228 ! The remaining code in this block are merely constraints to ensure a highly 00229 ! accurate water budget: most of the time this code will not have any 00230 ! effect on the soil water profile. 00231 ! If sublimation causes all of the remaining ice to be extracted, remove 00232 ! some of the liquid (a correction): NOTE that latent heating already accounted 00233 ! for in sublimation term, so no need to alter soil temperature. 00234 ! 00235 ZEXCESS(:) = MAX(0.0, - PWGI(:,1)) 00236 PWG (:,1) = PWG (:,1) - ZEXCESS(:) 00237 PWGI(:,1) = PWGI(:,1) + ZEXCESS(:) 00238 ZEXCESSFC(:,1)= ZEXCESSFC(:,1) - ZEXCESS(:) 00239 ! 00240 ! 4. Prevent some possible problems 00241 ! ------------------------------ 00242 ! 00243 ! If sublimation is negative (condensation), make sure ice does not 00244 ! exceed maximum possible. If it does, then put excess ice into layer below: 00245 ! This correction should rarely if ever cause any ice accumulation in the 00246 ! sub-surface layer: this is especially true of deeper layers but it is 00247 ! accounted for none-the-less. 00248 ! 00249 DO JL=1,INL 00250 DO JJ=1,INI 00251 IDEPTH=KWG_LAYER(JJ) 00252 IF(JL<=IDEPTH)THEN 00253 ZEXCESS(JJ) = MAX(0.0, PWGI(JJ,JL) - (PWSATZ(JJ,JL)-XWGMIN) ) 00254 PWGI(JJ,JL) = PWGI(JJ,JL) - ZEXCESS(JJ) 00255 ZEXCESSFC(JJ,JL) = ZEXCESSFC(JJ,JL) + ZEXCESS(JJ) 00256 IF(JL<IDEPTH)THEN 00257 PWGI(JJ,JL+1) = PWGI(JJ,JL+1) + ZEXCESS(JJ)*(PDZG(JJ,JL)/PDZG(JJ,JL+1)) 00258 ZEXCESSFC(JJ,JL+1)= ZEXCESSFC(JJ,JL+1) - ZEXCESS(JJ)*(PDZG(JJ,JL)/PDZG(JJ,JL+1)) 00259 ELSE 00260 PWGI_EXCESS(JJ) = ZEXCESS(JJ)*PDZG(JJ,IDEPTH)*XRHOLW/PTSTEP 00261 ENDIF 00262 ENDIF 00263 ENDDO 00264 ENDDO 00265 ! 00266 ! Prevent keeping track of very small numbers for ice content: (melt it) 00267 ! and conserve energy: 00268 ! 00269 DO JL=1,INL 00270 DO JJ=1,INI 00271 IDEPTH=KWG_LAYER(JJ) 00272 IF(JL<=IDEPTH.AND.PWGI(JJ,JL)>0.0.AND.PWGI(JJ,JL)<1.0E-6)THEN 00273 PWG (JJ,JL) = PWG(JJ,JL) + PWGI(JJ,JL) 00274 ZEXCESSFC(JJ,JL) = ZEXCESSFC(JJ,JL) + PWGI(JJ,JL) 00275 PWGI (JJ,JL) = 0.0 00276 ENDIF 00277 PTG (JJ,JL) = PTG(JJ,JL) - ZEXCESSFC(JJ,JL)*XLMTT*XRHOLW/PSOILHCAPZ(JJ,JL) 00278 ENDDO 00279 ENDDO 00280 ! 00281 ! 00282 IF (LHOOK) CALL DR_HOOK('ICE_SOILDIF',1,ZHOOK_HANDLE) 00283 ! 00284 !------------------------------------------------------------------------------- 00285 ! 00286 END SUBROUTINE ICE_SOILDIF
1.8.0
