|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE SNOW3L(HSNOWRES, TPTIME, OGLACIER, HIMPLICIT_WIND, & 00003 PPEW_A_COEF, PPEW_B_COEF, & 00004 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, & 00005 PSNOWSWE,PSNOWRHO,PSNOWHEAT,PSNOWALB, & 00006 PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE, & 00007 PTSTEP,PPS,PSR,PRR,PPSN3L, & 00008 PTA,PTG,PSW_RAD,PQA,PVMOD,PLW_RAD, PRHOA, & 00009 PUREF,PEXNS,PEXNA,PDIRCOSZW, & 00010 PZREF,PZ0,PZ0EFF,PZ0H,PALB, & 00011 PSOILCOND,PD_G, & 00012 PSNOWLIQ,PSNOWTEMP,PSNOWDZ, & 00013 PTHRUFAL,PGRNDFLUX,PEVAPCOR,PRNSNOW,PHSNOW,PGFLUXSNOW, & 00014 PHPSNOW,PLES3L,PLEL3L,PEVAP,PRI, & 00015 PEMISNOW,PCDSNOW,PUSTAR,PCHSNOW,PSNOWHMASS, & 00016 PPERMSNOWFRAC, PZENITH, PXLAT, PXLON ) 00017 ! ########################################################################## 00018 ! 00019 !!**** *SNOW3L* 00020 !! 00021 !! PURPOSE 00022 !! ------- 00023 ! 00024 ! 3-Layer snow scheme option (Boone and Etchevers 1999) 00025 ! For shallow snow cover, Default method of Douville et al. (1995) 00026 ! used with this option: Model "turns on" when snow sufficiently 00027 ! deep/above some preset critical snow depth. 00028 ! 00029 ! 00030 !!** METHOD 00031 !! ------ 00032 ! 00033 ! Direct calculation 00034 ! 00035 !! EXTERNAL 00036 !! -------- 00037 ! 00038 ! None 00039 !! 00040 !! IMPLICIT ARGUMENTS 00041 !! ------------------ 00042 !! 00043 !! 00044 !! 00045 !! REFERENCE 00046 !! --------- 00047 !! 00048 !! ISBA-ES: Boone and Etchevers (2001) 00049 !! ISBA: Belair (1995) 00050 !! ISBA: Noilhan and Planton (1989) 00051 !! ISBA: Noilhan and Mahfouf (1996) 00052 !! 00053 !! AUTHOR 00054 !! ------ 00055 !! A. Boone * Meteo-France * 00056 !! 00057 !! MODIFICATIONS 00058 !! ------------- 00059 !! Original 7/99 00060 !! Modified by A.Boone 05/02 (code, not physics) 00061 !! Modified by A.Boone 11/04 i) maximum density limit imposed (although 00062 !! rarely if ever reached), ii) check to 00063 !! see if upermost layer completely sublimates 00064 !! during a timestep (as snowpack becomes vanishly 00065 !! thin), iii) impose maximum grain size limit 00066 !! in radiation transmission computation. 00067 !! 00068 !! Modified by B. Decharme (03/2009): Consistency with Arpege permanent 00069 !! snow/ice treatment (LGLACIER for alb) 00070 !! Modified by A. Boone (04/2010): Implicit coupling and replace Qsat and DQsat 00071 !! by Qsati and DQsati, respectively. 00072 !! Modified by E. Brun (08/2012): Mass conservation in SNOW3LEVAPGONE 00073 !! Modified by B. Decharme (08/2012): Loop optimization 00074 !! Modified by B. Decharme (09/2012): New wind implicitation 00075 !! Modified by E. Brun (10/2012): Bug in vertical snow redistribution 00076 !! 00077 !------------------------------------------------------------------------------- 00078 ! 00079 !* 0. DECLARATIONS 00080 ! ------------ 00081 ! 00082 USE MODD_CSTS, ONLY : XTT, XRHOLW, XLMTT, XCL 00083 ! 00084 USE MODE_SNOW3L 00085 ! 00086 USE MODD_TYPE_DATE_SURF, ONLY: DATE_TIME 00087 ! 00088 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00089 USE PARKIND1 ,ONLY : JPRB 00090 ! 00091 IMPLICIT NONE 00092 ! 00093 ! 00094 !* 0.1 declarations of arguments 00095 ! 00096 REAL, INTENT(IN) :: PTSTEP 00097 ! PTSTEP = time step of the integration 00098 TYPE(DATE_TIME), INTENT(IN) :: TPTIME ! current date and time 00099 ! 00100 CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES 00101 ! HSNOWRES = ISBA-SNOW3L turbulant exchange option 00102 ! 'DEF' = Default: Louis (ISBA: Noilhan and Mahfouf 1996) 00103 ! 'RIL' = Limit Richarson number under very stable 00104 ! conditions (currently testing) 00105 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 00106 ! initialise WGI=WSAT, 00107 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 00108 ! False = No specific treatment 00109 ! 00110 CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option 00111 ! ! 'OLD' = direct 00112 ! ! 'NEW' = Taylor serie, order 1 00113 ! 00114 REAL, DIMENSION(:), INTENT(IN) :: PPS, PTA, PSW_RAD, PQA, 00115 PVMOD, PLW_RAD, PSR, PRR 00116 ! PSW_RAD = incoming solar radiation (W/m2) 00117 ! PLW_RAD = atmospheric infrared radiation (W/m2) 00118 ! PRR = rain rate [kg/(m2 s)] 00119 ! PSR = snow rate (SWE) [kg/(m2 s)] 00120 ! PTA = atmospheric temperature at level za (K) 00121 ! PVMOD = modulus of the wind parallel to the orography (m/s) 00122 ! PPS = surface pressure 00123 ! PQA = atmospheric specific humidity 00124 ! at level za 00125 ! 00126 REAL, DIMENSION(:), INTENT(IN) :: PTG, PSOILCOND, PD_G, PPSN3L 00127 ! PTG = Surface soil temperature (effective 00128 ! temperature the of layer lying below snow) 00129 ! PSOILCOND = soil thermal conductivity [W/(m K)] 00130 ! PD_G = Assumed first soil layer thickness (m) 00131 ! Used to calculate ground/snow heat flux 00132 ! PPSN3L = snow fraction 00133 ! 00134 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PUREF, PEXNS, PEXNA, PDIRCOSZW, PRHOA, PZ0, PZ0EFF, 00135 PALB, PZ0H, PPERMSNOWFRAC 00136 ! PZ0EFF = roughness length for momentum 00137 ! PZ0 = grid box average roughness length 00138 ! PZ0H = grid box average roughness length for heat 00139 ! PZREF = reference height of the first 00140 ! atmospheric level 00141 ! PUREF = reference height of the wind 00142 ! PRHOA = air density 00143 ! PEXNS = Exner function at surface 00144 ! PEXNA = Exner function at lowest atmos level 00145 ! PDIRCOSZW = Cosinus of the angle between the 00146 ! normal to the surface and the vertical 00147 ! PALB = soil/vegetation albedo 00148 ! PPERMSNOWFRAC = fraction of permanet snow/ice 00149 ! 00150 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, 00151 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, 00152 PPEQ_B_COEF 00153 ! PPEW_A_COEF = wind coefficient (m2s/kg) 00154 ! PPEW_B_COEF = wind coefficient (m/s) 00155 ! PPET_A_COEF = A-air temperature coefficient 00156 ! PPET_B_COEF = B-air temperature coefficient 00157 ! PPEQ_A_COEF = A-air specific humidity coefficient 00158 ! PPEQ_B_COEF = B-air specific humidity coefficient 00159 ! 00160 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWALB 00161 ! PSNOWALB = Prognostic surface snow albedo 00162 ! (does not include anything but 00163 ! the actual snow cover) 00164 ! 00165 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWHEAT, PSNOWRHO, PSNOWSWE 00166 ! PSNOWHEAT = Snow layer(s) heat content (J/m2) 00167 ! PSNOWRHO = Snow layer(s) averaged density (kg/m3) 00168 ! PSNOWSWE = Snow layer(s) liquid Water Equivalent (SWE:kg m-2) 00169 ! 00170 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST 00171 ! PSNOWGRAN1 = Snow layers grain feature 1 00172 ! PSNOWGRAN2 = Snow layer grain feature 2 00173 ! PSNOWHIST = Snow layer grain historical 00174 ! parameter (only for non 00175 ! dendritic snow) 00176 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWAGE ! Snow grain age 00177 ! 00178 REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWLIQ, PSNOWTEMP, PSNOWDZ 00179 ! PSNOWLIQ = Snow layer(s) liquid water content (m) 00180 ! PSNOWTEMP = Snow layer(s) temperature (m) 00181 ! PSNOWDZ = Snow layer(s) thickness (m) 00182 ! 00183 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL, PGRNDFLUX, PEVAPCOR 00184 ! PTHRUFAL = rate that liquid water leaves snow pack: 00185 ! paritioned into soil infiltration/runoff 00186 ! by ISBA [kg/(m2 s)] 00187 ! PGRNDFLUX = soil/snow interface heat flux (W/m2) 00188 ! PEVAPCOR = evaporation/sublimation correction term: 00189 ! extract any evaporation exceeding the 00190 ! actual snow cover (as snow vanishes) 00191 ! and apply it as a surface soil water 00192 ! sink. [kg/(m2 s)] 00193 ! 00194 REAL, DIMENSION(:), INTENT(OUT) :: PRNSNOW, PHSNOW, PGFLUXSNOW, PLES3L, PLEL3L, 00195 PHPSNOW, PCDSNOW, PUSTAR, PEVAP 00196 ! PLES3L = evaporation heat flux from snow (W/m2) 00197 ! PLEL3L = sublimation (W/m2) 00198 ! PHPSNOW = heat release from rainfall (W/m2) 00199 ! PRNSNOW = net radiative flux from snow (W/m2) 00200 ! PHSNOW = sensible heat flux from snow (W/m2) 00201 ! PGFLUXSNOW = net heat flux from snow (W/m2) 00202 ! PCDSNOW = drag coefficient for momentum over snow 00203 ! PUSTAR = friction velocity over snow (m/s) 00204 ! PEVAP = total evaporative flux (kg/m2/s) 00205 ! 00206 REAL, DIMENSION(:), INTENT(OUT) :: PCHSNOW, PEMISNOW, PSNOWHMASS 00207 ! PEMISNOW = snow surface emissivity 00208 ! PCHSNOW = drag coefficient for heat over snow 00209 ! PSNOWHMASS = heat content change due to mass 00210 ! changes in snowpack (J/m2): for budget 00211 ! calculations only. 00212 ! 00213 REAL, DIMENSION(:), INTENT(OUT) :: PRI 00214 ! PRI = Ridcharson number 00215 ! 00216 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle 00217 REAL, DIMENSION(:), INTENT(IN) :: PXLAT,PXLON ! LAT/LON after packing 00218 ! 00219 !* 0.2 declarations of local variables 00220 ! 00221 INTEGER :: JJ, JI ! Loop control 00222 ! 00223 INTEGER :: INI ! number of point 00224 INTEGER :: INLVLS ! number of snow layers 00225 ! 00226 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWTEMP, ZSCAP, ZSNOWDZN, ZSCOND, 00227 ZRADSINK 00228 ! ZSNOWTEMP = Snow layer(s) averaged temperature (K) 00229 ! ZSCAP = Snow layer(s) heat capacity [J/(K m3)] 00230 ! ZSNOWDZN = Updated snow layer thicknesses (m) 00231 ! ZSCOND = Snow layer(s) thermal conducivity [W/(m K)] 00232 ! ZRADSINK = Snow solar Radiation source terms (W/m2) 00233 ! 00234 REAL, DIMENSION(SIZE(PTA)) :: ZSNOW, ZSFCFRZ, ZTSTERM1, ZTSTERM2, 00235 ZCT, ZRA, ZSNOWFLUX, ZSNOWTEMPO1 00236 ! ZSNOW = Total snow depth (m) 00237 ! ZCT = inverse of the product of snow heat capacity 00238 ! and layer thickness [(m2 K)/J] 00239 ! ZRA = Surface aerodynamic resistance 00240 ! ZTSTERM1,ZTSTERM2 = Surface energy budget coefficients 00241 ! ZSNOWFLUX = heat flux between 1st and 2nd snow layers: 00242 ! used during surface melting (W/m2) 00243 ! ZSNOWTEMPO1= value of uppermost snow temperature 00244 ! before time integration (K) 00245 ! 00246 LOGICAL, DIMENSION(SIZE(PTA)) :: GSFCMELT 00247 ! GSFCMELT = FLAG if surface melt is occurring, used 00248 ! for surface albedo calculation. 00249 ! 00250 REAL, DIMENSION(SIZE(PTA)) :: ZRSRA, ZDQSAT, ZQSAT, ZRADXS, ZMELTXS, ZLIQHEATXS, 00251 ZLWUPSNOW, ZGRNDFLUXO 00252 ! ZRSRA = air density over aerodynamic resistance 00253 ! ZDQSAT = derrivative of saturation specific humidity 00254 ! ZQSAT = saturation specific humidity 00255 ! ZRADXS = shortwave radiation absorbed by soil surface 00256 ! (for thin snow sover) (W m-2) 00257 ! ZMELTXS = excess energy for snowmelt for vanishingly 00258 ! thin snowpacks: used to heat underlying surface 00259 ! in order to maintain energy conservation (W m-2) 00260 ! ZLIQHEATXS = excess snowpack heating for vanishingly thin 00261 ! snow cover: add energy to snow/ground heat 00262 ! flux (W m-2) 00263 ! ZLWUPSNOW = upwelling longwave raaditive flux (W m-2) 00264 ! ZGRNDFLUXO= snow-ground flux before correction (W m-2) 00265 ! This is used simply to test if snow 00266 ! completely melts during a timestep. 00267 ! 00268 REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2_IC, ZTA_IC, ZQA_IC, 00269 ZPET_A_COEF_T, ZPEQ_A_COEF_T, ZPET_B_COEF_T, ZPEQ_B_COEF_T 00270 ! ZUSTAR2_IC = implicit lowest atmospheric level friction (m2/s2) 00271 ! ZTA_IC = implicit lowest atmospheric level air temperature 00272 ! ZQA_IC = implicit lowest atmospheric level specific humidity 00273 ! ZPET_A_COEF_T = transformed A-air temperature coefficient 00274 ! ZPET_B_COEF_T = transformed B-air temperature coefficient 00275 ! ZPEQ_A_COEF_T = transformed A-air specific humidity coefficient 00276 ! ZPEQ_B_COEF_T = transformed B-air specific humidity coefficient 00277 ! 00278 REAL, PARAMETER :: ZSNOWDZMIN = 0.0001 ! (m) 00279 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00280 ! ZSNOWDZMIN = minimum snow layer thickness for 00281 ! thermal calculations. Used to prevent 00282 ! numerical problems as snow becomes 00283 ! vanishingly thin. 00284 ! - - --------------------------------------------------- 00285 ! 00286 ! 0. Initialization 00287 ! -------------- 00288 ! NOTE that snow layer thickness is used throughout this code: SWE 00289 ! is only used to diagnose the thickness at the beginning of this routine 00290 ! and it is updated at the end of this routine. 00291 ! 00292 IF (LHOOK) CALL DR_HOOK('SNOW3L',0,ZHOOK_HANDLE) 00293 PSNOWDZ(:,:) = PSNOWSWE(:,:)/PSNOWRHO(:,:) 00294 ! 00295 INI = SIZE(PSNOWSWE(:,:),1) 00296 INLVLS = SIZE(PSNOWSWE(:,:),2) ! total snow layers 00297 ! 00298 ZUSTAR2_IC = 0.0 00299 ZTA_IC = 0.0 00300 ZQA_IC = 0.0 00301 ! 00302 !* 1. Snow total depth 00303 ! ---------------- 00304 ! 00305 ZSNOW(:) = 0. 00306 DO JJ=1,SIZE(PSNOWDZ(:,:),2) 00307 DO JI=1,INI 00308 ZSNOW(JI) = ZSNOW(JI) + PSNOWDZ(JI,JJ) ! m 00309 ENDDO 00310 END DO 00311 ! 00312 ! 00313 !* 2. Snowfall 00314 ! -------- 00315 ! Caluclate uppermost density and thickness changes due to snowfall, 00316 ! and add heat content of falling snow 00317 ! 00318 CALL SNOW3LFALL(PTSTEP,OGLACIER,PSR,PTA,PVMOD,ZSNOW,PSNOWRHO,PSNOWDZ, & 00319 PSNOWHEAT,PSNOWHMASS,PSNOWALB,PPERMSNOWFRAC) 00320 ! 00321 ! 00322 !* 3. Update grid/discretization 00323 ! -------------------------- 00324 ! Reset grid to conform to model specifications: 00325 ! 00326 CALL SNOW3LGRID(ZSNOWDZN,ZSNOW) 00327 ! 00328 ! Mass/Heat redistribution: 00329 ! 00330 CALL SNOW3LTRANSF(ZSNOW,PSNOWDZ,ZSNOWDZN,PSNOWRHO,PSNOWHEAT) 00331 ! 00332 ! 00333 !* 4. Liquid water content and snow temperature 00334 ! ----------------------------------------- 00335 ! 00336 ! First diagnose snow temperatures and liquid 00337 ! water portion of the snow from snow heat content: 00338 ! 00339 ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO) 00340 ! 00341 ZSNOWTEMP(:,:) = XTT + ( ((PSNOWHEAT(:,:)/PSNOWDZ(:,:)) & 00342 + XLMTT*PSNOWRHO(:,:))/ZSCAP(:,:) ) 00343 ! 00344 PSNOWLIQ(:,:) = MAX(0.0,ZSNOWTEMP(:,:)-XTT)*ZSCAP(:,:)* & 00345 PSNOWDZ(:,:)/(XLMTT*XRHOLW) 00346 ! 00347 ZSNOWTEMP(:,:) = MIN(XTT,ZSNOWTEMP(:,:)) 00348 ! 00349 ! 00350 !* 5. Snow Compaction 00351 ! --------------- 00352 ! Calculate snow density: compaction/aging: density increases 00353 ! 00354 CALL SNOW3LCOMPACTN(PTSTEP,PSNOWRHO,PSNOWDZ,ZSNOWTEMP,ZSNOW) 00355 ! 00356 ! Update snow heat content (J/m2): 00357 ! 00358 ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO) 00359 PSNOWHEAT(:,:) = PSNOWDZ(:,:)*( ZSCAP(:,:)*(ZSNOWTEMP(:,:)-XTT) & 00360 - XLMTT*PSNOWRHO(:,:) ) + XLMTT*XRHOLW*PSNOWLIQ(:,:) 00361 ! 00362 ! 00363 !* 6. Solar radiation transmission 00364 ! ----------------------------- 00365 ! 00366 ! Heat source (-sink) term due to shortwave 00367 ! radiation transmission within the snowpack: 00368 ! 00369 CALL SNOW3LRAD(ZSNOWDZMIN,PSW_RAD,PSNOWALB,PSNOWDZ,PSNOWRHO, & 00370 PALB,ZRADSINK,ZRADXS) 00371 ! 00372 ! 00373 !* 7. Heat transfer and surface energy budget 00374 ! --------------------------------------- 00375 ! Snow thermal conductivity: 00376 ! 00377 CALL SNOW3LTHRM(PSNOWRHO,ZSCOND,ZSNOWTEMP,PPS) 00378 ! 00379 ! Precipitation heating term: 00380 ! Rainfall renders it's heat to the snow when it enters 00381 ! the snowpack: 00382 ! NOTE: for now set to zero because we'd need to remove this heat from 00383 ! the atmosphere to conserve energy. This can be done, but should 00384 ! be tested within an atmos model first probably... 00385 ! 00386 !PHPSNOW(:) = PRR(:)*XCL*(MAX(XTT,PTA(:))-XTT) ! (W/m2) 00387 PHPSNOW(:) = 0.0 00388 ! 00389 ! Surface Energy Budget calculations using ISBA linearized form 00390 ! and standard ISBA turbulent transfer formulation 00391 ! 00392 CALL SNOW3LEBUD(HSNOWRES, HIMPLICIT_WIND, & 00393 PPEW_A_COEF, PPEW_B_COEF, & 00394 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, & 00395 ZSNOWDZMIN, & 00396 PZREF,ZSNOWTEMP(:,1),PSNOWRHO(:,1),PSNOWLIQ(:,1),ZSCAP(:,1), & 00397 ZSCOND(:,1),ZSCOND(:,2), & 00398 PUREF,PEXNS,PEXNA,PDIRCOSZW,PVMOD, & 00399 PLW_RAD,PSW_RAD,PTA,PQA,PPS,PTSTEP, & 00400 PSNOWDZ(:,1),PSNOWDZ(:,2),PSNOWALB,PZ0,PZ0EFF,PZ0H, & 00401 ZSFCFRZ,ZRADSINK(:,1),PHPSNOW, & 00402 ZCT,PEMISNOW,PRHOA,ZTSTERM1,ZTSTERM2,ZRA,PCDSNOW,PCHSNOW, & 00403 ZQSAT, ZDQSAT, ZRSRA, ZUSTAR2_IC, PRI, & 00404 ZPET_A_COEF_T,ZPEQ_A_COEF_T,ZPET_B_COEF_T,ZPEQ_B_COEF_T ) 00405 ! 00406 ! 00407 ! Heat transfer: simple diffusion along the thermal gradient 00408 ! 00409 ZSNOWTEMPO1(:) = ZSNOWTEMP(:,1) ! save surface snow temperature before update 00410 ! 00411 CALL SNOW3LSOLVT(PTSTEP,ZSNOWDZMIN,PSNOWDZ,ZSCOND,ZSCAP,PTG, & 00412 PSOILCOND,PD_G,ZRADSINK,ZCT,ZTSTERM1,ZTSTERM2, & 00413 ZPET_A_COEF_T,ZPEQ_A_COEF_T,ZPET_B_COEF_T,ZPEQ_B_COEF_T, & 00414 ZTA_IC,ZQA_IC,PGRNDFLUX,ZGRNDFLUXO,ZSNOWTEMP,ZSNOWFLUX ) 00415 ! 00416 ! 00417 !* 8. Surface fluxes 00418 ! -------------- 00419 ! 00420 CALL SNOW3LFLUX(ZSNOWTEMP(:,1),PSNOWDZ(:,1),PEXNS,PEXNA, & 00421 ZUSTAR2_IC, & 00422 PTSTEP,PSNOWALB,PSW_RAD,PEMISNOW,ZLWUPSNOW,PLW_RAD, & 00423 ZTA_IC,ZSFCFRZ,ZQA_IC,PHPSNOW, & 00424 ZSNOWTEMPO1,ZSNOWFLUX,ZCT,ZRADSINK(:,1), & 00425 ZQSAT,ZDQSAT,ZRSRA, & 00426 PRNSNOW,PHSNOW,PGFLUXSNOW,PLES3L,PLEL3L,PEVAP, & 00427 PUSTAR,GSFCMELT ) 00428 ! 00429 ! 00430 !* 9. Snow melt 00431 ! --------- 00432 ! 00433 ! First Test to see if snow pack vanishes during this time step: 00434 ! 00435 CALL SNOW3LGONE(PTSTEP,PPSN3L,PLEL3L,PLES3L,PSNOWRHO, & 00436 PSNOWHEAT,ZRADSINK(:,INLVLS),PEVAPCOR,PTHRUFAL,PGRNDFLUX, & 00437 PGFLUXSNOW,ZGRNDFLUXO,PSNOWDZ,PSNOWLIQ,ZSNOWTEMP,ZRADXS) 00438 ! 00439 ! Add radiation not absorbed by snow to soil/vegetation interface flux 00440 ! (for thin snowpacks): 00441 ! 00442 PGRNDFLUX(:) = PGRNDFLUX(:) + ZRADXS(:) 00443 ! 00444 ! For "normal" melt: transform excess heat content into snow liquid: 00445 ! 00446 CALL SNOW3LMELT(PTSTEP,ZSCAP,ZSNOWTEMP,PSNOWDZ,PSNOWRHO, & 00447 PSNOWLIQ,ZMELTXS) 00448 ! 00449 ! Add any excess heat for melting to underlying surface 00450 ! (for vanishingly thin snowpacks): 00451 ! 00452 PGRNDFLUX(:) = PGRNDFLUX(:) + ZMELTXS(:) 00453 ! 00454 ! 00455 !* 10. Snow water flow and refreezing 00456 ! ------------------------------ 00457 ! Liquid water vertical transfer and possible snowpack runoff 00458 ! And refreezing/freezing of meltwater/rainfall (ripening of the snow) 00459 ! 00460 CALL SNOW3LREFRZ(PTSTEP,PRR,PSNOWRHO,ZSNOWTEMP,PSNOWDZ,PSNOWLIQ,PTHRUFAL) 00461 ! 00462 ! 00463 !* 11. Snow Evaporation/Sublimation mass updates: 00464 ! ------------------------------------------ 00465 ! 00466 CALL SNOW3LEVAPN(PPSN3L,PLES3L,PLEL3L,PTSTEP,ZSNOWTEMP(:,1),PSNOWRHO(:,1), & 00467 PSNOWDZ(:,1),PSNOWLIQ(:,1),PEVAPCOR) 00468 ! 00469 ! If all snow in uppermost layer evaporates/sublimates, re-distribute 00470 ! grid (below could be evoked for vanishingly thin snowpacks): 00471 ! 00472 CALL SNOW3LEVAPGONE(PSNOWHEAT,PSNOWDZ,PSNOWRHO,ZSNOWTEMP,PSNOWLIQ) 00473 ! 00474 ! 00475 !* 12. Update surface albedo: 00476 ! ---------------------- 00477 ! Snow clear sky albedo: 00478 ! 00479 CALL SNOW3LALB(PSNOWALB,PTSTEP,PSNOWLIQ(:,1),PSNOWDZ(:,1),PSNOWRHO(:,1), & 00480 GSFCMELT,PSR,PPERMSNOWFRAC) 00481 ! 00482 ! 00483 !* 13. Update snow heat content: 00484 ! ------------------------- 00485 ! Update the heat content (variable stored each time step) 00486 ! using current snow temperature and liquid water content: 00487 ! 00488 ! First, make check to make sure heat content not too large 00489 ! (this can result due to signifigant heating of thin snowpacks): 00490 ! add any excess heat to ground flux: 00491 ! 00492 DO JJ=1,INLVLS 00493 DO JI=1,INI 00494 ZLIQHEATXS(JI) = MAX(0.0,PSNOWLIQ(JI,JJ)*XRHOLW-0.10*PSNOWDZ(JI,JJ)*PSNOWRHO(JI,JJ))*XLMTT/PTSTEP 00495 PSNOWLIQ(JI,JJ)= PSNOWLIQ(JI,JJ) - ZLIQHEATXS(JI)*PTSTEP/(XRHOLW*XLMTT) 00496 PSNOWLIQ(JI,JJ)= MAX(0.0, PSNOWLIQ(JI,JJ)) 00497 PGRNDFLUX(JI) = PGRNDFLUX(JI) + ZLIQHEATXS(JI) 00498 ENDDO 00499 ENDDO 00500 ! 00501 PSNOWTEMP(:,:) = ZSNOWTEMP(:,:) 00502 ! 00503 ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO) 00504 ! 00505 PSNOWHEAT(:,:) = PSNOWDZ(:,:)*( ZSCAP(:,:)*(PSNOWTEMP(:,:)-XTT) & 00506 - XLMTT*PSNOWRHO(:,:) ) + XLMTT*XRHOLW*PSNOWLIQ(:,:) 00507 ! 00508 PSNOWSWE(:,:) = PSNOWDZ(:,:)*PSNOWRHO(:,:) 00509 ! 00510 ! 00511 IF (LHOOK) CALL DR_HOOK('SNOW3L',1,ZHOOK_HANDLE) 00512 CONTAINS 00513 ! 00514 ! 00515 ! 00516 !#################################################################### 00517 !#################################################################### 00518 !#################################################################### 00519 SUBROUTINE SNOW3LFALL(PTSTEP,OGLACIER,PSR,PTA,PVMOD,PSNOW,PSNOWRHO,PSNOWDZ, & 00520 PSNOWHEAT,PSNOWHMASS,PSNOWALB,PPERMSNOWFRAC) 00521 ! 00522 !! PURPOSE 00523 !! ------- 00524 ! Calculate changes to snowpack resulting from snowfall. 00525 ! Update mass and heat content of uppermost layer. 00526 ! 00527 ! 00528 USE MODD_CSTS, ONLY : XLMTT, XTT, XCI 00529 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN, XANSMAX, XAGLAMAX 00530 ! 00531 USE MODE_SNOW3L 00532 ! 00533 IMPLICIT NONE 00534 ! 00535 !* 0.1 declarations of arguments 00536 ! 00537 REAL, INTENT(IN) :: PTSTEP 00538 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 00539 ! initialise WGI=WSAT, 00540 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 00541 ! False = No specific treatment 00542 ! 00543 REAL, DIMENSION(:), INTENT(IN) :: PSR, PTA, PVMOD, PPERMSNOWFRAC 00544 ! 00545 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOW, PSNOWALB 00546 ! 00547 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ, PSNOWHEAT 00548 ! 00549 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWHMASS 00550 ! 00551 ! 00552 !* 0.2 declarations of local variables 00553 ! 00554 INTEGER :: JJ, JI 00555 ! 00556 INTEGER :: INI 00557 INTEGER :: INLVLS 00558 ! 00559 REAL, DIMENSION(SIZE(PTA)) :: ZSNOWFALL, ZRHOSNEW, 00560 ZSNOW, ZSNOWTEMP, 00561 ZSNOWFALL_DELTA, ZSCAP, ZANSMAX 00562 ! 00563 ! ISBA-ES CROCUS (Pahaut 1976): snowfall density coefficients: 00564 ! 00565 REAL, PARAMETER :: ZSNOWFALL_A_SN = 109.0 ! kg/m3 00566 REAL, PARAMETER :: ZSNOWFALL_B_SN = 6.0 ! kg/(m3 K) 00567 REAL, PARAMETER :: ZSNOWFALL_C_SN = 26.0 ! kg/(m7/2 s1/2) 00568 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00569 ! 00570 !------------------------------------------------------------------------------- 00571 ! 00572 ! 0. Initialize: 00573 ! ------------------ 00574 ! 00575 IF (LHOOK) CALL DR_HOOK('SNOW3LFALL',0,ZHOOK_HANDLE) 00576 ZRHOSNEW(:) = XRHOSMIN_ES 00577 ZSNOWFALL(:) = 0.0 00578 ZSCAP(:) = 0.0 00579 ZSNOW(:) = PSNOW(:) 00580 INI = SIZE(PSNOWDZ(:,:),1) 00581 INLVLS = SIZE(PSNOWDZ(:,:),2) 00582 ! 00583 PSNOWHMASS(:) = 0.0 00584 ! 00585 ! 00586 ! 1. Incorporate snowfall into snowpack: 00587 ! -------------------------------------- 00588 ! 00589 ! 00590 ! Heat content of newly fallen snow (J/m2): 00591 ! NOTE for now we assume the snowfall has 00592 ! the temperature of the snow surface upon reaching the snow. 00593 ! This is done as opposed to using the air temperature since 00594 ! this flux is quite small and has little to no impact 00595 ! on the time scales of interest. If we use the above assumption 00596 ! then, then the snowfall advective heat flux is zero. 00597 ! 00598 ZSNOWTEMP(:) = XTT 00599 ! 00600 WHERE (PSR(:) > 0.0 .AND. PSNOWDZ(:,1)>0.) 00601 ZSCAP(:) = SNOW3LSCAP(PSNOWRHO(:,1)) 00602 ZSNOWTEMP(:) = XTT + (PSNOWHEAT(:,1) + & 00603 XLMTT*PSNOWRHO(:,1)*PSNOWDZ(:,1))/ & 00604 (ZSCAP(:)*MAX(XSNOWDMIN/INLVLS,PSNOWDZ(:,1))) 00605 ZSNOWTEMP(:) = MIN(XTT, ZSNOWTEMP(:)) 00606 END WHERE 00607 ! 00608 WHERE (PSR(:) > 0.0) 00609 ! 00610 PSNOWHMASS(:) = PSR(:)*(XCI*(ZSNOWTEMP(:)-XTT)-XLMTT)*PTSTEP 00611 ! 00612 ! Snowfall density: Following CROCUS (Pahaut 1976) 00613 ! 00614 ZRHOSNEW(:) = MAX(XRHOSMIN_ES, ZSNOWFALL_A_SN + ZSNOWFALL_B_SN*(PTA(:)-XTT)+ & 00615 ZSNOWFALL_C_SN*SQRT(PVMOD(:))) 00616 ! 00617 ! 00618 ! Augment total pack depth: 00619 ! 00620 ZSNOWFALL(:) = PSR(:)*PTSTEP/ZRHOSNEW(:) ! snowfall thickness (m) 00621 ! 00622 PSNOW(:) = PSNOW(:) + ZSNOWFALL(:) 00623 ! 00624 ! Fresh snowfall changes the snowpack 00625 ! density, increases the total liquid water 00626 ! equivalent: in uppermost snow layer: 00627 ! 00628 PSNOWRHO(:,1) = (PSNOWDZ(:,1)*PSNOWRHO(:,1) + ZSNOWFALL(:)*ZRHOSNEW(:))/ & 00629 (PSNOWDZ(:,1)+ZSNOWFALL(:)) 00630 PSNOWDZ(:,1) = PSNOWDZ(:,1) + ZSNOWFALL(:) 00631 ! 00632 ! Add energy of snowfall to snowpack: 00633 ! Update heat content (J/m2) (therefore the snow temperature 00634 ! and liquid content): 00635 ! 00636 PSNOWHEAT(:,1) = PSNOWHEAT(:,1) + PSNOWHMASS(:) 00637 ! 00638 END WHERE 00639 ! 00640 ! 00641 ! 2. Case of new snowfall on a previously snow-free surface: 00642 ! ---------------------------------------------------------- 00643 ! 00644 ! When snow first falls on a surface devoid of snow, 00645 ! redistribute the snow mass throughout the 3 layers: 00646 ! (temperature already set in the calling routine 00647 ! for this case), also, intialize the albedo: 00648 ! 00649 IF(OGLACIER)THEN 00650 ZANSMAX(:) = XAGLAMAX * PPERMSNOWFRAC(:) + XANSMAX * (1.0-PPERMSNOWFRAC(:)) 00651 ELSE 00652 ZANSMAX(:) = XANSMAX 00653 ENDIF 00654 ! 00655 ZSNOWFALL_DELTA(:) = 0.0 00656 WHERE(ZSNOW(:) == 0.0 .AND. PSR(:) > 0.0) 00657 ZSNOWFALL_DELTA(:) = 1.0 00658 PSNOWALB(:) = ZANSMAX(:) 00659 END WHERE 00660 ! 00661 DO JJ=1,INLVLS 00662 DO JI=1,INI 00663 PSNOWDZ(JI,JJ) = ZSNOWFALL_DELTA(JI)*(ZSNOWFALL(JI) /INLVLS) + & 00664 (1.0-ZSNOWFALL_DELTA(JI))*PSNOWDZ(JI,JJ) 00665 ! 00666 PSNOWHEAT(JI,JJ) = ZSNOWFALL_DELTA(JI)*(PSNOWHMASS(JI)/INLVLS) + & 00667 (1.0-ZSNOWFALL_DELTA(JI))*PSNOWHEAT(JI,JJ) 00668 ! 00669 PSNOWRHO(JI,JJ) = ZSNOWFALL_DELTA(JI)*ZRHOSNEW(JI) + & 00670 (1.0-ZSNOWFALL_DELTA(JI))*PSNOWRHO(JI,JJ) 00671 ENDDO 00672 ENDDO 00673 ! 00674 IF (LHOOK) CALL DR_HOOK('SNOW3LFALL',1,ZHOOK_HANDLE) 00675 ! 00676 ! 00677 END SUBROUTINE SNOW3LFALL 00678 !#################################################################### 00679 !#################################################################### 00680 !#################################################################### 00681 SUBROUTINE SNOW3LCOMPACTN(PTSTEP,PSNOWRHO,PSNOWDZ, & 00682 PSNOWTEMP,PSNOW ) 00683 ! 00684 !! PURPOSE 00685 !! ------- 00686 ! Snow compaction due to overburden and settling. 00687 ! Mass is unchanged: layer thickness is reduced 00688 ! in proportion to density increases. Method 00689 ! of Anderson (1976): see Loth and Graf, 1993, 00690 ! J. of Geophys. Res., 98, 10,451-10,464. 00691 ! 00692 ! 00693 USE MODD_CSTS, ONLY : XTT, XG 00694 USE MODD_SNOW_PAR, ONLY : XRHOSMAX_ES 00695 ! 00696 IMPLICIT NONE 00697 ! 00698 !* 0.1 declarations of arguments 00699 ! 00700 REAL, INTENT(IN) :: PTSTEP 00701 ! 00702 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP 00703 ! 00704 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ 00705 ! 00706 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW 00707 ! 00708 ! 00709 !* 0.2 declarations of local variables 00710 ! 00711 INTEGER :: JJ, JI 00712 ! 00713 INTEGER :: INI 00714 INTEGER :: INLVLS 00715 ! 00716 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHO2, ZSETTLE, 00717 ZVISCOCITY, ZSMASS, 00718 ZWSNOWDZ 00719 ! 00720 REAL, DIMENSION(SIZE(PSNOW)) :: ZSMASSC 00721 ! 00722 ! ISBA-ES Compaction/Settling Coefficients from Loth and Graf (1993) 00723 ! following Mellor (1964) and Anderson (1976): 00724 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 00725 ! 00726 REAL, PARAMETER :: ZSNOWCMPCT_RHOD = 150.0 ! (kg/m3) 00727 REAL, PARAMETER :: ZSNOWCMPCT_ACM = 2.8e-6 ! (1/s) 00728 REAL, PARAMETER :: ZSNOWCMPCT_BCM = 0.04 ! (1/K) 00729 REAL, PARAMETER :: ZSNOWCMPCT_CCM = 0.460 ! (m3/kg) 00730 REAL, PARAMETER :: ZSNOWCMPCT_V0 = 3.7e7 ! (Pa/s) 00731 REAL, PARAMETER :: ZSNOWCMPCT_VT = 0.081 ! (1/K) 00732 REAL, PARAMETER :: ZSNOWCMPCT_VR = 0.018 ! (m3/kg) 00733 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00734 ! 00735 !------------------------------------------------------------------------------- 00736 ! 00737 ! 0. Initialization: 00738 ! ------------------ 00739 ! 00740 IF (LHOOK) CALL DR_HOOK('SNOW3LCOMPACTN',0,ZHOOK_HANDLE) 00741 ZSNOWRHO2(:,:) = PSNOWRHO(:,:) 00742 ZSETTLE(:,:) = ZSNOWCMPCT_ACM 00743 ZVISCOCITY(:,:) = ZSNOWCMPCT_V0 00744 INI = SIZE(PSNOWDZ(:,:),1) 00745 INLVLS = SIZE(PSNOWDZ(:,:),2) 00746 ! 00747 ! 00748 ! 1. Cumulative snow mass (kg/m2): 00749 ! -------------------------------- 00750 ! 00751 ZWSNOWDZ(:,:) = PSNOWDZ(:,:)*PSNOWRHO(:,:) 00752 ZSMASSC(:) = 0.0 00753 DO JJ=1,INLVLS 00754 DO JI=1,INI 00755 ZSMASS(JI,JJ) = ZSMASSC(JI) + ZWSNOWDZ(JI,JJ) 00756 ZSMASSC(JI) = ZSMASSC(JI) + ZWSNOWDZ(JI,JJ) 00757 ENDDO 00758 ENDDO 00759 ! 00760 ! 2. Compaction/Settling 00761 ! ---------------------- 00762 ! Compaction/settling if density below upper limit 00763 ! (compaction is generally quite small above ~ 500 kg m-3): 00764 ! 00765 WHERE(PSNOWRHO(:,:) < XRHOSMAX_ES) 00766 ! 00767 ! First calculate settling due to freshly fallen snow: 00768 ! 00769 ZSETTLE(:,:) = ZSNOWCMPCT_ACM*EXP( & 00770 -ZSNOWCMPCT_BCM*(XTT-PSNOWTEMP(:,:)) & 00771 -ZSNOWCMPCT_CCM*MAX(0.0, & 00772 PSNOWRHO(:,:)-ZSNOWCMPCT_RHOD)) 00773 ! 00774 ! Snow viscocity: 00775 ! 00776 ZVISCOCITY(:,:) = ZSNOWCMPCT_V0*EXP( ZSNOWCMPCT_VT*(XTT-PSNOWTEMP(:,:)) + & 00777 ZSNOWCMPCT_VR*PSNOWRHO(:,:) ) 00778 ! 00779 ! 00780 ! Calculate snow density: compaction from weight/over-burden 00781 ! Anderson 1976 method: 00782 ! 00783 ! 00784 ZSNOWRHO2(:,:) = PSNOWRHO(:,:) + PSNOWRHO(:,:)*PTSTEP*( & 00785 (XG*ZSMASS(:,:)/ZVISCOCITY(:,:)) & 00786 + ZSETTLE(:,:) ) 00787 ! 00788 ! 00789 ! Conserve mass by decreasing grid thicknesses in response 00790 ! to density increases 00791 ! 00792 PSNOWDZ(:,:) = PSNOWDZ(:,:)*(PSNOWRHO(:,:)/ZSNOWRHO2(:,:)) 00793 ! 00794 ! 00795 END WHERE 00796 ! 00797 ! 3. Update total snow depth and density profile: 00798 ! ----------------------------------------------- 00799 ! 00800 ! Compaction/augmentation of total snowpack depth 00801 ! 00802 PSNOW(:) = 0. 00803 DO JJ=1,INLVLS 00804 DO JI=1,INI 00805 PSNOW(JI) = PSNOW(JI) + PSNOWDZ(JI,JJ) 00806 ENDDO 00807 ENDDO 00808 ! 00809 ! Update density (kg m-3): 00810 ! 00811 PSNOWRHO(:,:) = ZSNOWRHO2(:,:) 00812 IF (LHOOK) CALL DR_HOOK('SNOW3LCOMPACTN',1,ZHOOK_HANDLE) 00813 ! 00814 !------------------------------------------------------------------------------- 00815 ! 00816 END SUBROUTINE SNOW3LCOMPACTN 00817 !#################################################################### 00818 !#################################################################### 00819 !#################################################################### 00820 SUBROUTINE SNOW3LTRANSF(PSNOW,PSNOWDZ,PSNOWDZN, & 00821 PSNOWRHO,PSNOWHEAT ) 00822 ! 00823 !! PURPOSE 00824 !! ------- 00825 ! Snow mass,heat and characteristics redistibution in case of 00826 ! grid resizing. Total mass and heat content of the overall snowpack 00827 ! unchanged/conserved within this routine. 00828 ! Same method as in Crocus 00829 ! 00830 USE MODD_SNOW_PAR, ONLY : XSNOWCRITD 00831 ! 00832 IMPLICIT NONE 00833 ! 00834 ! 00835 !* 0.1 declarations of arguments 00836 ! 00837 REAL, DIMENSION(: ), INTENT(IN) :: PSNOW 00838 ! 00839 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZN 00840 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT 00841 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO 00842 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ 00843 ! 00844 !* 0.2 declarations of local variables 00845 ! 00846 INTEGER :: JI, JL, JLO 00847 ! 00848 INTEGER :: INI 00849 INTEGER :: INLVLS 00850 ! 00851 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHON 00852 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEATN 00853 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZTOP_NEW 00854 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZBOT_NEW 00855 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHOO 00856 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEATO 00857 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWDZO 00858 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZTOP_OLD 00859 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZBOT_OLD 00860 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEAN, ZSNOWHEAO 00861 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZMASTOTN, ZMASTOTO 00862 ! 00863 REAL, DIMENSION(SIZE(PSNOW)) :: ZPSNOW_OLD, ZPSNOW_NEW 00864 REAL, DIMENSION(SIZE(PSNOW)) :: ZSUMHEAT, ZSUMSWE, ZSNOWMIX_DELTA 00865 ! 00866 REAL :: ZPROPOR 00867 ! 00868 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00869 ! 00870 !------------------------------------------------------------------------------- 00871 ! 00872 ! 0. Initialization: 00873 ! ------------------ 00874 ! 00875 ! 00876 IF (LHOOK) CALL DR_HOOK('SNOW3LTRANSF',0,ZHOOK_HANDLE) 00877 ! 00878 INI = SIZE(PSNOWRHO,1) 00879 INLVLS = SIZE(PSNOWRHO,2) 00880 ! 00881 ZPSNOW_NEW(:) = 0.0 00882 ZPSNOW_OLD(:) = PSNOW(:) 00883 ! 00884 DO JL=1,INLVLS 00885 DO JI=1,INI 00886 ZPSNOW_NEW(JI)=ZPSNOW_NEW(JI)+PSNOWDZN(JI,JL) 00887 ENDDO 00888 ENDDO 00889 ! 00890 ! initialization of variables describing the initial snowpack 00891 ! 00892 ZSNOWDZO (:,:) = PSNOWDZ (:,:) 00893 ZSNOWRHOO (:,:) = PSNOWRHO (:,:) 00894 ZSNOWHEATO(:,:) = PSNOWHEAT(:,:) 00895 ! 00896 ! 1. Calculate vertical grid limits (m): 00897 ! -------------------------------------- 00898 ! 00899 ZSNOWZTOP_OLD(:,1) = ZPSNOW_OLD(:) 00900 ZSNOWZTOP_NEW(:,1) = ZPSNOW_NEW(:) 00901 ZSNOWZBOT_OLD(:,1) = ZSNOWZTOP_OLD(:,1)-ZSNOWDZO(:,1) 00902 ZSNOWZBOT_NEW(:,1) = ZSNOWZTOP_NEW(:,1)-PSNOWDZN(:,1) 00903 ! 00904 DO JL=2,INLVLS!-1 00905 DO JI=1,INI 00906 ZSNOWZTOP_OLD(JI,JL) = ZSNOWZBOT_OLD(JI,JL-1) 00907 ZSNOWZTOP_NEW(JI,JL) = ZSNOWZBOT_NEW(JI,JL-1) 00908 ZSNOWZBOT_OLD(JI,JL) = ZSNOWZTOP_OLD(JI,JL )-ZSNOWDZO(JI,JL) 00909 ZSNOWZBOT_NEW(JI,JL) = ZSNOWZTOP_NEW(JI,JL )-PSNOWDZN(JI,JL) 00910 ENDDO 00911 ENDDO 00912 !-------------------------------Ã supprimer------------------------------------------------ 00913 ! Check consistency 00914 DO JI=1,INI 00915 IF (ABS(ZSNOWZBOT_OLD(JI,INLVLS)) > 0.00001) write (*,*) 'Error bottom OLD, pt:',JI 00916 IF (ABS(ZSNOWZBOT_NEW(JI,INLVLS)) > 0.00001) write (*,*) 'Error bottom NEW, pt:',JI 00917 ENDDO 00918 !-------------------------------Ã supprimer------------------------------------------------ 00919 ZSNOWZBOT_OLD(:,INLVLS)=0.0 00920 ZSNOWZBOT_NEW(:,INLVLS)=0.0 00921 ! 00922 ! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing: 00923 ! -------------------------------------------------------------------- 00924 ! 00925 ! loop over the new snow layers 00926 ! Summ or avergage of the constituting quantities of the old snow layers 00927 ! which are totally or partially inserted in the new snow layer 00928 ! 00929 ZSNOWHEAN(:,:)=0.0 00930 ZMASTOTN (:,:)=0.0 00931 ! 00932 DO JL=1,INLVLS 00933 DO JLO=1, INLVLS 00934 DO JI=1,INI 00935 IF((ZSNOWZTOP_OLD(JI,JLO)>ZSNOWZBOT_NEW(JI,JL)).AND.(ZSNOWZBOT_OLD(JI,JLO)<ZSNOWZTOP_NEW(JI,JL)))THEN 00936 ! 00937 ZPROPOR = (MIN(ZSNOWZTOP_OLD(JI,JLO), ZSNOWZTOP_NEW(JI,JL)) & 00938 - MAX(ZSNOWZBOT_OLD(JI,JLO), ZSNOWZBOT_NEW(JI,JL)))& 00939 / ZSNOWDZO(JI,JLO) 00940 ! 00941 ZMASTOTN (JI,JL)=ZMASTOTN (JI,JL)+ZPROPOR*ZSNOWRHOO (JI,JLO)*ZSNOWDZO(JI,JLO) 00942 ZSNOWHEAN(JI,JL)=ZSNOWHEAN(JI,JL)+ZPROPOR*ZSNOWHEATO(JI,JLO) 00943 ! 00944 ENDIF 00945 ENDDO 00946 ENDDO 00947 ENDDO 00948 ! 00949 ! the new layer inherits from the weighted average properties of the old ones 00950 ! heat and mass 00951 ! 00952 ZSNOWHEATN(:,:)= ZSNOWHEAN(:,:) 00953 ZSNOWRHON (:,:)= ZMASTOTN (:,:)/PSNOWDZN(:,:) 00954 ! 00955 ! 4. Vanishing or very thin snowpack check: 00956 ! ----------------------------------------- 00957 ! 00958 ! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous): 00959 ! this avoids problems related to heat and mass exchange for 00960 ! thin layers during heavy snowfall or signifigant melt: one 00961 ! new/old layer can exceed the thickness of several old/new layers. 00962 ! Therefore, mix (conservative): 00963 ! 00964 ZSUMHEAT(:) = 0.0 00965 ZSUMSWE(:) = 0.0 00966 ZSNOWMIX_DELTA(:) = 0.0 00967 ! 00968 DO JL=1,INLVLS 00969 DO JI=1,INI 00970 IF(PSNOW(JI) < XSNOWCRITD)THEN 00971 ZSUMHEAT (JI) = ZSUMHEAT(JI) + PSNOWHEAT(JI,JL) 00972 ZSUMSWE (JI) = ZSUMSWE (JI) + PSNOWRHO (JI,JL)*PSNOWDZ(JI,JL) 00973 ZSNOWMIX_DELTA(JI) = 1.0 00974 ENDIF 00975 ENDDO 00976 ENDDO 00977 ! 00978 ! Heat and mass are evenly distributed vertically: 00979 ! heat and mass (density and thickness) are constant 00980 ! in profile: 00981 ! 00982 DO JL=1,INLVLS 00983 DO JI=1,INI 00984 ZSNOWHEATN(JI,JL) = ZSNOWMIX_DELTA(JI)*(ZSUMHEAT(JI)/INLVLS) + & 00985 (1.0-ZSNOWMIX_DELTA(JI))*ZSNOWHEATN(JI,JL) 00986 ! 00987 PSNOWDZN(JI,JL) = ZSNOWMIX_DELTA(JI)*(PSNOW(JI)/INLVLS) + & 00988 (1.0-ZSNOWMIX_DELTA(JI))*PSNOWDZN(JI,JL) 00989 ! 00990 ZSNOWRHON(JI,JL) = ZSNOWMIX_DELTA(JI)*(ZSUMSWE(JI)/PSNOW(JI)) + & 00991 (1.0-ZSNOWMIX_DELTA(JI))*ZSNOWRHON(JI,JL) 00992 ENDDO 00993 ENDDO 00994 ! 00995 !-------------------------------Ã supprimer------------------------------------------------ 00996 ! check of consistency between new and old snowpacks 00997 ! 00998 ZSNOWHEAN(:,:)=0. 00999 ZSNOWHEAO(:,:)=0. 01000 ZMASTOTN (:,:)=0. 01001 ZMASTOTO (:,:)=0. 01002 ! 01003 ZPSNOW_NEW(:)=0. 01004 ZPSNOW_OLD(:)=0. 01005 01006 DO JL=1,INLVLS 01007 DO JI=1,INI 01008 ! 01009 ZSNOWHEAN (JI,1) = ZSNOWHEAN (JI,1) + ZSNOWHEATN(JI,JL) 01010 ZMASTOTN (JI,1) = ZMASTOTN (JI,1) + ZSNOWRHON (JI,JL)*PSNOWDZN(JI,JL) 01011 ZPSNOW_NEW(JI ) = ZPSNOW_NEW(JI ) + PSNOWDZN (JI,JL) 01012 ! 01013 ZSNOWHEAO (JI,1) =ZSNOWHEAO (JI,1) + ZSNOWHEATO(JI,JL) 01014 ZMASTOTO (JI,1) =ZMASTOTO (JI,1) + ZSNOWRHOO (JI,JL)*ZSNOWDZO(JI,JL) 01015 ZPSNOW_OLD(JI ) =ZPSNOW_OLD(JI ) + ZSNOWDZO (JI,JL) 01016 ! 01017 ENDDO 01018 ENDDO 01019 DO JI=1,INI 01020 if (abs(ZSNOWHEAN (JI,1)-ZSNOWHEAO (JI,1))>0.0001.OR. & 01021 ABS(ZMASTOTN (JI,1)-ZMASTOTO (JI,1))>0.0001.OR. & 01022 ABS(ZPSNOW_NEW(JI )-ZPSNOW_OLD(JI ))> 0.0001 )THEN 01023 write(*,*) 'pt:',JI,'Warning diff', ZSNOWHEAN(JI,1)-ZSNOWHEAO(JI,1),ZMASTOTN(JI,1)-ZMASTOTO(JI,1),ZPSNOW_NEW(JI)-ZPSNOW_OLD(JI) 01024 ENDIF 01025 ENDDO 01026 !-------------------------------Ã supprimer------------------------------------------------ 01027 ! 01028 ! 01029 ! 5. Update mass (density and thickness) and heat: 01030 ! ------------------------------------------------ 01031 ! 01032 PSNOWDZ (:,:) = PSNOWDZN (:,:) 01033 PSNOWRHO (:,:) = ZSNOWRHON (:,:) 01034 PSNOWHEAT(:,:) = ZSNOWHEATN(:,:) 01035 ! 01036 IF (LHOOK) CALL DR_HOOK('SNOW3LTRANSF',1,ZHOOK_HANDLE) 01037 ! 01038 ! 01039 !------------------------------------------------------------------------------- 01040 ! 01041 END SUBROUTINE SNOW3LTRANSF 01042 !#################################################################### 01043 !#################################################################### 01044 !#################################################################### 01045 SUBROUTINE SNOW3LALB(PALBEDOSC,PTSTEP,PSNOWLIQ,PSNOWDZ, & 01046 PSNOWRHO,OSFCMELT,PSR,PPERMSNOWFRAC) 01047 ! 01048 !! PURPOSE 01049 !! ------- 01050 ! Calculate the snow surface albedo. Use the method of 01051 ! Douville et al. 1995, Clim. Dyn., 12, 21-35. but instead 01052 ! of using whether or not there is snowmelt to determine 01053 ! albedo decrease rate, use a continuous formualtion in which 01054 ! weights are calculated as a function of the degree of saturation of the 01055 ! uppermost layer (with respect to liquid water content). 01056 ! 01057 ! 01058 USE MODD_CSTS, ONLY : XDAY 01059 USE MODD_SNOW_PAR, ONLY : XWCRN, XANSMAX, XANSMIN, XANS_TODRY, & 01060 XSNOWDMIN, XANS_T, XAGLAMIN, XAGLAMAX 01061 ! 01062 ! 01063 USE MODE_SNOW3L 01064 ! 01065 IMPLICIT NONE 01066 ! 01067 !* 0.1 declarations of arguments 01068 ! 01069 REAL, INTENT(IN) :: PTSTEP 01070 ! 01071 REAL, DIMENSION(:), INTENT(IN) :: PSNOWLIQ,PSNOWDZ,PSNOWRHO,PSR,PPERMSNOWFRAC 01072 ! 01073 LOGICAL, DIMENSION(:), INTENT(IN) :: OSFCMELT 01074 ! 01075 REAL, DIMENSION(:), INTENT(INOUT) :: PALBEDOSC 01076 ! 01077 !* 0.2 declarations of local variables 01078 ! 01079 REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZPCPALB, ZALBDRY, ZANSMIN, ZANSMAX, 01080 ZALBWET, ZWHOLDMAX, ZFRACLIQ 01081 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01082 ! 01083 !------------------------------------------------------------------------------- 01084 ! 01085 ! 0. Initialize: 01086 ! ------------------ 01087 ! 01088 IF (LHOOK) CALL DR_HOOK('SNOW3LALB',0,ZHOOK_HANDLE) 01089 ZFRACLIQ(:) = 0.0 01090 ZALBWET(:) = 0.0 01091 ZWHOLDMAX(:) = 0.0 01092 ! 01093 IF(OGLACIER)THEN 01094 ZANSMIN(:) = XAGLAMIN * PPERMSNOWFRAC(:) + XANSMIN * (1.0-PPERMSNOWFRAC(:)) 01095 ZANSMAX(:) = XAGLAMAX * PPERMSNOWFRAC(:) + XANSMAX * (1.0-PPERMSNOWFRAC(:)) 01096 ELSE 01097 ZANSMIN(:) = XANSMIN 01098 ZANSMAX(:) = XANSMAX 01099 ENDIF 01100 ! 01101 ! 01102 ! 1. Albedo change (increase) due to snowfall: 01103 ! -------------------------------------------- 01104 ! 01105 ZPCPALB(:) = (PSR(:)*PTSTEP/XWCRN)*(ZANSMAX(:) - ZANSMIN(:)) 01106 ! 01107 ! 2. dry snow albedo rate of change: 01108 ! ---------------------------------- 01109 ! 01110 ZALBDRY(:) = PALBEDOSC(:) - XANS_TODRY*(PTSTEP/XDAY) 01111 ! 01112 ! 01113 WHERE(OSFCMELT) 01114 ! 01115 ! 3. wet snow albedo rate of change: 01116 ! ---------------------------------- 01117 ! 01118 ZALBWET(:) = (PALBEDOSC(:) - ZANSMIN(:))*EXP(-XANS_T*PTSTEP/XDAY) & 01119 + ZANSMIN(:) 01120 ! 01121 ! Calculate albedo rate based on liquid water 01122 ! content of the first snow layer: During melt, if liquid 01123 ! content is at the maximum, then snow albedo 01124 ! decreases the fastest: 01125 ! 01126 ZWHOLDMAX(:) = SNOW3LHOLD(PSNOWRHO,PSNOWDZ) 01127 ! 01128 ZFRACLIQ(:) = MIN(1.0, PSNOWLIQ(:)/MAX(XSNOWDMIN,ZWHOLDMAX)) 01129 ! 01130 END WHERE 01131 ! 01132 ! 4. Updated Albedo: 01133 ! ------------------ 01134 ! 01135 PALBEDOSC(:) = ZFRACLIQ(:)*ZALBWET(:) + (1.-ZFRACLIQ(:))*ZALBDRY(:) & 01136 + ZPCPALB(:) 01137 ! 01138 PALBEDOSC(:) = MIN(ZANSMAX(:), MAX(ZANSMIN(:),PALBEDOSC(:))) 01139 IF (LHOOK) CALL DR_HOOK('SNOW3LALB',1,ZHOOK_HANDLE) 01140 ! 01141 !------------------------------------------------------------------------------- 01142 ! 01143 END SUBROUTINE SNOW3LALB 01144 !#################################################################### 01145 !#################################################################### 01146 !#################################################################### 01147 SUBROUTINE SNOW3LRAD(PSNOWDZMIN, PSW_RAD, PSNOWALB, PSNOWDZ, & 01148 PSNOWRHO, PALB, PRADSINK, PRADXS ) 01149 ! 01150 !! PURPOSE 01151 !! ------- 01152 ! Calculate the transmission of shortwave (solar) radiation 01153 ! through the snowpack (using a form of Beer's Law: exponential 01154 ! decay of radiation with increasing snow depth). 01155 ! 01156 IMPLICIT NONE 01157 ! 01158 !* 0.1 declarations of arguments 01159 ! 01160 REAL, INTENT(IN) :: PSNOWDZMIN 01161 ! 01162 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD, PSNOWALB, PALB 01163 ! 01164 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWDZ 01165 ! 01166 REAL, DIMENSION(:), INTENT(OUT) :: PRADXS 01167 ! 01168 REAL, DIMENSION(:,:), INTENT(OUT) :: PRADSINK 01169 ! 01170 ! 01171 !* 0.2 declarations of local variables 01172 ! 01173 INTEGER :: JJ, JI 01174 ! 01175 INTEGER :: INI 01176 INTEGER :: INLVLS 01177 ! 01178 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZRADTOT 01179 ! 01180 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOW, 01181 ZDSGRAIN, ZSXNU, 01182 ZCOEF, ZSNOWDZ 01183 ! 01184 ! ISBA-ES Radiation extinction coefficients: (see Loth and Graf 1993): 01185 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 01186 ! 01187 REAL, PARAMETER :: ZSNOWRAD_CVEXT = 3.8e-3 ! [(m5/2)/kg] 01188 REAL, PARAMETER :: ZSNOWRAD_AGRAIN = 1.6e-4 ! (m) 01189 REAL, PARAMETER :: ZSNOWRAD_BGRAIN = 1.1e-13 ! (m13/kg4) 01190 REAL, PARAMETER :: ZDSGRAIN_MAX = 2.796e-3 ! m 01191 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01192 !------------------------------------------------------------------------------- 01193 ! 01194 ! 0. Initialization: 01195 ! ------------------ 01196 ! 01197 IF (LHOOK) CALL DR_HOOK('SNOW3LRAD',0,ZHOOK_HANDLE) 01198 INI = SIZE(PSNOWDZ(:,:),1) 01199 INLVLS = SIZE(PSNOWDZ(:,:),2) 01200 ! 01201 ! 01202 ! 1. Vanishingly thin snowpack check: 01203 ! ----------------------------------- 01204 ! For vanishingly thin snowpacks, much of the radiation 01205 ! can pass through snowpack into underlying soil, making 01206 ! a large (albeit temporary) thermal gradient: by imposing 01207 ! a minimum thickness, this increases the radiation absorbtion 01208 ! for vanishingly thin snowpacks. 01209 ! 01210 ZSNOWDZ(:,:) = MAX(PSNOWDZMIN, PSNOWDZ(:,:)) 01211 ! 01212 ! 01213 ! 2. Extinction of net shortwave radiation 01214 ! ---------------------------------------- 01215 ! Fn of snow depth and density (Loth and Graf 1993: 01216 ! SNOWCVEXT => from Bohren and Barkstrom 1974 01217 ! SNOWAGRAIN and SNOWBGRAIN=> from Jordan 1976) 01218 ! 01219 ! Calculate grid levels: 01220 ! 01221 ZSNOW(:,1) = ZSNOWDZ(:,1) 01222 DO JJ=2,INLVLS 01223 DO JI=1,INI 01224 ZSNOW(JI,JJ) = ZSNOW(JI,JJ-1) + ZSNOWDZ(JI,JJ) 01225 ENDDO 01226 ENDDO 01227 ! 01228 ! Snow grain size: 01229 ! 01230 ZDSGRAIN(:,:) = MIN(ZDSGRAIN_MAX, ZSNOWRAD_AGRAIN + ZSNOWRAD_BGRAIN*(PSNOWRHO(:,:)**4)) 01231 ! 01232 ! NOTE: as this is a "sink" term for heat, it 01233 ! as assigned a negative value. 01234 ! 01235 ZSXNU(:,:) = ZSNOWRAD_CVEXT*PSNOWRHO(:,:)/SQRT(ZDSGRAIN(:,:)) 01236 ZCOEF(:,:) = EXP(-ZSXNU(:,:)*ZSNOW(:,:)) 01237 ! 01238 ! 3. Radiation at each level: (W/m2) 01239 ! ---------------------------------- 01240 ! 01241 DO JJ=1,INLVLS 01242 DO JI=1,INI 01243 PRADSINK(JI,JJ) = -PSW_RAD(JI)*(1.-PSNOWALB(JI))*ZCOEF(JI,JJ) 01244 ENDDO 01245 ENDDO 01246 ! 01247 ! For thin snow packs, radiation might reach base of 01248 ! snowpack...so we influence this amount with sfc albedo 01249 ! and (outside of this routine) add any excess heat 01250 ! to underlying soil: 01251 ! 01252 PRADSINK(:,INLVLS) = PRADSINK(:,INLVLS)*(1.0-PALB(:)) 01253 ! 01254 ! 4. Excess radiation not absorbed by snowpack (W/m2): 01255 ! ---------------------------------------------------- 01256 ! 01257 ZRADTOT(:) = PRADSINK(:,1)+(1.-PSNOWALB(:))*PSW_RAD(:) 01258 DO JJ=2,INLVLS 01259 DO JI=1,INI 01260 ZRADTOT(JI) = ZRADTOT(JI) + PRADSINK(JI,JJ)-PRADSINK(JI,JJ-1) 01261 ENDDO 01262 ENDDO 01263 ! 01264 PRADXS(:) = (1.-PSNOWALB(:))*PSW_RAD(:) - ZRADTOT(:) 01265 IF (LHOOK) CALL DR_HOOK('SNOW3LRAD',1,ZHOOK_HANDLE) 01266 ! 01267 !------------------------------------------------------------------------------- 01268 ! 01269 END SUBROUTINE SNOW3LRAD 01270 !#################################################################### 01271 !#################################################################### 01272 !#################################################################### 01273 SUBROUTINE SNOW3LTHRM(PSNOWRHO,PSCOND,PSNOWTEMP,PPS) 01274 ! 01275 !! PURPOSE 01276 !! ------- 01277 ! Calculate snow thermal conductivity from 01278 ! Sun et al. 1999, J. of Geophys. Res., 104, 19587-19579 01279 ! (vapor) and Anderson, 1976, NOAA Tech. Rep. NWS 19 (snow). 01280 ! 01281 ! 01282 USE MODD_CSTS,ONLY : XP00 01283 ! 01284 IMPLICIT NONE 01285 ! 01286 !* 0.1 declarations of arguments 01287 ! 01288 REAL, DIMENSION(:), INTENT(IN) :: PPS 01289 ! 01290 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP, PSNOWRHO 01291 ! 01292 REAL, DIMENSION(:,:), INTENT(OUT) :: PSCOND 01293 ! 01294 ! 01295 !* 0.2 declarations of local variables 01296 ! 01297 INTEGER :: JJ, JI 01298 ! 01299 INTEGER :: INI 01300 INTEGER :: INLVLS 01301 ! 01302 ! ISBA-ES Thermal conductivity coefficients from Anderson (1976): 01303 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 01304 ! 01305 REAL, PARAMETER :: ZSNOWTHRMCOND1 = 0.02 ! [W/(m K)] 01306 REAL, PARAMETER :: ZSNOWTHRMCOND2 = 2.5E-6 ! [W m5/(kg2 K)] 01307 ! 01308 ! ISBA-ES Thermal conductivity: Implicit vapor diffn effects 01309 ! (sig only for new snow OR high altitudes) 01310 ! from Sun et al. (1999): based on data from Jordan (1991) 01311 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 01312 ! 01313 REAL, PARAMETER :: ZSNOWTHRMCOND_AVAP = -0.06023 ! [W/(m K)] 01314 REAL, PARAMETER :: ZSNOWTHRMCOND_BVAP = -2.5425 ! (W/m) 01315 REAL, PARAMETER :: ZSNOWTHRMCOND_CVAP = -289.99 ! (K) 01316 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01317 !------------------------------------------------------------------------------- 01318 ! 01319 INI = SIZE(PSNOWRHO(:,:),1) 01320 INLVLS = SIZE(PSNOWRHO(:,:),2) 01321 ! 01322 ! 1. Snow thermal conductivity 01323 ! ---------------------------- 01324 ! 01325 IF (LHOOK) CALL DR_HOOK('SNOW3LTHRM',0,ZHOOK_HANDLE) 01326 DO JJ=1,INLVLS 01327 DO JI=1,INI 01328 ! 01329 PSCOND(JI,JJ) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*PSNOWRHO(JI,JJ)*PSNOWRHO(JI,JJ)) & 01330 + MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(PSNOWTEMP(JI,JJ) & 01331 + ZSNOWTHRMCOND_CVAP)))*(XP00/PPS(JI))) 01332 ! 01333 ENDDO 01334 ENDDO 01335 IF (LHOOK) CALL DR_HOOK('SNOW3LTHRM',1,ZHOOK_HANDLE) 01336 ! 01337 !------------------------------------------------------------------------------- 01338 ! 01339 END SUBROUTINE SNOW3LTHRM 01340 !#################################################################### 01341 !#################################################################### 01342 !#################################################################### 01343 SUBROUTINE SNOW3LEBUD(HSNOWRES, HIMPLICIT_WIND, & 01344 PPEW_A_COEF, PPEW_B_COEF, & 01345 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, & 01346 PSNOWDZMIN, & 01347 PZREF,PTS,PSNOWRHO,PSNOWLIQ,PSCAP,PSCOND1,PSCOND2, & 01348 PUREF,PEXNS,PEXNA,PDIRCOSZW,PVMOD, & 01349 PLW_RAD,PSW_RAD,PTA,PQA,PPS,PTSTEP, & 01350 PSNOWDZ1,PSNOWDZ2,PALBT,PZ0,PZ0EFF,PZ0H, & 01351 PSFCFRZ,PRADSINK,PHPSNOW, & 01352 PCT,PEMIST,PRHOA,PTSTERM1,PTSTERM2,PRA,PCDSNOW,PCHSNOW, & 01353 PQSAT,PDQSAT,PRSRA,PUSTAR2_IC,PRI, & 01354 PPET_A_COEF_T,PPEQ_A_COEF_T,PPET_B_COEF_T,PPEQ_B_COEF_T ) 01355 ! 01356 !! PURPOSE 01357 !! ------- 01358 ! Calculate surface energy budget linearization (terms) and turbulent 01359 ! exchange coefficients/resistance between surface and atmosphere. 01360 ! (Noilhan and Planton 1989; Giordani 1993; Noilhan and Mahfouf 1996) 01361 ! 01362 USE MODD_SURF_PAR, ONLY : XUNDEF 01363 USE MODD_CSTS, ONLY : XCPD, XRHOLW, XSTEFAN, XLVTT, XLSTT 01364 USE MODD_SNOW_PAR, ONLY : X_RI_MAX, XEMISSN 01365 ! 01366 USE MODE_THERMOS 01367 ! 01368 USE MODI_SURFACE_RI 01369 USE MODI_SURFACE_AERO_COND 01370 USE MODI_SURFACE_CD 01371 ! 01372 ! 01373 IMPLICIT NONE 01374 ! 01375 !* 0.1 declarations of arguments 01376 ! 01377 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN 01378 ! 01379 CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES ! type of sfc resistance 01380 ! DEFAULT=Louis (1979), standard ISBA 01381 ! method. Option to limit Ri number 01382 ! for very srtable conditions 01383 ! 01384 CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option 01385 ! ! 'OLD' = direct 01386 ! ! 'NEW' = Taylor serie, order 1 01387 ! 01388 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, 01389 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, 01390 PPEQ_B_COEF 01391 ! PPEW_A_COEF = wind coefficient (m2s/kg) 01392 ! PPEW_B_COEF = wind coefficient (m/s) 01393 ! PPET_A_COEF = A-air temperature coefficient 01394 ! PPET_B_COEF = B-air temperature coefficient 01395 ! PPEQ_A_COEF = A-air specific humidity coefficient 01396 ! PPEQ_B_COEF = B-air specific humidity coefficient 01397 ! 01398 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PTS, PSNOWDZ1, PSNOWDZ2, 01399 PRADSINK, PSNOWRHO, PSNOWLIQ, PSCAP, 01400 PSCOND1, PSCOND2, 01401 PZ0, PHPSNOW, 01402 PALBT, PZ0EFF, PZ0H 01403 ! 01404 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD, PLW_RAD, PTA, PQA, PPS, PRHOA 01405 ! 01406 REAL, DIMENSION(:), INTENT(IN) :: PUREF, PEXNS, PEXNA, PDIRCOSZW, PVMOD 01407 ! 01408 REAL, DIMENSION(:), INTENT(OUT) :: PTSTERM1, PTSTERM2, PEMIST, PRA, 01409 PCT, PSFCFRZ, PCDSNOW, PCHSNOW, 01410 PQSAT, PDQSAT, PRSRA 01411 ! 01412 REAL, DIMENSION(:), INTENT(OUT) :: PUSTAR2_IC, 01413 PPET_A_COEF_T, PPEQ_A_COEF_T, 01414 PPET_B_COEF_T, PPEQ_B_COEF_T 01415 ! 01416 REAL, DIMENSION(:), INTENT(OUT) :: PRI 01417 ! 01418 !* 0.2 declarations of local variables 01419 ! 01420 REAL, DIMENSION(SIZE(PTS)) :: ZAC, ZRI, 01421 ZSCONDA, ZA, ZB, ZC, 01422 ZCDN, ZSNOWDZM1, ZSNOWDZM2, 01423 ZVMOD, ZUSTAR2, ZTS3, ZLVT, 01424 Z_CCOEF 01425 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01426 ! 01427 !------------------------------------------------------------------------------- 01428 ! 01429 ! 1. New saturated specific humidity and derrivative: 01430 ! --------------------------------------------------- 01431 ! 01432 IF (LHOOK) CALL DR_HOOK('SNOW3LEBUD',0,ZHOOK_HANDLE) 01433 ! 01434 ZRI (:) = XUNDEF 01435 ! 01436 PQSAT (:) = QSATI(PTS(:),PPS(:)) 01437 PDQSAT(:) = DQSATI(PTS(:),PPS(:),PQSAT(:)) 01438 ! 01439 ! 01440 ! 2. Surface properties: 01441 ! ---------------------- 01442 ! It might be of interest to use snow-specific roughness 01443 ! or a temperature dependence on emissivity: 01444 ! but for now, use ISBA defaults. 01445 ! 01446 PEMIST(:) = XEMISSN 01447 ! 01448 ! 2. Computation of resistance and drag coefficient 01449 ! ------------------------------------------------- 01450 ! 01451 CALL SURFACE_RI(PTS, PQSAT, PEXNS, PEXNA, PTA, PQA, & 01452 PZREF, PUREF, PDIRCOSZW, PVMOD, ZRI ) 01453 ! 01454 ! Simple adaptation of method by Martin and Lejeune (1998) 01455 ! to apply a lower limit to turbulent transfer coef 01456 ! by defining a maximum Richarson number for stable 01457 ! conditions: 01458 ! 01459 IF(HSNOWRES=='RIL') ZRI(:) = MIN(X_RI_MAX, ZRI(:)) 01460 ! 01461 PRI(:)=ZRI(:) 01462 ! 01463 ! Surface aerodynamic resistance for heat transfers 01464 ! 01465 CALL SURFACE_AERO_COND(ZRI, PZREF, PUREF, PVMOD, PZ0, PZ0H, ZAC, PRA, PCHSNOW) 01466 ! 01467 ! For atmospheric model coupling: 01468 ! 01469 CALL SURFACE_CD(ZRI, PZREF, PUREF, PZ0EFF, PZ0H, PCDSNOW, ZCDN) 01470 ! 01471 PRSRA(:) = PRHOA(:) / PRA(:) 01472 ! 01473 ! 01474 ! Modify flux-form implicit coupling coefficients: 01475 ! - wind components: 01476 ! 01477 IF(HIMPLICIT_WIND=='OLD')THEN 01478 ! old implicitation 01479 ZUSTAR2(:) = ( PCDSNOW(:)*PVMOD(:)*PPEW_B_COEF(:)) / & 01480 (1.0-PRHOA(:)*PCDSNOW(:)*PVMOD(:)*PPEW_A_COEF(:)) 01481 ELSE 01482 ! new implicitation 01483 ZUSTAR2(:) = (PCDSNOW(:)*PVMOD(:)*(2.*PPEW_B_COEF(:)-PVMOD(:))) & 01484 / (1.0-2.0*PRHOA(:)*PCDSNOW(:)*PVMOD(:)*PPEW_A_COEF(:)) 01485 ENDIF 01486 ! 01487 ZVMOD(:) = PRHOA(:)*PPEW_A_COEF(:)*ZUSTAR2(:) + PPEW_B_COEF(:) 01488 ZVMOD(:) = MAX(ZVMOD(:),0.) 01489 ! 01490 WHERE(PPEW_A_COEF(:)/= 0.) 01491 ZUSTAR2(:) = MAX( ( ZVMOD(:) - PPEW_B_COEF(:) ) / (PRHOA(:)*PPEW_A_COEF(:)), 0.) 01492 ENDWHERE 01493 ! 01494 ! implicit wind friction 01495 ZUSTAR2(:) = MAX(ZUSTAR2(:),0.) 01496 ! 01497 PUSTAR2_IC(:) = ZUSTAR2(:) 01498 ! 01499 ! 3. Calculate linearized surface energy budget components: 01500 ! --------------------------------------------------------- 01501 ! To prevent numerical difficulties for very thin snow 01502 ! layers, limit the grid "thinness": this is important as 01503 ! layers become vanishing thin: 01504 ! 01505 ZSNOWDZM1(:) = MAX(PSNOWDZ1(:), PSNOWDZMIN) 01506 ZSNOWDZM2(:) = MAX(PSNOWDZ2(:), PSNOWDZMIN) 01507 ! 01508 ! Surface thermal inertia: 01509 ! 01510 PCT(:) = 1.0/(PSCAP(:)*ZSNOWDZM1(:)) 01511 ! 01512 ! Fraction of surface frozen (sublimation) with the remaining 01513 ! fraction being liquid (evaporation): 01514 ! 01515 PSFCFRZ(:) = 1.0 - PSNOWLIQ(:)*XRHOLW/(ZSNOWDZM1(:)*PSNOWRHO(:)) 01516 ! 01517 ! Thermal conductivity between uppermost and lower snow layers: 01518 ! 01519 ZSCONDA(:) = (ZSNOWDZM1(:)*PSCOND1(:) + ZSNOWDZM2(:)*PSCOND2(:))/ & 01520 (ZSNOWDZM1(:)+ZSNOWDZM2(:)) 01521 ! 01522 ! Transform implicit coupling coefficients: 01523 ! Note, surface humidity is 100% over snow. 01524 ! 01525 ! - specific humidity: 01526 ! 01527 Z_CCOEF(:) = 1.0 - PPEQ_A_COEF(:)*PRSRA(:) 01528 ! 01529 PPEQ_A_COEF_T(:) = - PPEQ_A_COEF(:)*PRSRA(:)*PDQSAT(:)/Z_CCOEF(:) 01530 ! 01531 PPEQ_B_COEF_T(:) = ( PPEQ_B_COEF(:) - PPEQ_A_COEF(:)*PRSRA(:)*(PQSAT(:) - & 01532 PDQSAT(:)*PTS(:)) )/Z_CCOEF(:) 01533 ! 01534 ! - air temperature: 01535 ! (assumes A and B correspond to potential T): 01536 ! 01537 Z_CCOEF(:) = (1.0 - PPET_A_COEF(:)*PRSRA(:))/PEXNA(:) 01538 ! 01539 PPET_A_COEF_T(:) = - PPET_A_COEF(:)*PRSRA(:)/(PEXNS(:)*Z_CCOEF(:)) 01540 ! 01541 PPET_B_COEF_T(:) = PPET_B_COEF(:)/Z_CCOEF(:) 01542 ! 01543 ! 01544 ! Energy budget solution terms: 01545 ! 01546 ZTS3(:) = PEMIST(:) * XSTEFAN * PTS(:)**3 01547 ZLVT(:) = (1.-PSFCFRZ(:))*XLVTT + PSFCFRZ(:)*XLSTT 01548 ! 01549 ZA(:) = 1. / PTSTEP + PCT(:) * (4. * ZTS3(:) + & 01550 PRSRA(:) * ZLVT(:) * (PDQSAT(:) - PPEQ_A_COEF_T(:)) & 01551 + PRSRA(:) * XCPD * ( (1./PEXNS(:))-(PPET_A_COEF_T(:)/PEXNA(:)) ) & 01552 + (2*ZSCONDA(:)/(ZSNOWDZM2(:)+ZSNOWDZM1(:))) ) 01553 ! 01554 ZB(:) = 1. / PTSTEP + PCT(:) * (3. * ZTS3(:) + & 01555 PRSRA(:) * PDQSAT(:) * ZLVT(:) ) 01556 ! 01557 ZC(:) = PCT(:) * (PRSRA(:) * XCPD * PPET_B_COEF_T(:)/PEXNA(:) + PSW_RAD(:) * & 01558 (1. - PALBT(:)) + PEMIST(:)*PLW_RAD(:) - PRSRA(:) * & 01559 ZLVT(:) * (PQSAT(:)-PPEQ_B_COEF_T(:)) & 01560 + PHPSNOW(:) + PRADSINK(:) ) 01561 ! 01562 ! 01563 ! Coefficients needed for implicit solution 01564 ! of linearized surface energy budget: 01565 ! 01566 PTSTERM2(:) = 2*ZSCONDA(:)*PCT(:)/(ZA(:)*(ZSNOWDZM2(:)+ZSNOWDZM1(:))) 01567 ! 01568 PTSTERM1(:) = (PTS(:)*ZB(:) + ZC(:))/ZA(:) 01569 IF (LHOOK) CALL DR_HOOK('SNOW3LEBUD',1,ZHOOK_HANDLE) 01570 ! 01571 !------------------------------------------------------------------------------- 01572 ! 01573 END SUBROUTINE SNOW3LEBUD 01574 !#################################################################### 01575 !#################################################################### 01576 !#################################################################### 01577 SUBROUTINE SNOW3LSOLVT(PTSTEP,PSNOWDZMIN, & 01578 PSNOWDZ,PSCOND,PSCAP,PTG, & 01579 PSOILCOND,PD_G, & 01580 PRADSINK,PCT,PTERM1,PTERM2, & 01581 PPET_A_COEF_T,PPEQ_A_COEF_T, & 01582 PPET_B_COEF_T,PPEQ_B_COEF_T, & 01583 PTA_IC, PQA_IC, & 01584 PGBAS,PGBASO,PSNOWTEMP,PSNOWFLUX ) 01585 ! 01586 !! PURPOSE 01587 !! ------- 01588 ! This subroutine solves the 1-d diffusion of 'ZSNOWTEMP' using a 01589 ! layer averaged set of equations which are time differenced 01590 ! using the backward-difference scheme (implicit). 01591 ! The eqs are solved rapidly by taking advantage of the 01592 ! fact that the matrix is tridiagonal. This is a very 01593 ! general routine and can be used for the 1-d diffusion of any 01594 ! quantity as long as the diffusity is not a function of the 01595 ! quantity being diffused. Aaron Boone 8-98. Soln to the eq: 01596 ! 01597 ! c dQ d k dQ dS 01598 ! -- = -- -- - -- 01599 ! dt dx dx dx 01600 ! 01601 ! where k = k(x) (thermal conductivity), c = c(x) (heat capacity) 01602 ! as an eg. for temperature/heat eq. S is a sink (-source) term. 01603 ! Diffusivity is k/c 01604 ! 01605 ! 01606 USE MODD_CSTS,ONLY : XTT 01607 ! 01608 USE MODI_TRIDIAG_GROUND 01609 ! 01610 IMPLICIT NONE 01611 ! 01612 !* 0.1 declarations of arguments 01613 ! 01614 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN 01615 ! 01616 REAL, DIMENSION(:), INTENT(IN) :: PTG, PSOILCOND, PD_G, 01617 PCT, PTERM1, PTERM2 01618 01619 ! 01620 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ, PSCOND, PSCAP, 01621 PRADSINK 01622 ! 01623 REAL, DIMENSION(:), INTENT(IN) :: PPET_A_COEF_T, PPEQ_A_COEF_T, 01624 PPET_B_COEF_T, PPEQ_B_COEF_T 01625 ! 01626 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP 01627 ! 01628 REAL, DIMENSION(:), INTENT(OUT) :: PGBAS, PGBASO, PSNOWFLUX, 01629 PTA_IC, PQA_IC 01630 ! 01631 ! 01632 !* 0.2 declarations of local variables 01633 ! 01634 ! 01635 INTEGER :: JJ, JI 01636 ! 01637 INTEGER :: INI 01638 INTEGER :: INLVLS 01639 ! 01640 REAL, DIMENSION(SIZE(PTG)) :: ZSNOWTEMP_DELTA 01641 ! 01642 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSNOWTEMP, ZDTERM, ZCTERM, 01643 ZFRCV, ZAMTRX, ZBMTRX, 01644 ZCMTRX 01645 ! 01646 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZWORK1, ZWORK2, ZDZDIF, 01647 ZSNOWDZM 01648 ! 01649 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)-1) :: ZSNOWTEMP_M, 01650 ZFRCV_M, ZAMTRX_M, 01651 ZBMTRX_M, ZCMTRX_M 01652 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01653 !------------------------------------------------------------------------------- 01654 ! 01655 ! 0. Initialize: 01656 ! ------------------ 01657 ! 01658 IF (LHOOK) CALL DR_HOOK('SNOW3LSOLVT',0,ZHOOK_HANDLE) 01659 ZSNOWTEMP(:,:) = PSNOWTEMP(:,:) 01660 INI = SIZE(PSNOWDZ(:,:),1) 01661 INLVLS = SIZE(PSNOWDZ(:,:),2) 01662 ! 01663 ! 01664 ! 1. Calculate tri-diagnoal matrix coefficients: 01665 ! ---------------------------------------------- 01666 ! For heat transfer, assume a minimum grid 01667 ! thickness (to prevent numerical 01668 ! problems for very thin snow cover): 01669 ! 01670 ZSNOWDZM(:,:) = MAX(PSNOWDZ(:,:), PSNOWDZMIN) 01671 ! 01672 ZWORK1(:,:) = ZSNOWDZM(:,:)*PSCOND(:,:) 01673 DO JJ=1,INLVLS-1 01674 DO JI=1,INI 01675 ZDZDIF(JI,JJ) = ZSNOWDZM(JI,JJ) + ZSNOWDZM(JI,JJ+1) 01676 ZWORK2(JI,JJ) = ZSNOWDZM(JI,JJ+1)*PSCOND(JI,JJ+1) 01677 ENDDO 01678 ENDDO 01679 ZDZDIF(:,INLVLS) = ZSNOWDZM(:,INLVLS) + PD_G(:) 01680 ZWORK2(:,INLVLS) = PD_G(:)*PSOILCOND(:) 01681 ! 01682 ZDTERM(:,:) = 2.0*(ZWORK1(:,:)+ZWORK2(:,:))/(ZDZDIF(:,:)*ZDZDIF(:,:)) 01683 ! 01684 ZCTERM(:,:) = PSCAP(:,:)*ZSNOWDZM(:,:)/PTSTEP 01685 ! 01686 ! 2. Set up tri-diagonal matrix 01687 ! ----------------------------- 01688 ! 01689 ! Upper BC 01690 ! 01691 ZAMTRX(:,1) = 0.0 01692 ZBMTRX(:,1) = 1./(PCT(:)*PTSTEP) 01693 ZCMTRX(:,1) = -PTERM2(:)*ZBMTRX(:,1) 01694 ZFRCV(:,1) = PTERM1(:)*ZBMTRX(:,1) 01695 ! 01696 ! 01697 ! Interior Grid 01698 ! 01699 DO JJ=2,INLVLS-1 01700 DO JI=1,INI 01701 ZAMTRX(JI,JJ) = -ZDTERM(JI,JJ-1) 01702 ZBMTRX(JI,JJ) = ZCTERM(JI,JJ) + ZDTERM(JI,JJ-1) + ZDTERM(JI,JJ) 01703 ZCMTRX(JI,JJ) = -ZDTERM(JI,JJ) 01704 ZFRCV (JI,JJ) = ZCTERM(JI,JJ)*PSNOWTEMP(JI,JJ) - (PRADSINK(JI,JJ-1)-PRADSINK(JI,JJ)) 01705 ENDDO 01706 ENDDO 01707 ! 01708 !Lower BC 01709 ! 01710 ZAMTRX(:,INLVLS) = -ZDTERM(:,INLVLS-1) 01711 ZBMTRX(:,INLVLS) = ZCTERM(:,INLVLS) + ZDTERM(:,INLVLS-1) + & 01712 ZDTERM(:,INLVLS) 01713 ZCMTRX(:,INLVLS) = 0.0 01714 ZFRCV(:,INLVLS) = ZCTERM(:,INLVLS)*PSNOWTEMP(:,INLVLS) + & 01715 ZDTERM(:,INLVLS)*PTG(:) & 01716 - (PRADSINK(:,INLVLS-1)-PRADSINK(:,INLVLS)) 01717 ! 01718 ! - - ------------------------------------------------- 01719 ! 01720 ! 4. Compute solution vector 01721 ! -------------------------- 01722 ! 01723 CALL TRIDIAG_GROUND(ZAMTRX,ZBMTRX,ZCMTRX,ZFRCV,ZSNOWTEMP) 01724 ! 01725 ! Heat flux between surface and 2nd snow layers: (W/m2) 01726 ! 01727 PSNOWFLUX(:) = ZDTERM(:,1)*(ZSNOWTEMP(:,1) - ZSNOWTEMP(:,2)) 01728 ! 01729 ! 01730 ! 5. Snow melt case 01731 ! ----------------- 01732 ! If melting in uppermost layer, assume surface layer 01733 ! temperature at freezing point and re-evaluate lower 01734 ! snowpack temperatures. This is done as it is most likely 01735 ! most signigant melting will occur within a time step in surface layer. 01736 ! Surface energy budget (and fluxes) will 01737 ! be re-calculated (outside of this routine): 01738 ! 01739 ZAMTRX_M(:,1) = 0.0 01740 ZBMTRX_M(:,1) = ZCTERM(:,2) + ZDTERM(:,1) + ZDTERM(:,2) 01741 ZCMTRX_M(:,1) = -ZDTERM(:,2) 01742 ZFRCV_M(:,1) = ZCTERM(:,2)*PSNOWTEMP(:,2) + XTT*ZDTERM(:,1) - & 01743 (PRADSINK(:,1)-PRADSINK(:,2)) 01744 ! 01745 DO JJ=2,INLVLS-1 01746 DO JI=1,INI 01747 ZAMTRX_M (JI,JJ) = ZAMTRX (JI,JJ+1) 01748 ZBMTRX_M (JI,JJ) = ZBMTRX (JI,JJ+1) 01749 ZCMTRX_M (JI,JJ) = ZCMTRX (JI,JJ+1) 01750 ZFRCV_M (JI,JJ) = ZFRCV (JI,JJ+1) 01751 ZSNOWTEMP_M(JI,JJ) = PSNOWTEMP(JI,JJ+1) 01752 ENDDO 01753 ENDDO 01754 ! 01755 CALL TRIDIAG_GROUND(ZAMTRX_M,ZBMTRX_M,ZCMTRX_M,ZFRCV_M,ZSNOWTEMP_M) 01756 ! 01757 ! If melting for 2 consecuative time steps, then replace current T-profile 01758 ! with one assuming T=Tf in surface layer: 01759 ! 01760 ZSNOWTEMP_DELTA(:) = 0.0 01761 ! 01762 WHERE(ZSNOWTEMP(:,1) > XTT .AND. PSNOWTEMP(:,1) == XTT) 01763 PSNOWFLUX(:) = ZDTERM(:,1)*(XTT - ZSNOWTEMP_M(:,1)) 01764 ZSNOWTEMP_DELTA(:) = 1.0 01765 END WHERE 01766 ! 01767 DO JJ=2,INLVLS 01768 DO JI=1,INI 01769 ZSNOWTEMP(JI,JJ) = ZSNOWTEMP_DELTA(JI)*ZSNOWTEMP_M(JI,JJ-1) & 01770 + (1.0-ZSNOWTEMP_DELTA(JI))*ZSNOWTEMP(JI,JJ) 01771 ENDDO 01772 ENDDO 01773 ! 01774 ! 01775 ! 6. Lower boundary flux: 01776 ! ----------------------- 01777 ! NOTE: evaluate this term assuming the snow layer 01778 ! can't exceed the freezing point as this adjustment 01779 ! is made in melting routine. Then must adjust temperature 01780 ! to conserve energy: 01781 ! 01782 PGBASO(:) = ZDTERM(:,INLVLS)*(ZSNOWTEMP(:,INLVLS) -PTG(:)) 01783 PGBAS(:) = ZDTERM(:,INLVLS)*(MIN(XTT,ZSNOWTEMP(:,INLVLS))-PTG(:)) 01784 ! 01785 ZSNOWTEMP(:,INLVLS) = ZSNOWTEMP(:,INLVLS) + (PGBASO(:)-PGBAS(:))/ZCTERM(:,INLVLS) 01786 ! 01787 ! 7. Update temperatute profile in time: 01788 ! -------------------------------------- 01789 ! 01790 PSNOWTEMP(:,:) = ZSNOWTEMP(:,:) 01791 ! 01792 ! 01793 ! 8. Compute new (implicit) air T and specific humidity 01794 ! ----------------------------------------------------- 01795 ! 01796 PTA_IC(:) = PPET_B_COEF_T(:) + PPET_A_COEF_T(:)* PSNOWTEMP(:,1) 01797 ! 01798 PQA_IC(:) = PPEQ_B_COEF_T(:) + PPEQ_A_COEF_T(:)* PSNOWTEMP(:,1) 01799 IF (LHOOK) CALL DR_HOOK('SNOW3LSOLVT',1,ZHOOK_HANDLE) 01800 ! 01801 ! 01802 END SUBROUTINE SNOW3LSOLVT 01803 !#################################################################### 01804 !#################################################################### 01805 !#################################################################### 01806 SUBROUTINE SNOW3LMELT(PTSTEP,PSCAP,PSNOWTEMP,PSNOWDZ, & 01807 PSNOWRHO,PSNOWLIQ,PMELTXS ) 01808 ! 01809 ! 01810 !! PURPOSE 01811 !! ------- 01812 ! Calculate snow melt (resulting from surface fluxes, ground fluxes, 01813 ! or internal shortwave radiation absorbtion). It is used to 01814 ! augment liquid water content, maintain temperatures 01815 ! at or below freezing, and possibly reduce the mass 01816 ! or compact the layer(s). 01817 ! 01818 ! 01819 USE MODD_CSTS,ONLY : XTT, XLMTT, XRHOLW, XRHOLI 01820 ! 01821 USE MODE_SNOW3L 01822 ! 01823 IMPLICIT NONE 01824 ! 01825 !* 0.1 declarations of arguments 01826 ! 01827 REAL, INTENT(IN) :: PTSTEP 01828 ! 01829 REAL, DIMENSION(:,:), INTENT(IN) :: PSCAP 01830 ! 01831 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWRHO, 01832 PSNOWLIQ 01833 ! 01834 REAL, DIMENSION(:), INTENT(OUT) :: PMELTXS 01835 ! 01836 ! 01837 !* 0.2 declarations of local variables 01838 ! 01839 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZPHASE, ZCMPRSFACT, 01840 ZSNOWLWE, ZWHOLDMAX, 01841 ZSNOWMELT, ZSNOWTEMP, 01842 ZMELTXS 01843 ! 01844 INTEGER :: JWRK, JI ! loop counter 01845 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01846 !------------------------------------------------------------------------------- 01847 ! 01848 ! 0. Initialize: 01849 ! --------------------------- 01850 ! 01851 IF (LHOOK) CALL DR_HOOK('SNOW3LMELT',0,ZHOOK_HANDLE) 01852 ZPHASE(:,:) = 0.0 01853 ZCMPRSFACT(:,:) = 0.0 01854 ZSNOWLWE(:,:) = 0.0 01855 ZWHOLDMAX(:,:) = 0.0 01856 ZSNOWMELT(:,:) = 0.0 01857 ZSNOWTEMP(:,:) = 0.0 01858 ZMELTXS(:,:) = 0.0 01859 ! 01860 ! 1. Determine amount of melt in each layer: 01861 ! ------------------------------------------ 01862 ! 01863 WHERE(PSNOWDZ > 0.0) 01864 ! 01865 ! Total Liquid equivalent water content of snow (m): 01866 ! 01867 ZSNOWLWE(:,:) = PSNOWRHO(:,:)*PSNOWDZ(:,:)/XRHOLW 01868 ! 01869 ! Melt snow if excess energy and snow available: 01870 ! Phase change (J/m2) 01871 ! 01872 ZPHASE(:,:) = MIN(PSCAP(:,:)*MAX(0.0, PSNOWTEMP(:,:) - XTT)* & 01873 PSNOWDZ(:,:), & 01874 MAX(0.0,ZSNOWLWE(:,:)-PSNOWLIQ(:,:))*XLMTT*XRHOLW) 01875 ! 01876 ! 01877 ! Update snow liq water content and temperature if melting: 01878 ! liquid water available for next layer from melting of snow 01879 ! which is assumed to be leaving the current layer (m): 01880 ! 01881 ZSNOWMELT(:,:) = ZPHASE(:,:)/(XLMTT*XRHOLW) 01882 ! 01883 ! Cool off snow layer temperature due to melt: 01884 ! 01885 ZSNOWTEMP(:,:) = PSNOWTEMP(:,:) - ZPHASE(:,:)/(PSCAP(:,:)*PSNOWDZ(:,:)) 01886 ! 01887 PSNOWTEMP(:,:) = MIN(XTT, ZSNOWTEMP(:,:)) 01888 ! 01889 ZMELTXS(:,:) = (ZSNOWTEMP(:,:)-PSNOWTEMP(:,:))*PSCAP(:,:)*PSNOWDZ(:,:) 01890 ! 01891 ! Loss of snowpack depth: (m) and liquid equiv (m): 01892 ! Compression factor for melt loss: this decreases 01893 ! layer thickness and increases density thereby leaving 01894 ! total SWE constant. NOTE: All melt water 01895 ! in excess of the holding capacity is NOT used 01896 ! for compression, rather it decreases the layer 01897 ! thickness ONLY, causing a loss of SWE (outside 01898 ! of this routine). 01899 ! 01900 ZWHOLDMAX(:,:) = SNOW3LHOLD(PSNOWRHO,PSNOWDZ) 01901 ! 01902 ZCMPRSFACT(:,:) = (ZSNOWLWE(:,:)-MIN(PSNOWLIQ(:,:)+ZSNOWMELT(:,:), & 01903 ZWHOLDMAX(:,:)))/ & 01904 (ZSNOWLWE(:,:)-MIN(PSNOWLIQ(:,:),ZWHOLDMAX(:,:))) 01905 ! 01906 PSNOWDZ(:,:) = PSNOWDZ(:,:)*ZCMPRSFACT(:,:) 01907 PSNOWRHO(:,:) = ZSNOWLWE(:,:)*XRHOLW/PSNOWDZ(:,:) 01908 ! 01909 ! Make sure maximum density is not surpassed! If it is, lower the density 01910 ! and increase the snow thickness accordingly: 01911 01912 ZCMPRSFACT(:,:) = MAX(XRHOLI, PSNOWRHO(:,:))/XRHOLI 01913 PSNOWDZ(:,:) = PSNOWDZ(:,:)*ZCMPRSFACT(:,:) 01914 PSNOWRHO(:,:) = ZSNOWLWE(:,:)*XRHOLW/PSNOWDZ(:,:) 01915 ! 01916 ! 01917 ! 2. Add snow melt to current snow liquid water content: 01918 ! ------------------------------------------------------ 01919 ! 01920 PSNOWLIQ(:,:) = PSNOWLIQ(:,:) + ZSNOWMELT(:,:) 01921 ! 01922 END WHERE 01923 ! 01924 ! 3. Excess heat from melting 01925 ! --------------------------- 01926 ! use it to warm underlying ground/vegetation layer to conserve energy 01927 ! 01928 PMELTXS(:) = 0. 01929 DO JWRK = 1, SIZE(ZMELTXS,2) 01930 DO JI = 1, SIZE(ZMELTXS,1) 01931 PMELTXS(JI) = PMELTXS(JI) + ZMELTXS(JI,JWRK) 01932 ENDDO 01933 ENDDO 01934 PMELTXS(:) = PMELTXS(:) / PTSTEP ! (W/m2) 01935 IF (LHOOK) CALL DR_HOOK('SNOW3LMELT',1,ZHOOK_HANDLE) 01936 ! 01937 ! 01938 ! 01939 END SUBROUTINE SNOW3LMELT 01940 !#################################################################### 01941 !#################################################################### 01942 !#################################################################### 01943 SUBROUTINE SNOW3LREFRZ(PTSTEP,PRR, & 01944 PSNOWRHO,PSNOWTEMP,PSNOWDZ,PSNOWLIQ, & 01945 PTHRUFAL ) 01946 ! 01947 ! 01948 !! PURPOSE 01949 !! ------- 01950 ! Calculate any freezing/refreezing of liquid water in the snowpack. 01951 ! Also, calculate liquid water transmission and snow runoff. 01952 ! Water flow causes layer thickness reduction and can cause 01953 ! rapid densification of a layer. 01954 ! 01955 ! 01956 USE MODD_CSTS, ONLY : XTT, XLMTT, XRHOLW 01957 USE MODD_SNOW_PAR, ONLY : XSNOWDMIN 01958 ! 01959 USE MODE_SNOW3L 01960 ! 01961 IMPLICIT NONE 01962 ! 01963 !* 0.1 declarations of arguments 01964 ! 01965 REAL, INTENT(IN) :: PTSTEP 01966 ! 01967 REAL, DIMENSION(:), INTENT(IN) :: PRR 01968 ! 01969 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWLIQ, PSNOWRHO 01970 ! 01971 REAL, DIMENSION(:), INTENT(INOUT) :: PTHRUFAL 01972 ! 01973 ! 01974 ! 01975 !* 0.2 declarations of local variables 01976 ! 01977 INTEGER :: JJ, JI 01978 ! 01979 INTEGER :: INI 01980 INTEGER :: INLVLS 01981 ! 01982 REAL, DIMENSION(SIZE(PRR)) :: ZPCPXS, ZTOTWCAP, ZRAINFALL 01983 ! 01984 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZPHASE, ZFLOWLIQ, 01985 ZSNOWLIQ, ZSNOWRHO, 01986 ZWHOLDMAX, ZSNOWDZ, 01987 ZSNOWTEMP, ZSCAP 01988 ! 01989 REAL, DIMENSION(SIZE(PSNOWRHO,1),0:SIZE(PSNOWRHO,2)) :: ZFLOWLIQT 01990 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01991 ! 01992 !------------------------------------------------------------------------------- 01993 ! 01994 ! 0. Initialize: 01995 ! -------------- 01996 ! 01997 IF (LHOOK) CALL DR_HOOK('SNOW3LREFRZ',0,ZHOOK_HANDLE) 01998 ZSNOWRHO(:,:) = PSNOWRHO(:,:) 01999 ZSNOWLIQ(:,:) = PSNOWLIQ(:,:) 02000 ZSNOWTEMP(:,:) = PSNOWTEMP(:,:) 02001 INI = SIZE(PSNOWDZ(:,:),1) 02002 INLVLS = SIZE(PSNOWDZ(:,:),2) 02003 ! 02004 ! 02005 ! 1. Refreeze due to heat transfer 02006 ! ----------------------------- 02007 ! Freeze liquid water in any layers which cooled due 02008 ! to heat transfer. 02009 ! 02010 ZSCAP(:,:) = SNOW3LSCAP(ZSNOWRHO) 02011 ZPHASE(:,:) = MIN(ZSCAP(:,:)* & 02012 MAX(0.0, XTT - ZSNOWTEMP(:,:))*PSNOWDZ(:,:), & 02013 ZSNOWLIQ(:,:)*XLMTT*XRHOLW) 02014 ! 02015 ! 02016 ! Warm layer and reduce liquid if freezing occurs: 02017 ! 02018 ZSNOWDZ(:,:) = MAX(XSNOWDMIN/INLVLS, PSNOWDZ(:,:)) 02019 ! 02020 PSNOWTEMP(:,:) = ZSNOWTEMP(:,:) + ZPHASE(:,:)/(ZSCAP(:,:)*ZSNOWDZ(:,:)) 02021 ! 02022 ! Reduce liquid portion if freezing occurs: 02023 ! 02024 ZSNOWLIQ(:,:) = ZSNOWLIQ(:,:) - ( (PSNOWTEMP(:,:)-ZSNOWTEMP(:,:))* & 02025 ZSCAP(:,:)*ZSNOWDZ(:,:)/(XLMTT*XRHOLW) ) 02026 ! 02027 ! 02028 ! 2. Reduce thickness due to snowmelt in excess of holding capacity 02029 ! -------------------------------------------------------------- 02030 ! Any water in excess of the 02031 ! Maximum holding space for liquid water 02032 ! amount is drained into next layer down. 02033 ! Loss of water due to snowmelt causes a reduction 02034 ! in snow layer mass by a reduction in thickness. 02035 ! 02036 ZWHOLDMAX(:,:) = SNOW3LHOLD(PSNOWRHO,PSNOWDZ) 02037 ! 02038 ZFLOWLIQ(:,:) = MAX(0.,ZSNOWLIQ(:,:)-ZWHOLDMAX(:,:)) 02039 ! 02040 ZSNOWLIQ(:,:) = ZSNOWLIQ(:,:) - ZFLOWLIQ(:,:) 02041 ! 02042 ZSNOWDZ(:,:) = PSNOWDZ(:,:) - ZFLOWLIQ(:,:)*XRHOLW/ZSNOWRHO(:,:) 02043 ! 02044 ZSNOWDZ(:,:) = MAX(0.0, ZSNOWDZ(:,:)) ! to prevent possible very small 02045 ! negative values (machine prescision 02046 ! as snow vanishes 02047 ! 02048 ! 02049 ! 3. Liquid water flow: liquid precipitation and meltwater 02050 ! ----------------------------------------------------- 02051 ! 02052 ! Rainfall flowing into uppermost snow layer: 02053 ! If rainfall is excessive enough (or layers thin enough) 02054 ! it is simply routed directly to runoff: First calculate 02055 ! the total snow pack available liquid water holding capacity: 02056 ! 02057 ZTOTWCAP(:) = 0. 02058 DO JJ=1,INLVLS 02059 DO JI=1,INI 02060 ZTOTWCAP(JI) = ZTOTWCAP(JI) + ZWHOLDMAX(JI,JJ) 02061 ENDDO 02062 ENDDO 02063 ! 02064 ! Rain entering snow (m): 02065 ! 02066 ZRAINFALL(:) = PRR(:)*PTSTEP/XRHOLW ! rainfall (m) 02067 ! 02068 ZFLOWLIQT(:,0)= MIN(ZRAINFALL(:),ZTOTWCAP(:)) 02069 ! 02070 ! Rain assumed to directly pass through the pack to runoff (m): 02071 ! 02072 ZPCPXS(:) = ZRAINFALL(:) - ZFLOWLIQT(:,0) 02073 ! 02074 DO JJ=1,INLVLS 02075 DO JI=1,INI 02076 ZFLOWLIQT(JI,JJ) = ZFLOWLIQ(JI,JJ) 02077 ENDDO 02078 ENDDO 02079 ! 02080 ! 02081 ! Thickness is maintained during water through-flow, 02082 ! so that mass transfer is represented by 02083 ! density changes: NOTE a maximum density 02084 ! is assumed (XRHOSMAX_ES) so that all flow 02085 ! which would result in densities exceeding 02086 ! this limit are routed to next layer down. 02087 ! First test for saturation, then 02088 ! rout excess water down to next layer down 02089 ! and repeat calculation. Net gain in liquid (mass) is 02090 ! translated into a density increase: 02091 ! 02092 ZFLOWLIQ(:,:) = 0.0 ! clear this array for work 02093 PSNOWLIQ(:,:) = ZSNOWLIQ(:,:) ! reset liquid water content 02094 ! 02095 DO JJ=1,INLVLS 02096 DO JI=1,INI 02097 ZSNOWLIQ(JI,JJ) = ZSNOWLIQ(JI,JJ) + ZFLOWLIQT(JI,JJ-1) 02098 ZFLOWLIQ(JI,JJ) = MAX(0.0, ZSNOWLIQ(JI,JJ)-ZWHOLDMAX(JI,JJ)) 02099 ZSNOWLIQ(JI,JJ) = ZSNOWLIQ(JI,JJ) - ZFLOWLIQ(JI,JJ) 02100 ZSNOWRHO(JI,JJ) = ZSNOWRHO(JI,JJ) + (ZSNOWLIQ(JI,JJ) - PSNOWLIQ(JI,JJ))* & 02101 XRHOLW/MAX(XSNOWDMIN/INLVLS,ZSNOWDZ(JI,JJ)) 02102 ZFLOWLIQT(JI,JJ) = ZFLOWLIQT(JI,JJ) + ZFLOWLIQ(JI,JJ) 02103 ENDDO 02104 ENDDO 02105 ! 02106 ! 02107 ! Any remaining throughflow after freezing is available to 02108 ! the soil for infiltration or surface runoff (m). 02109 ! I.E. This is the amount of water leaving the snowpack: 02110 ! Rate water leaves the snowpack [kg/(m2 s)]: 02111 ! 02112 PTHRUFAL(:) = PTHRUFAL(:) + ZFLOWLIQT(:,INLVLS) 02113 ! 02114 ! Add excess rain (rain which flowed directly through the snow 02115 ! due to saturation): 02116 ! 02117 PTHRUFAL(:) = (PTHRUFAL(:) + ZPCPXS(:))*XRHOLW/PTSTEP 02118 ! 02119 ! 02120 ! 4. Update thickness and density and any freezing: 02121 ! ---------------------------------------------- 02122 ! 02123 PSNOWDZ(:,:) = ZSNOWDZ(:,:) 02124 PSNOWRHO(:,:) = ZSNOWRHO(:,:) 02125 PSNOWLIQ(:,:) = ZSNOWLIQ(:,:) 02126 IF (LHOOK) CALL DR_HOOK('SNOW3LREFRZ',1,ZHOOK_HANDLE) 02127 ! 02128 !------------------------------------------------------------------------------- 02129 ! 02130 END SUBROUTINE SNOW3LREFRZ 02131 !#################################################################### 02132 !#################################################################### 02133 !#################################################################### 02134 SUBROUTINE SNOW3LFLUX(PSNOWTEMP,PSNOWDZ,PEXNS,PEXNA, & 02135 PUSTAR2_IC, & 02136 PTSTEP,PALBT,PSW_RAD,PEMIST,PLWUPSNOW, & 02137 PLW_RAD,PTA,PSFCFRZ,PQA,PHPSNOW, & 02138 PSNOWTEMPO1,PSNOWFLUX,PCT,PRADSINK, & 02139 PQSAT,PDQSAT,PRSRA, & 02140 PRN,PH,PGFLUX,PLES3L,PLEL3L,PEVAP, & 02141 PUSTAR,OSFCMELT ) 02142 ! 02143 ! 02144 !! PURPOSE 02145 !! ------- 02146 ! Calculate the surface fluxes (atmospheric/surface). 02147 ! (Noilhan and Planton 1989; Noilhan and Mahfouf 1996) 02148 ! 02149 ! 02150 USE MODD_CSTS,ONLY : XSTEFAN, XCPD, XLSTT, XLVTT, XTT 02151 ! 02152 USE MODE_THERMOS 02153 ! 02154 IMPLICIT NONE 02155 ! 02156 !* 0.1 declarations of arguments 02157 ! 02158 REAL, INTENT(IN) :: PTSTEP 02159 ! 02160 REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ, PSNOWTEMPO1, PSNOWFLUX, PCT, 02161 PRADSINK, PEXNS, PEXNA 02162 ! 02163 REAL, DIMENSION(:), INTENT(IN) :: PALBT, PSW_RAD, PEMIST, PLW_RAD, 02164 PTA, PSFCFRZ, PQA, 02165 PHPSNOW, PQSAT, PDQSAT, PRSRA, 02166 PUSTAR2_IC 02167 ! 02168 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWTEMP 02169 ! 02170 REAL, DIMENSION(:), INTENT(OUT) :: PRN, PH, PGFLUX, PLES3L, PLEL3L, 02171 PEVAP, PLWUPSNOW, PUSTAR 02172 ! 02173 LOGICAL, DIMENSION(:), INTENT(OUT) :: OSFCMELT 02174 ! 02175 ! 02176 !* 0.2 declarations of local variables 02177 ! 02178 REAL, DIMENSION(SIZE(PSNOWDZ)) :: ZEVAPC, ZLE, ZSNOWTEMP, ZSMSNOW, ZGFLUX, 02179 ZDELTAT, ZSNOWTO3 02180 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02181 ! 02182 !------------------------------------------------------------------------------- 02183 ! 02184 ! 0. Initialize: 02185 ! -------------- 02186 ! 02187 IF (LHOOK) CALL DR_HOOK('SNOW3LFLUX',0,ZHOOK_HANDLE) 02188 ZSNOWTEMP(:) = PSNOWTEMP(:) 02189 ZLE(:) = 0.0 02190 ZSMSNOW(:) = 0.0 02191 ZGFLUX(:) = 0.0 02192 ! 02193 OSFCMELT(:) = .FALSE. 02194 ! 02195 ZSNOWTO3(:) = PSNOWTEMPO1(:) ** 3 ! to save some CPU time, store this 02196 ! 02197 ! 1. Flux calculations when melt not occuring at surface (W/m2): 02198 ! -------------------------------------------------------------- 02199 ! 02200 ! 02201 ZDELTAT(:) = PSNOWTEMP(:) - PSNOWTEMPO1(:) ! surface T time change: 02202 ! 02203 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * ZSNOWTO3(:)*( PSNOWTEMPO1(:) + 4.* ZDELTAT(:) ) 02204 ! 02205 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02206 ! 02207 PH(:) = PRSRA(:) * XCPD * (PSNOWTEMP(:)/PEXNS(:) - PTA(:)/PEXNA(:)) 02208 ! 02209 ZEVAPC(:) = PRSRA(:) * ( (PQSAT(:) - PQA(:)) + PDQSAT(:)*ZDELTAT(:) ) 02210 ! 02211 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02212 ! 02213 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02214 ! 02215 ZLE(:) = PLES3L(:) + PLEL3L(:) 02216 ! 02217 PGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02218 ! 02219 ! 02220 ! 2. Initial melt adjustment 02221 ! -------------------------- 02222 ! If energy avalabile to melt snow, then recalculate fluxes 02223 ! at the freezing point and add residual heat to layer 02224 ! average heat. 02225 ! 02226 ! A) If temperature change is > 0 and passes freezing point this timestep, 02227 ! then recalculate fluxes at freezing point and excess energy 02228 ! will be used outside of this routine to change snow heat content: 02229 ! 02230 WHERE (PSNOWTEMP > XTT .AND. PSNOWTEMPO1 < XTT) 02231 ! 02232 OSFCMELT(:)= .TRUE. 02233 ! 02234 ZDELTAT(:) = XTT - PSNOWTEMPO1(:) 02235 ! 02236 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * ZSNOWTO3(:)*( PSNOWTEMPO1(:) + 4.* ZDELTAT(:) ) 02237 ! 02238 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02239 ! 02240 PH(:) = PRSRA(:) * XCPD * (XTT/PEXNS(:) - PTA(:)/PEXNA(:)) 02241 ! 02242 ZEVAPC(:) = PRSRA(:) * ( (PQSAT(:) - PQA(:)) + PDQSAT(:)*ZDELTAT(:) ) 02243 ! 02244 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02245 ! 02246 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02247 ! 02248 ZLE(:) = PLES3L(:) + PLEL3L(:) 02249 ! 02250 ZGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02251 ! 02252 ZSMSNOW(:) = PGFLUX(:) - ZGFLUX(:) 02253 ! 02254 PGFLUX(:) = ZGFLUX(:) 02255 ! 02256 ! This will be used to change heat content of snow: 02257 ! 02258 ZSNOWTEMP(:) = PSNOWTEMP(:) - ZSMSNOW(:)*PTSTEP*PCT(:) 02259 ! 02260 END WHERE 02261 ! 02262 ! 3. Ongoing melt adjustment: explicit solution 02263 ! --------------------------------------------- 02264 ! If temperature change is 0 and at freezing point this timestep, 02265 ! then recalculate fluxes and surface temperature *explicitly* 02266 ! as this is *exact* for snow at freezing point (Brun, Martin) 02267 ! 02268 WHERE(PSNOWTEMP(:) > XTT .AND. PSNOWTEMPO1(:) >= XTT) 02269 ! 02270 OSFCMELT(:) = .TRUE. 02271 ! 02272 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * (XTT ** 4) 02273 ! 02274 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02275 ! 02276 PH(:) = PRSRA(:) * XCPD * (XTT/PEXNS(:) - PTA(:)/PEXNA(:)) 02277 ! 02278 ZEVAPC(:) = PRSRA(:) * (PQSAT(:) - PQA(:)) 02279 ! 02280 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02281 ! 02282 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02283 ! 02284 ZLE(:) = PLES3L(:) + PLEL3L(:) 02285 ! 02286 PGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02287 ! 02288 ZSNOWTEMP(:) = XTT + PTSTEP*PCT(:)*(PGFLUX(:) + PRADSINK(:) - PSNOWFLUX(:)) 02289 ! 02290 END WHERE 02291 ! 02292 ! 4. Update surface temperature: 02293 ! ------------------------------ 02294 ! 02295 PSNOWTEMP(:) = ZSNOWTEMP(:) 02296 ! 02297 ! 5. Final evaporative flux (kg/m2/s) 02298 ! 02299 PEVAP(:) = ZEVAPC(:) 02300 ! 02301 ! 6. Friction velocity 02302 ! -------------------- 02303 ! 02304 PUSTAR(:) = SQRT(PUSTAR2_IC(:)) 02305 ! 02306 IF (LHOOK) CALL DR_HOOK('SNOW3LFLUX',1,ZHOOK_HANDLE) 02307 ! 02308 !------------------------------------------------------------------------------- 02309 ! 02310 END SUBROUTINE SNOW3LFLUX 02311 !#################################################################### 02312 !#################################################################### 02313 !#################################################################### 02314 SUBROUTINE SNOW3LEVAPN(PPSN3L,PLES3L,PLEL3L,PTSTEP,PSNOWTEMP, & 02315 PSNOWRHO,PSNOWDZ,PSNOWLIQ,PEVAPCOR) 02316 ! 02317 ! 02318 !! PURPOSE 02319 !! ------- 02320 ! Remove mass from uppermost snow layer in response to 02321 ! evaporation (liquid) and sublimation. 02322 ! 02323 ! 02324 USE MODD_CSTS, ONLY : XLVTT, XRHOLW, XLSTT, XLMTT, XCI, XTT 02325 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN 02326 ! 02327 IMPLICIT NONE 02328 ! 02329 !* 0.1 declarations of arguments 02330 ! 02331 REAL, INTENT(IN) :: PTSTEP 02332 ! 02333 REAL, DIMENSION(:), INTENT(IN) :: PPSN3L 02334 ! 02335 REAL, DIMENSION(:), INTENT(IN) :: PLES3L, PLEL3L ! (W/m2) 02336 ! 02337 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ, PSNOWLIQ, 02338 PEVAPCOR, PSNOWTEMP 02339 ! 02340 !* 0.2 declarations of local variables 02341 ! 02342 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWEVAPS, ZSNOWEVAP, ZSNOWEVAPX, 02343 ZSNOWDZ, ZEVAPCOR, ZSNOWHEAT, ZSCAP 02344 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02345 ! ZEVAPCOR = for vanishingy thin snow cover, 02346 ! allow any excess evaporation 02347 ! to be extracted from the soil 02348 ! to maintain an accurate water 02349 ! balance [kg/(m2 s)] 02350 ! 02351 ! 02352 !------------------------------------------------------------------------------- 02353 ! 02354 ! 0. Initialize: 02355 ! -------------- 02356 ! 02357 IF (LHOOK) CALL DR_HOOK('SNOW3LEVAPN',0,ZHOOK_HANDLE) 02358 ZEVAPCOR(:) = 0.0 02359 ZSNOWEVAPS(:) = 0.0 02360 ZSNOWEVAP(:) = 0.0 02361 ZSNOWEVAPX(:) = 0.0 02362 ZSNOWDZ(:) = 0.0 02363 ZSCAP(:) = 0.0 02364 ZSNOWHEAT(:) = 0.0 02365 ! 02366 ! 02367 ! 02368 WHERE(PSNOWDZ > 0.0) 02369 ! 02370 ! 1. Evaporation of snow liquid water 02371 ! ----------------------------------- 02372 ! Evaporation reduces liq water equivalent content 02373 ! of snow pack either by reducing the snow density 02374 ! (evaporation) or the layer thickness (sublimation). 02375 ! Condensation does the reverse. 02376 ! 02377 ! Multiply evaporation components by snow fraction 02378 ! to be consistent with fluxes from the snow covered 02379 ! fraction of grid box 02380 ! 02381 ! Evaporation: 02382 ! Reduce density and liquid water content: 02383 ! 02384 ZSNOWEVAP(:) = PPSN3L(:)*PLEL3L(:)*PTSTEP/(XLVTT*XRHOLW) 02385 ZSNOWEVAPX(:) = MIN(ZSNOWEVAP(:),PSNOWLIQ(:)) 02386 ! 02387 ! 02388 PSNOWRHO(:) = PSNOWRHO(:) - ZSNOWEVAPX(:)*XRHOLW/PSNOWDZ(:) 02389 PSNOWLIQ(:) = PSNOWLIQ(:) - ZSNOWEVAPX(:) 02390 ! 02391 ! 02392 ! Budget check: If density drops below minimum, then extract the 02393 ! corresponding water mass from soil (for vanishingly thin snow covers): 02394 ! 02395 ZEVAPCOR(:) = MAX(0.0,XRHOSMIN_ES-PSNOWRHO(:))*PSNOWDZ(:)/PTSTEP 02396 PSNOWRHO(:) = MAX(XRHOSMIN_ES,PSNOWRHO(:)) 02397 ! 02398 ! 02399 ! Budget check: as last traces of liquid in snow evaporates, it is possible 02400 ! evaporation could exceed liquid present in snow. If this is the case, 02401 ! remove residual mass from solid snow in order to maintain high 02402 ! accuracy water balance: 02403 ! 02404 ZSNOWEVAPX(:) = MAX(0.0, ZSNOWEVAP(:) - ZSNOWEVAPX(:)) 02405 ZSNOWDZ(:) = PSNOWDZ(:) - ZSNOWEVAPX(:)*XRHOLW/PSNOWRHO(:) 02406 PSNOWDZ(:) = MAX(0.0, ZSNOWDZ(:)) 02407 ZEVAPCOR(:) = ZEVAPCOR(:) + MAX(0.0,-ZSNOWDZ(:))*PSNOWRHO(:)/PTSTEP 02408 ! 02409 !! 2. Sublimation of snow ice 02410 ! --------------------------- 02411 ! Reduce layer thickness and total snow depth 02412 ! if sublimation: add to correction term if potential 02413 ! sublimation exceeds available snow cover. 02414 ! 02415 ZSNOWEVAPS(:) = PPSN3L(:)*PLES3L(:)*PTSTEP/(XLSTT*PSNOWRHO(:)) 02416 ZSNOWDZ(:) = PSNOWDZ(:) - ZSNOWEVAPS(:) 02417 PSNOWDZ(:) = MAX(0.0, ZSNOWDZ(:)) 02418 ZEVAPCOR(:) = ZEVAPCOR(:) + MAX(0.0,-ZSNOWDZ(:))*PSNOWRHO(:)/PTSTEP 02419 ! 02420 ! Adjust heat content considering mass loss...NOTE, temperature does not 02421 ! change unless *possibly* if liquid present in the snowpack 02422 ! 02423 ZSCAP(:) = SNOW3LSCAP(PSNOWRHO) 02424 ZSNOWHEAT(:) = PSNOWDZ(:)*( ZSCAP(:)*(PSNOWTEMP(:)-XTT) & 02425 - XLMTT*PSNOWRHO(:) ) + XLMTT*XRHOLW*PSNOWLIQ(:) 02426 ! 02427 PSNOWTEMP(:) = XTT + ( ((ZSNOWHEAT(:)/MAX(XSNOWDMIN,PSNOWDZ(:))) & 02428 + XLMTT*PSNOWRHO(:))/ZSCAP(:) ) 02429 ! 02430 PSNOWLIQ(:) = MAX(0.0,PSNOWTEMP(:)-XTT)*ZSCAP(:)* & 02431 PSNOWDZ(:)/(XLMTT*XRHOLW) 02432 ! 02433 PSNOWTEMP(:) = MIN(XTT,PSNOWTEMP(:)) 02434 ! 02435 END WHERE 02436 ! 02437 ! 3. Update evaporation correction term: 02438 ! -------------------------------------- 02439 ! 02440 PEVAPCOR(:) = PEVAPCOR(:) + ZEVAPCOR(:) 02441 IF (LHOOK) CALL DR_HOOK('SNOW3LEVAPN',1,ZHOOK_HANDLE) 02442 ! 02443 !------------------------------------------------------------------------------- 02444 ! 02445 END SUBROUTINE SNOW3LEVAPN 02446 !#################################################################### 02447 !#################################################################### 02448 !#################################################################### 02449 SUBROUTINE SNOW3LGONE(PTSTEP,PPSN3L,PLEL3L,PLES3L,PSNOWRHO, & 02450 PSNOWHEAT,PRADSINK,PEVAPCOR,PTHRUFAL,PGRNDFLUX, & 02451 PGFLUXSNOW,PGRNDFLUXO,PSNOWDZ,PSNOWLIQ,PSNOWTEMP, & 02452 PRADXS) 02453 ! 02454 ! 02455 ! 02456 !! PURPOSE 02457 !! ------- 02458 ! Account for the case when the last trace of snow melts 02459 ! during a time step: ensure mass and heat balance of 02460 ! snow AND underlying surface. 02461 ! 02462 ! 02463 USE MODD_CSTS,ONLY : XTT, XLSTT, XLVTT 02464 ! 02465 IMPLICIT NONE 02466 ! 02467 !* 0.1 declarations of arguments 02468 ! 02469 REAL, INTENT(IN) :: PTSTEP 02470 ! 02471 REAL, DIMENSION(:), INTENT(IN) :: PPSN3L, PLEL3L, PLES3L, PGFLUXSNOW, 02472 PRADSINK, PGRNDFLUXO 02473 ! 02474 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWHEAT 02475 ! 02476 REAL, DIMENSION(:), INTENT(INOUT) :: PGRNDFLUX, PRADXS 02477 ! 02478 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWLIQ, PSNOWTEMP 02479 ! 02480 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL ! melt water [kg/(m2 s)] 02481 ! 02482 REAL, DIMENSION(:), INTENT(OUT) :: PEVAPCOR ! [kg/(m2 s)] 02483 ! PEVAPCOR = for vanishingy thin snow cover, 02484 ! allow any excess evaporation 02485 ! to be extracted from the soil 02486 ! to maintain an accurate water 02487 ! balance. 02488 ! 02489 !* 0.2 declarations of local variables 02490 ! 02491 INTEGER :: JJ, JI 02492 ! 02493 INTEGER :: INI 02494 INTEGER :: INLVLS 02495 ! 02496 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWHEATC, ZSNOWGONE_DELTA 02497 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02498 ! 02499 !------------------------------------------------------------------------------- 02500 ! 02501 ! 0. Initialize: 02502 ! -------------- 02503 ! 02504 IF (LHOOK) CALL DR_HOOK('SNOW3LGONE',0,ZHOOK_HANDLE) 02505 PEVAPCOR(:) = 0.0 02506 PTHRUFAL(:) = 0.0 02507 INI = SIZE(PSNOWDZ(:,:),1) 02508 INLVLS = SIZE(PSNOWDZ(:,:),2) 02509 ZSNOWHEATC(:) = 0. 02510 DO JJ=1,INLVLS 02511 DO JI=1,INI 02512 ZSNOWHEATC(JI) = ZSNOWHEATC(JI) + PSNOWHEAT(JI,JJ) ! total heat content (J m-2) 02513 ENDDO 02514 ENDDO 02515 ZSNOWGONE_DELTA(:) = 1.0 02516 ! 02517 ! 1. Simple test to see if snow vanishes: 02518 ! --------------------------------------- 02519 ! If so, set thicknesses (and therefore mass and heat) and liquid content 02520 ! to zero, and adjust fluxes of water, evaporation and heat into underlying 02521 ! surface. Note, test with flux computed *before* correction since this represents 02522 ! actual inflow of heat from below (as heat content correction owing to a corrected 02523 ! flux has not yet been done: here we compare to pre-corrected heat content). 02524 ! 02525 WHERE(PGFLUXSNOW(:) - PGRNDFLUXO(:) + PRADSINK(:) >= (-ZSNOWHEATC(:)/PTSTEP) ) 02526 PGRNDFLUX(:) = PGFLUXSNOW(:) + (ZSNOWHEATC(:)/PTSTEP) 02527 PEVAPCOR(:) = PPSN3L(:)*((PLEL3L(:)/XLVTT) + (PLES3L(:)/(XLSTT*PSNOWRHO(:,1)))) 02528 PRADXS(:) = 0.0 02529 ZSNOWGONE_DELTA(:) = 0.0 ! FLAG...if=0 then snow vanishes, else=1 02530 END WHERE 02531 ! 02532 PTHRUFAL(:) = 0. 02533 DO JJ=1,INLVLS 02534 DO JI=1,INI 02535 PTHRUFAL(JI) = PTHRUFAL(JI) + (1.0-ZSNOWGONE_DELTA(JI))*PSNOWRHO(JI,JJ)*PSNOWDZ(JI,JJ)/PTSTEP 02536 ENDDO 02537 END DO 02538 ! 02539 ! 2. Final update of snow state 02540 ! ----------------------------- 02541 ! (either still present or not): 02542 ! 02543 DO JJ=1,INLVLS 02544 DO JI=1,INI 02545 PSNOWDZ (JI,JJ) = PSNOWDZ (JI,JJ)*ZSNOWGONE_DELTA(JI) 02546 PSNOWLIQ (JI,JJ) = PSNOWLIQ (JI,JJ)*ZSNOWGONE_DELTA(JI) 02547 PSNOWTEMP(JI,JJ) = (1.0-ZSNOWGONE_DELTA(JI))*XTT + PSNOWTEMP(JI,JJ)*ZSNOWGONE_DELTA(JI) 02548 ENDDO 02549 ENDDO 02550 IF (LHOOK) CALL DR_HOOK('SNOW3LGONE',1,ZHOOK_HANDLE) 02551 ! 02552 ! 02553 END SUBROUTINE SNOW3LGONE 02554 !#################################################################### 02555 !#################################################################### 02556 !#################################################################### 02557 SUBROUTINE SNOW3LEVAPGONE(PSNOWHEAT,PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ) 02558 ! 02559 !! PURPOSE 02560 !! ------- 02561 ! 02562 ! If all snow in uppermost layer evaporates/sublimates, re-distribute 02563 ! grid (below assumes very thin snowpacks so layer-thicknesses are 02564 ! constant). 02565 ! 02566 ! 02567 USE MODD_CSTS, ONLY : XTT, XRHOLW, XLMTT 02568 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN, XRHOSMAX_ES 02569 ! 02570 IMPLICIT NONE 02571 ! 02572 !* 0.1 declarations of arguments 02573 ! 02574 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO ! snow density profile (kg/m3) 02575 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ ! snow layer thickness profile (m) 02576 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT ! snow heat content/enthalpy (J/m2) 02577 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP ! snow temperature profile (K) 02578 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWLIQ ! snow liquid water profile (m) 02579 ! 02580 !* 0.2 declarations of local variables 02581 ! 02582 INTEGER :: JJ, JI 02583 ! 02584 INTEGER :: INI 02585 INTEGER :: INLVLS 02586 ! 02587 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOWHEAT_1D ! total heat content (J/m2) 02588 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOW ! total snow depth (m) 02589 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZMASS ! total mass 02590 ! 02591 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSCAP ! Snow layer heat capacity (J/K/m3) 02592 ! 02593 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02594 ! 02595 !------------------------------------------------------------------------------- 02596 ! 02597 ! Initialize: 02598 ! 02599 IF (LHOOK) CALL DR_HOOK('SNOW3LEVAPGONE',0,ZHOOK_HANDLE) 02600 INI = SIZE(PSNOWDZ,1) 02601 INLVLS = SIZE(PSNOWDZ,2) 02602 ! 02603 ! First, determine where uppermost snow layer has completely 02604 ! evaporated/sublimated (as it becomes thin): 02605 ! 02606 ZSNOWHEAT_1D(:) = 0. 02607 ZSNOW(:) = 0. 02608 ZMASS(:) = 0. 02609 ! 02610 ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO(:,:)) 02611 ! 02612 DO JJ=2,INLVLS 02613 DO JI=1,INI 02614 IF(PSNOWDZ(JI,1) == 0.0)THEN 02615 ZSNOWHEAT_1D(JI) = ZSNOWHEAT_1D(JI) + XLMTT*XRHOLW*PSNOWLIQ(JI,JJ) & 02616 + PSNOWDZ(JI,JJ)*(ZSCAP(JI,JJ)*(PSNOWTEMP(JI,JJ)-XTT) & 02617 - XLMTT*PSNOWRHO(JI,JJ) ) 02618 ZSNOW (JI) = ZSNOW(JI) + PSNOWDZ(JI,JJ) 02619 ZMASS (JI) = ZMASS(JI) + PSNOWDZ(JI,JJ)*PSNOWRHO(JI,JJ) 02620 ENDIF 02621 ENDDO 02622 ENDDO 02623 ! 02624 ! Where uppermost snow layer has vanished, redistribute vertical 02625 ! snow mass and heat profiles (and associated quantities): 02626 ! 02627 DO JJ=1,INLVLS 02628 DO JI=1,INI 02629 IF(ZSNOW(JI)/= 0.0)THEN 02630 ZSNOW (JI) = MAX(0.5*XSNOWDMIN,ZSNOW(JI)) 02631 PSNOWDZ (JI,JJ) = ZSNOW(JI)/REAL(INLVLS) 02632 PSNOWHEAT(JI,JJ) = ZSNOWHEAT_1D(JI)/REAL(INLVLS) 02633 PSNOWRHO (JI,JJ) = ZMASS (JI)/ZSNOW(JI) 02634 ENDIF 02635 ENDDO 02636 ENDDO 02637 ! 02638 ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO(:,:)) 02639 ! 02640 DO JJ=1,INLVLS 02641 DO JI=1,INI 02642 IF(ZSNOW(JI)/= 0.0)THEN 02643 PSNOWTEMP(JI,JJ) = XTT + ( ((PSNOWHEAT(JI,JJ)/PSNOWDZ(JI,JJ)) & 02644 + XLMTT*PSNOWRHO(JI,JJ))/ZSCAP(JI,JJ) ) 02645 PSNOWLIQ (JI,JJ) = MAX(0.0,PSNOWTEMP(JI,JJ)-XTT)*ZSCAP(JI,JJ)* & 02646 PSNOWDZ(JI,JJ)/(XLMTT*XRHOLW) 02647 PSNOWTEMP(JI,JJ) = MIN(XTT,PSNOWTEMP(JI,JJ)) 02648 ENDIF 02649 ENDDO 02650 ENDDO 02651 IF (LHOOK) CALL DR_HOOK('SNOW3LEVAPGONE',1,ZHOOK_HANDLE) 02652 ! 02653 END SUBROUTINE SNOW3LEVAPGONE 02654 ! 02655 ! 02656 ! 02657 END SUBROUTINE SNOW3L
1.8.0
