|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ########################################################################## 00002 SUBROUTINE SNOWCRO(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 !!**** *SNOWCRO* 00020 !! 00021 !! PURPOSE 00022 !! ------- 00023 ! 00024 ! Detailed snowpack scheme Crocus, computationnally based on the 00025 ! 3-Layer snow scheme option (Boone and Etchevers 1999) 00026 ! For shallow snow cover, Default method of Douville et al. (1995) 00027 ! used with this option: Model "turns on" when snow sufficiently 00028 ! deep/above some preset critical snow depth. 00029 ! 00030 ! 00031 ! 00032 ! 00033 !!** METHOD 00034 !! ------ 00035 ! 00036 ! Direct calculation 00037 ! 00038 !! EXTERNAL 00039 !! -------- 00040 ! 00041 ! None 00042 !! 00043 !! IMPLICIT ARGUMENTS 00044 !! ------------------ 00045 !! 00046 !! 00047 !! 00048 !! REFERENCE 00049 !! --------- 00050 !! 00051 !! ISBA: Belair (1995) 00052 !! ISBA: Noilhan and Planton (1989) 00053 !! ISBA: Noilhan and Mahfouf (1996) 00054 !! ISBA-ES: Boone and Etchevers (2001) 00055 !! Crocus : Brun et al., 1989 (J. Glaciol.) 00056 !! Crocus : Brun et al., 1992 (J. Glaciol.) 00057 !! Crocus : Vionnet et al., in prep (Geosci. Mod. Devel. Discuss.) 00058 !! 00059 !! 00060 !! AUTHOR 00061 !! ------ 00062 !! A. Boone * Meteo-France * 00063 !! V. Vionnet * Meteo-France * 00064 !! E. Brun * Meteo-France * 00065 !! 00066 !! MODIFICATIONS 00067 !! ------------- 00068 !! Original 7/99 00069 !! Modified by A.Boone 05/02 (code, not physics) 00070 !! Modified by A.Boone 11/04 i) maximum density limit imposed (although 00071 !! rarely if ever reached), ii) check to 00072 !! see if upermost layer completely sublimates 00073 !! during a timestep (as snowpack becomes vanishly 00074 !! thin), iii) impose maximum grain size limit 00075 !! in radiation transmission computation. 00076 !! 00077 !! Modified by B. Decharme (03/2009): Consistency with Arpege permanent 00078 !! snow/ice treatment (LGLACIER for alb) 00079 !! Modified by A. Boone (04/2010): Implicit coupling and replace Qsat and DQsat 00080 !! by Qsati and DQsati, respectively. 00081 !! Modified by E. Brun, V. Vionnet, S. Morin (05/2011): 00082 !! Addition of Crocus processes and 00083 !! parametrizations to 00084 !! the SNOW-3L code. This includes the dynamic handling 00085 !! of snow layers and the inclusion of snow metamorphism 00086 !! rules similar to the original Crocus implementation. 00087 !! Modified by B. Decharme (09/2012): New wind implicitation 00088 !! 00089 !! Modified by M. Lafaysse (07/2012) : 00090 !! * Albedo and roughness parametrizations 00091 !! for surface ice over glaciers 00092 !! MODIF 2012-10-03 : don't modify roughness if implicit coupling 00093 !! (test PPEW_A_COEF == 0. ) 00094 !! * SNOWCROALB is now called by SNOWCRORAD to remove duplicated code 00095 !! * Parameters for albedo are moved to modd_snow_par 00096 !! * PSNOWAGE is stored as an age 00097 !! (days since snowfall) and not as a date 00098 !! to allow spinup simulations 00099 !! * New rules for optimal discretization of very thick snowpacks 00100 !! * Optional outputs for debugging 00101 !! 00102 !! Modified by E. Brun and M. Lafaysse (07/2012) : 00103 !! * Implement sublimation in SNOWDRIFT 00104 !! * Flag in namelist to activate SNOWDRIFT and SNOWDRIFT_SUBLIM 00105 !! Modified by E. Brun and M. Lafaysse (08/2012) : 00106 !! * UEPSI replaced by 0 in the if statement of case 1.3.3.2 (SNOWCROMETAMO) 00107 !! * If SNOWDRIFT is activated the wind do not modify grain types during snowfall 00108 !! (redundant with snowdrift) 00109 !! Modified by E. Brun (24/09/2012) : 00110 !! * Correction coupling coefficient for specific humidity in SNOWCROEBUD 00111 !! * PSFCFRZ(:) = 1.0 for systematic solid/vapor latent fluxes in SNOWCROEBUD 00112 00113 00114 00115 !------------------------------------------------------------------------------- 00116 ! 00117 !* 0. DECLARATIONS 00118 ! ------------ 00119 USE MODD_CSTS, ONLY : XTT, XRHOLW, XLMTT,XLSTT,XLVTT, XCL, XCI, XPI 00120 ! 00121 USE MODE_SNOW3L 00122 ! 00123 USE MODD_SNOW_PAR, ONLY : XZ0ICEZ0SNOW, XRHOTHRESHOLD_ICE, LSNOWDRIFT 00124 ! 00125 USE MODD_SNOW_METAMO 00126 ! 00127 USE MODD_TYPE_DATE_SURF, ONLY: DATE_TIME 00128 ! 00129 USE MODI_ABOR1_SFX 00130 ! 00131 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00132 USE PARKIND1 ,ONLY : JPRB 00133 ! 00134 ! this module is not used anymore 00135 ! USE MODI_GREGODSTRATI 00136 00137 USE MODE_CRODEBUG 00138 00139 IMPLICIT NONE 00140 ! 00141 ! 00142 !* 0.1 declarations of arguments 00143 ! 00144 REAL, INTENT(IN) :: PTSTEP 00145 ! PTSTEP = time step of the integration 00146 TYPE(DATE_TIME), INTENT(IN) :: TPTIME ! current date and time 00147 ! 00148 CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES 00149 ! HSNOWRES = ISBA-SNOW3L turbulant exchange option 00150 ! 'DEF' = Default: Louis (ISBA: Noilhan and Mahfouf 1996) 00151 ! 'RIL' = Limit Richarson number under very stable 00152 ! conditions (currently testing) 00153 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 00154 ! initialise WGI=WSAT, 00155 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 00156 ! False = No specific treatment 00157 ! 00158 CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option 00159 ! ! 'OLD' = direct 00160 ! ! 'NEW' = Taylor serie, order 1 00161 ! 00162 REAL, DIMENSION(:), INTENT(IN) :: PPS, PTA, PSW_RAD, PQA, 00163 PVMOD, PLW_RAD, PSR, PRR 00164 ! PSW_RAD = incoming solar radiation (W/m2) 00165 ! PLW_RAD = atmospheric infrared radiation (W/m2) 00166 ! PRR = rain rate [kg/(m2 s)] 00167 ! PSR = snow rate (SWE) [kg/(m2 s)] 00168 ! PTA = atmospheric temperature at level za (K) 00169 ! PVMOD = modulus of the wind parallel to the orography (m/s) 00170 ! PPS = surface pressure 00171 ! PQA = atmospheric specific humidity 00172 ! at level za 00173 ! 00174 REAL, DIMENSION(:), INTENT(IN) :: PTG, PSOILCOND, PD_G, PPSN3L 00175 ! PTG = Surface soil temperature (effective 00176 ! temperature the of layer lying below snow) 00177 ! PSOILCOND = soil thermal conductivity [W/(m K)] 00178 ! PD_G = Assumed first soil layer thickness (m) 00179 ! Used to calculate ground/snow heat flux 00180 ! PPSN3L = snow fraction 00181 ! 00182 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PUREF, PEXNS, PEXNA, PDIRCOSZW, PRHOA, PZ0, PZ0EFF, 00183 PALB, PZ0H, PPERMSNOWFRAC 00184 ! PZ0EFF = roughness length for momentum 00185 ! PZ0 = grid box average roughness length 00186 ! PZ0H = grid box average roughness length for heat 00187 ! PZREF = reference height of the first 00188 ! atmospheric level 00189 ! PUREF = reference height of the wind 00190 ! PRHOA = air density 00191 ! PEXNS = Exner function at surface 00192 ! PEXNA = Exner function at lowest atmos level 00193 ! PDIRCOSZW = Cosinus of the angle between the 00194 ! normal to the surface and the vertical 00195 ! PALB = soil/vegetation albedo 00196 ! PPERMSNOWFRAC = fraction of permanet snow/ice 00197 ! 00198 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, 00199 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, 00200 PPEQ_B_COEF 00201 ! PPEW_A_COEF = wind coefficient (m2s/kg) 00202 ! PPEW_B_COEF = wind coefficient (m/s) 00203 ! PPET_A_COEF = A-air temperature coefficient 00204 ! PPET_B_COEF = B-air temperature coefficient 00205 ! PPEQ_A_COEF = A-air specific humidity coefficient 00206 ! PPEQ_B_COEF = B-air specific humidity coefficient 00207 ! 00208 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWALB 00209 ! PSNOWALB = Prognostic surface snow albedo 00210 ! (does not include anything but 00211 ! the actual snow cover) 00212 ! 00213 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWHEAT, PSNOWRHO, PSNOWSWE 00214 ! PSNOWHEAT = Snow layer(s) heat content (J/m2) 00215 ! PSNOWRHO = Snow layer(s) averaged density (kg/m3) 00216 ! PSNOWSWE = Snow layer(s) liquid Water Equivalent (SWE:kg m-2) 00217 ! 00218 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST 00219 ! PSNOWGRAN1 = Snow layers grain feature 1 00220 ! PSNOWGRAN2 = Snow layer grain feature 2 00221 ! PSNOWHIST = Snow layer grain historical 00222 ! parameter (only for non 00223 ! dendritic snow) 00224 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWAGE ! Snow grain age 00225 ! 00226 REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWLIQ, PSNOWTEMP, PSNOWDZ 00227 ! PSNOWLIQ = Snow layer(s) liquid water content (m) 00228 ! PSNOWTEMP = Snow layer(s) temperature (m) 00229 ! PSNOWDZ = Snow layer(s) thickness (m) 00230 ! 00231 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL, PGRNDFLUX, PEVAPCOR 00232 ! PTHRUFAL = rate that liquid water leaves snow pack: 00233 ! paritioned into soil infiltration/runoff 00234 ! by ISBA [kg/(m2 s)] 00235 ! PGRNDFLUX = soil/snow interface heat flux (W/m2) 00236 ! PEVAPCOR = evaporation/sublimation correction term: 00237 ! extract any evaporation exceeding the 00238 ! actual snow cover (as snow vanishes) 00239 ! and apply it as a surface soil water 00240 ! sink. [kg/(m2 s)] 00241 ! 00242 REAL, DIMENSION(:), INTENT(OUT) :: PRNSNOW, PHSNOW, PGFLUXSNOW, PLES3L, PLEL3L, 00243 PHPSNOW, PCDSNOW, PUSTAR, PEVAP 00244 ! PLES3L = evaporation heat flux from snow (W/m2) 00245 ! PLEL3L = sublimation (W/m2) 00246 ! PHPSNOW = heat release from rainfall (W/m2) 00247 ! PRNSNOW = net radiative flux from snow (W/m2) 00248 ! PHSNOW = sensible heat flux from snow (W/m2) 00249 ! PGFLUXSNOW = net heat flux from snow (W/m2) 00250 ! PCDSNOW = drag coefficient for momentum over snow 00251 ! PUSTAR = friction velocity over snow (m/s) 00252 ! PEVAP = total evaporative flux (kg/m2/s) 00253 ! 00254 REAL, DIMENSION(:), INTENT(OUT) :: PCHSNOW, PEMISNOW, PSNOWHMASS 00255 ! PEMISNOW = snow surface emissivity 00256 ! PCHSNOW = drag coefficient for heat over snow 00257 ! PSNOWHMASS = heat content change due to mass 00258 ! changes in snowpack (J/m2): for budget 00259 ! calculations only. 00260 ! 00261 REAL, DIMENSION(:), INTENT(OUT) :: PRI 00262 ! PRI = Ridcharson number 00263 ! 00264 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle 00265 REAL, DIMENSION(:), INTENT(IN) :: PXLAT,PXLON ! LAT/LON after packing 00266 ! 00267 ! 00268 !* 0.2 declarations of local variables 00269 ! 00270 ! 00271 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWTEMP, ZSCAP, ZSNOWDZN, ZSCOND, 00272 ZRADSINK 00273 ! ZSNOWTEMP = Snow layer(s) averaged temperature (K) 00274 ! ZSCAP = Snow layer(s) heat capacity [J/(K m3)] 00275 ! ZSNOWDZN = Updated snow layer thicknesses (m) 00276 ! ZSCOND = Snow layer(s) thermal conducivity [W/(m K)] 00277 ! ZRADSINK = Snow solar Radiation source terms (W/m2) 00278 ! 00279 ! 00280 REAL, DIMENSION(SIZE(PTA)) :: ZSNOWBIS 00281 ! ZSNOWBIS = Total snow depth after snowfall 00282 ! 00283 REAL, DIMENSION(SIZE(PTA)) :: ZSNOW, ZSFCFRZ, ZTSTERM1, ZTSTERM2, 00284 ZCT, ZRA, ZSNOWFLUX, ZSNOWTEMPO1 00285 ! ZSNOW = Total snow depth (m) 00286 ! ZCT = inverse of the product of snow heat capacity 00287 ! and layer thickness [(m2 K)/J] 00288 ! ZRA = Surface aerodynamic resistance 00289 ! ZTSTERM1,ZTSTERM2 = Surface energy budget coefficients 00290 ! ZSNOWFLUX = heat flux between 1st and 2nd snow layers: 00291 ! used during surface melting (W/m2) 00292 ! ZSNOWTEMPO1= value of uppermost snow temperature 00293 ! before time integration (K) 00294 ! 00295 LOGICAL, DIMENSION(SIZE(PTA)) :: GSNOWFALL,GMODIF_MAILLAGE 00296 ! GSNOWFALL = FLAG if snowfall exceed PSNOW/10, used for 00297 ! grid updating. 00298 ! 00299 REAL, DIMENSION(SIZE(PTA)) :: ZRSRA, ZDQSAT, ZQSAT, ZRADXS, ZLIQHEATXS, 00300 ZLWUPSNOW 00301 ! ZRSRA = air density over aerodynamic resistance 00302 ! ZDQSAT = derrivative of saturation specific humidity 00303 ! ZQSAT = saturation specific humidity 00304 ! ZRADXS = shortwave radiation absorbed by soil surface 00305 ! (for thin snow sover) (W m-2) 00306 ! ZLIQHEATXS = excess snowpack heating for vanishingly thin 00307 ! snow cover: add energy to snow/ground heat 00308 ! flux (W m-2) 00309 ! ZLWUPSNOW = upwelling longwave raaditive flux (W m-2) 00310 ! 00311 REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2_IC, ZTA_IC, ZQA_IC, 00312 ZPET_A_COEF_T, ZPEQ_A_COEF_T, ZPET_B_COEF_T, ZPEQ_B_COEF_T 00313 ! ZUSTAR2_IC = implicit lowest atmospheric level friction (m2/s2) 00314 ! ZTA_IC = implicit lowest atmospheric level air temperature 00315 ! ZQA_IC = implicit lowest atmospheric level specific humidity 00316 ! ZPET_A_COEF_T = transformed A-air temperature coefficient 00317 ! ZPET_B_COEF_T = transformed B-air temperature coefficient 00318 ! ZPEQ_A_COEF_T = transformed A-air specific humidity coefficient 00319 ! ZPEQ_B_COEF_T = transformed B-air specific humidity coefficient 00320 ! 00321 REAL, PARAMETER :: ZSNOWDZMIN = 0.0001 ! (m) 00322 ! ZSNOWDZMIN = minimum snow layer thickness for 00323 ! thermal calculations. Used to prevent 00324 ! numerical problems as snow becomes 00325 ! vanishingly thin. 00326 ! 00327 INTEGER :: JJ,JST ! looping indexes 00328 ! 00329 LOGICAL :: OCOND_GRAIN, OCOND_YEN 00330 REAL, DIMENSION(SIZE(PTA)) :: ZSNOWRHOF, ZSNOWDZF,ZSNOWGRAN1F, 00331 ZSNOWGRAN2F, ZSNOWHISTF 00332 REAL, DIMENSION(SIZE(PTA)) :: ZSNOWAGEF 00333 00334 ! New roughness lengths in case of glaciers without snow. 00335 REAL, DIMENSION(SIZE(PTA))::ZZ0_SNOWICE,ZZ0H_SNOWICE,ZZ0EFF_SNOWICE 00336 00337 !spectral albedo (3 bands for now) :: ready to output if necessary 00338 REAL, DIMENSION(SIZE(PSNOWRHO,1),3):: ZSPECTRALALBEDO 00339 00340 INTEGER, DIMENSION(SIZE(PTA)) :: INLVLS_USE ! varying number of effective layers 00341 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZWHOLDMAX 00342 INTEGER :: IPRINT ! gridpoint number to be printed 00343 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00344 00345 REAL :: ZTSTEPDAYS ! time step in days 00346 00347 LOGICAL::GCROINFOPRINT ! print daily informations 00348 LOGICAL::GCRODEBUGPRINT,GCRODEBUGDETAILSPRINT,GCRODEBUGPRINTATM ! print diagnostics for debugging 00349 LOGICAL::GCRODEBUGPRINTBALANCE 00350 INTEGER::IDEBUG 00351 00352 ! To stop simulation if mass/energy balances not closed 00353 LOGICAL,PARAMETER::LPSTOPBALANCE=.FALSE. 00354 00355 !To control and print eneregy balance 00356 REAL , DIMENSION(SIZE(PTA)) :: ZSUMMASS_INI,ZSUMHEAT_INI,ZSUMMASS_FIN,ZSUMHEAT_FIN 00357 ! 00358 IF (LHOOK) CALL DR_HOOK('SNOWCRO',0,ZHOOK_HANDLE) 00359 ! 00360 ! Look if we have to print snow profiles for debugging 00361 GCROINFOPRINT = LCRODAILYINFO .AND. (TPTIME%TIME ==0.0) 00362 ! 00363 IF (LCRODEBUG) THEN 00364 GCRODEBUGPRINT = (TPTIME%TDATE%YEAR*10000+TPTIME%TDATE%MONTH*100+TPTIME%TDATE%DAY & 00365 >= NTIMECRODEBUG) .AND. (TPTIME%TIME/3600.>=NHOURCRODEBUG) .AND. & 00366 (TPTIME%TDATE%YEAR*10000+TPTIME%TDATE%MONTH*100+TPTIME%TDATE%DAY & 00367 < NENDCRODEBUG) 00368 GCRODEBUGDETAILSPRINT = LCRODEBUGDETAILS.AND.GCRODEBUGPRINT 00369 GCRODEBUGPRINTATM = LCRODEBUGATM.AND.GCRODEBUGPRINT 00370 ELSE 00371 GCRODEBUGPRINT = .FALSE. 00372 GCRODEBUGDETAILSPRINT = .FALSE. 00373 GCRODEBUGPRINTATM = .FALSE. 00374 END IF 00375 ! 00376 ! Look if we have to compute and print energy balance control 00377 GCRODEBUGPRINTBALANCE = LCONTROLBALANCE .AND. & 00378 (TPTIME%TDATE%YEAR*10000+TPTIME%TDATE%MONTH*100+TPTIME%TDATE%DAY & 00379 >= NTIMECRODEBUG) .AND. (TPTIME%TIME/3600.>=NHOURCRODEBUG) .AND. & 00380 (TPTIME%TDATE%YEAR*10000+TPTIME%TDATE%MONTH*100+TPTIME%TDATE%DAY & 00381 < NENDCRODEBUG) 00382 ! 00383 IF (LCRODEBUG .OR. GCROINFOPRINT .OR. GCRODEBUGPRINTBALANCE) THEN 00384 IF (XLATCRODEBUG >= -90 .AND. XLONCRODEBUG >= -180.) THEN 00385 CALL GETPOINT_CRODEBUG(PXLAT,PXLON,IDEBUG) 00386 ELSE 00387 IDEBUG=NPOINTCRODEBUG 00388 END IF 00389 END IF 00390 ! 00391 ZUSTAR2_IC = 0.0 00392 ZTA_IC = 0.0 00393 ZQA_IC = 0.0 00394 ! 00395 OCOND_GRAIN=.TRUE. 00396 OCOND_YEN=.TRUE.!FALSE. !(if TRUE : use of the Yen (1981) thermal conductivity paramztrization ; 00397 ! otherwise, use the default ISBA-ES thermal conductivity parametrization) 00398 ! 00399 PGRNDFLUX=0. 00400 PSNOWHMASS=0. 00401 PHSNOW=0. 00402 PRNSNOW=0. 00403 PLES3L=0. 00404 PLEL3L=0. 00405 PHPSNOW=0. 00406 PEVAPCOR=0. 00407 PTHRUFAL=0. 00408 ! 00409 ! pour imprimer des diagnostics sur un des points 00410 IPRINT=1 00411 ! 00412 ! - - --------------------------------------------------- 00413 ! 00414 ! 0. Initialization 00415 ! -------------- 00416 ! NOTE that snow layer thickness is used throughout this code: SWE 00417 ! is only used to diagnose the thickness at the beginning of this routine 00418 ! and it is updated at the end of this routine. 00419 ! 00420 ! Initialization of the actual number of snow layers, total snow depth 00421 ! and layer thicknesses 00422 ! 00423 GSNOWFALL(:) = .FALSE. 00424 INLVLS_USE(:)=0 00425 DO JJ=1, SIZE(ZSNOW) 00426 DO JST=1,SIZE(PSNOWSWE(:,:),2) 00427 IF (PSNOWSWE(JJ,JST)>0.) THEN 00428 PSNOWDZ(JJ,JST) = PSNOWSWE(JJ,JST)/PSNOWRHO(JJ,JST) 00429 INLVLS_USE(JJ)= JST 00430 ELSE 00431 PSNOWDZ(JJ,JST) = 0. 00432 ENDIF 00433 ENDDO ! end loop snow layers 00434 ENDDO ! end loop grid points 00435 00436 ! Incrementation of snow layers age 00437 ZTSTEPDAYS=PTSTEP/86400. ! time step in days 00438 WHERE (PSNOWSWE >0) 00439 PSNOWAGE=PSNOWAGE+ZTSTEPDAYS 00440 ENDWHERE 00441 ! 00442 !Compute total SWE and heat for energy control 00443 IF (GCRODEBUGPRINTBALANCE) THEN 00444 DO JJ=1, SIZE(ZSNOW) 00445 ZSUMMASS_INI(JJ) = SUM(PSNOWSWE (JJ,1:INLVLS_USE(JJ))) 00446 ZSUMHEAT_INI(JJ) = SUM(PSNOWHEAT(JJ,1:INLVLS_USE(JJ))) 00447 ENDDO ! end loop grid points 00448 ENDIF 00449 ! 00450 ! Print of some simulation characteristics 00451 ! 00452 IF(GCROINFOPRINT) THEN 00453 CALL SNOWCROPRINTDATE() 00454 WRITE(*,FMT="(A12,I3,A12,I4)") & 00455 'nlayer:',INLVLS_USE(IDEBUG), ' nbpoints:', SIZE(ZSNOW) 00456 ! WRITE(*,*) 'PZ0H: ', PZ0H(IDEBUG) 00457 write(*,*) 'Snow fraction =',PPSN3L(IDEBUG) 00458 ENDIF 00459 ! 00460 IF (GCRODEBUGPRINT) THEN 00461 CALL SNOWCROPRINTDATE() 00462 CALL SNOWCROPRINTPROFILE("crocus initialization",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00463 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),PSNOWTEMP(IDEBUG,:),& 00464 PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),PSNOWGRAN1(IDEBUG,:),& 00465 PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00466 END IF 00467 00468 IF (GCRODEBUGPRINTATM) THEN 00469 CALL SNOWCROPRINTATM("forcing data :",PTA(IDEBUG),PQA(IDEBUG),PVMOD(IDEBUG),& 00470 PRR(IDEBUG),PSR(IDEBUG),PSW_RAD(IDEBUG),PLW_RAD(IDEBUG),& 00471 PTG(IDEBUG),PSOILCOND(IDEBUG),PD_G(IDEBUG),PPSN3L(IDEBUG)) 00472 END IF 00473 ! 00474 !* 1. Snow total depth 00475 ! ---------------- 00476 ! 00477 ZSNOW(:) = 0. 00478 DO JJ=1, SIZE(ZSNOW) 00479 ZSNOW(JJ)=SUM(PSNOWDZ(JJ,1:INLVLS_USE(JJ))) 00480 ENDDO 00481 ! 00482 ZSNOWBIS(:) = ZSNOW(:) 00483 ! 00484 ! 00485 !* 2. Snowfall 00486 ! -------- 00487 ! Calculate uppermost density and thickness changes due to snowfall, 00488 ! and add heat content of falling snow 00489 ! 00490 ! 00491 CALL SNOWNLFALL_UPGRID(TPTIME, OGLACIER, & 00492 PTSTEP,PSR,PTA,PVMOD,ZSNOWBIS,PSNOWRHO,PSNOWDZ, & 00493 PSNOWHEAT,PSNOWHMASS,PSNOWALB,PPERMSNOWFRAC, & 00494 PSNOWGRAN1,PSNOWGRAN2,GSNOWFALL, ZSNOWDZN,& 00495 ZSNOWRHOF, ZSNOWDZF,ZSNOWGRAN1F, ZSNOWGRAN2F, ZSNOWHISTF, & 00496 ZSNOWAGEF, GMODIF_MAILLAGE,INLVLS_USE,LSNOWDRIFT) 00497 ! 00498 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWFALL_UPGRID",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00499 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),PSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00500 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00501 ! 00502 ZSNOW(:) = ZSNOWBIS(:) 00503 ! 00504 !* 3. Update grid/discretization 00505 ! -------------------------- 00506 ! Reset grid to conform to model specifications: 00507 ! 00508 DO JJ=1, SIZE(ZSNOW) 00509 IF(GMODIF_MAILLAGE(JJ)) THEN 00510 ! 00511 CALL SNOWNLGRIDFRESH_1D(ZSNOW(JJ),PSNOWDZ(JJ,:),ZSNOWDZN(JJ,:), & 00512 PSNOWRHO(JJ,:),PSNOWHEAT(JJ,:),PSNOWGRAN1(JJ,:),PSNOWGRAN2(JJ,:), & 00513 PSNOWHIST(JJ,:),PSNOWAGE(JJ,:),GSNOWFALL(JJ), & 00514 ZSNOWRHOF(JJ),ZSNOWDZF(JJ),PSNOWHMASS(JJ),ZSNOWGRAN1F(JJ),& 00515 ZSNOWGRAN2F(JJ), ZSNOWHISTF(JJ),ZSNOWAGEF(JJ), & 00516 INLVLS_USE(JJ)) 00517 ENDIF 00518 ENDDO 00519 ! 00520 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWNLGRIDFRESH_1D",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00521 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),PSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00522 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00523 ! 00524 !* 4. Liquid water content and snow temperature 00525 ! ----------------------------------------- 00526 ! 00527 ! First diagnose snow temperatures and liquid 00528 ! water portion of the snow from snow heat content: 00529 ! update some snow layers parameters after new discretization 00530 ! 00531 DO JJ=1, SIZE(ZSNOW) 00532 ! active layers 00533 DO JST=1,INLVLS_USE(JJ) 00534 PSNOWSWE(JJ,JST) = PSNOWDZ(JJ,JST)*PSNOWRHO(JJ,JST) 00535 ! 00536 ZSCAP(JJ,JST) = SNOW3LSCAP(PSNOWRHO(JJ,JST)) 00537 ! 00538 ZSNOWTEMP(JJ,JST) = XTT + ( ((PSNOWHEAT(JJ,JST)/PSNOWDZ(JJ,JST)) & 00539 + XLMTT*PSNOWRHO(JJ,JST))/ZSCAP(JJ,JST) ) 00540 ! 00541 PSNOWLIQ(JJ,JST) = MAX(0.0,ZSNOWTEMP(JJ,JST)-XTT)*ZSCAP(JJ,JST)* & 00542 PSNOWDZ(JJ,JST)/(XLMTT*XRHOLW) 00543 ! 00544 ZSNOWTEMP(JJ,JST) = MIN(XTT,ZSNOWTEMP(JJ,JST)) 00545 ENDDO ! end loop active snow layers 00546 ! unactive layers 00547 DO JST=INLVLS_USE(JJ)+1,SIZE(PSNOWSWE,2) 00548 PSNOWSWE(JJ,JST)=0.0 00549 PSNOWRHO(JJ,JST)=999. 00550 PSNOWDZ(JJ,JST)=0. 00551 PSNOWGRAN1(JJ,JST)=0. 00552 PSNOWGRAN2(JJ,JST)=0. 00553 PSNOWHIST(JJ,JST)=0. 00554 PSNOWAGE(JJ,JST)=0. 00555 PSNOWHEAT(JJ,JST)=0. 00556 ZSNOWTEMP(JJ,JST)=XTT 00557 PSNOWLIQ(JJ,JST)=0. 00558 ENDDO ! end loop unactive snow layers 00559 ENDDO ! end loop grid points 00560 ! 00561 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE( 00562 "after liquid water/& &temperature diagnostic",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00563 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00564 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00565 ! 00566 ! 4.BIS Snow metamorphism 00567 ! ----------------- 00568 ! 00569 CALL SNOWCROMETAMO(PSNOWDZ,PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,ZSNOWTEMP, & 00570 PSNOWLIQ,PTSTEP,PSNOWSWE,INLVLS_USE) 00571 ! 00572 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROMETAMO",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00573 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00574 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00575 00576 !* 5. Snow Compaction 00577 ! --------------- 00578 ! Calculate snow density: compaction/aging: density increases 00579 ! 00580 CALL SNOWCROCOMPACTN(PTSTEP,PSNOWRHO,PSNOWDZ,ZSNOWTEMP,ZSNOW, & 00581 PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST, PSNOWLIQ,INLVLS_USE,& 00582 PDIRCOSZW) 00583 00584 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROCOMPACTN",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00585 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00586 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00587 00588 ! 00589 !* 5.1 Snow Compaction and Metamorphism due to snow drift 00590 ! --------------- 00591 ! PRINT*,LSNOWDRIFT 00592 IF (LSNOWDRIFT) THEN 00593 CALL SNOWDRIFT(PTSTEP, PVMOD, PSNOWRHO,PSNOWDZ, ZSNOW, & 00594 PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,INLVLS_USE,PTA,PQA,PPS,PRHOA) 00595 00596 ! old call before EB modifications for sublimation (to remove) 00597 ! CALL SNOWDRIFT(PTSTEP, PVMOD, PSNOWRHO,PSNOWDZ, ZSNOW, & 00598 ! PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,INLVLS_USE) 00599 00600 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWDRIFT",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00601 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00602 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00603 END IF 00604 ! 00605 ! Update snow heat content (J/m2) using dry density instead of total density: 00606 ! 00607 DO JJ=1, SIZE(ZSNOW) 00608 DO JST=1,INLVLS_USE(JJ) 00609 ZSCAP(JJ,JST) = SNOW3LSCAP(PSNOWRHO(JJ,JST)-PSNOWLIQ(JJ,JST)* & 00610 XRHOLW/MAX(PSNOWDZ(JJ,JST),ZSNOWDZMIN)) 00611 PSNOWHEAT(JJ,JST) = PSNOWDZ(JJ,JST)*( ZSCAP(JJ,JST)*(ZSNOWTEMP(JJ,JST)-XTT)- & 00612 XLMTT*PSNOWRHO(JJ,JST) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST) 00613 ENDDO ! end loop snow layers 00614 ENDDO ! end loop grid points 00615 ! 00616 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after update snow heat content",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00617 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00618 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00619 ! 00620 !* 6. Solar radiation transmission 00621 ! ----------------------------- 00622 ! 00623 ! Heat source (-sink) term due to shortwave 00624 ! radiation transmission within the snowpack: 00625 ! 00626 CALL SNOWCRORAD(TPTIME,OGLACIER,& 00627 PSW_RAD,PSNOWALB,PSNOWDZ,PSNOWRHO, & 00628 PALB,ZRADSINK,ZRADXS, & 00629 PSNOWGRAN1, PSNOWGRAN2, PSNOWAGE,PPS, & 00630 PZENITH, PPERMSNOWFRAC,INLVLS_USE ) 00631 ! 00632 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCRORAD",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00633 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00634 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00635 ! 00636 !* 7. Heat transfer and surface energy budget 00637 ! --------------------------------------- 00638 ! Snow thermal conductivity: 00639 ! 00640 CALL SNOWCROTHRM(PSNOWRHO,ZSCOND,ZSNOWTEMP,PPS,PSNOWLIQ,OCOND_GRAIN,OCOND_YEN) 00641 ! 00642 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROTHRM",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00643 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00644 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00645 ! 00646 ! Precipitation heating term: 00647 ! Rainfall renders it's heat to the snow when it enters 00648 ! the snowpack: 00649 ! 00650 PHPSNOW(:) = PRR(:)*XCL*(MAX(XTT,PTA(:))-XTT) ! (W/m2) 00651 ! 00652 ! Surface Energy Budget calculations using ISBA linearized form 00653 ! and standard ISBA turbulent transfer formulation 00654 ! 00655 00656 IF (ALL(PPEW_A_COEF==0.)) THEN 00657 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00658 ! Modif Matthieu Lafaysse for glaciers 00659 ! For surface ice, modify roughness lengths 00660 ! Only if not implicit coupling 00661 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00662 WHERE(PSNOWRHO(:,1)>XRHOTHRESHOLD_ICE) 00663 ZZ0_SNOWICE = PZ0*XZ0ICEZ0SNOW 00664 ZZ0H_SNOWICE = PZ0H*XZ0ICEZ0SNOW 00665 ZZ0EFF_SNOWICE = PZ0EFF*XZ0ICEZ0SNOW 00666 ELSEWHERE 00667 ZZ0_SNOWICE = PZ0 00668 ZZ0H_SNOWICE = PZ0H 00669 ZZ0EFF_SNOWICE = PZ0EFF 00670 ENDWHERE 00671 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00672 ELSE 00673 ZZ0_SNOWICE = PZ0 00674 ZZ0H_SNOWICE = PZ0H 00675 ZZ0EFF_SNOWICE = PZ0EFF 00676 END IF 00677 00678 CALL SNOWCROEBUD(HSNOWRES, HIMPLICIT_WIND, & 00679 PPEW_A_COEF, PPEW_B_COEF, & 00680 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, & 00681 ZSNOWDZMIN, & 00682 PZREF,ZSNOWTEMP(:,1),PSNOWRHO(:,1),PSNOWLIQ(:,1),ZSCAP(:,1), & 00683 ZSCOND(:,1),ZSCOND(:,2), & 00684 PUREF,PEXNS,PEXNA,PDIRCOSZW,PVMOD, & 00685 PLW_RAD,PSW_RAD,PTA,PQA,PPS,PTSTEP, & 00686 PSNOWDZ(:,1),PSNOWDZ(:,2),PSNOWALB,ZZ0_SNOWICE, & 00687 ZZ0EFF_SNOWICE,ZZ0H_SNOWICE, & 00688 ZSFCFRZ,ZRADSINK(:,1),PHPSNOW, & 00689 ZCT,PEMISNOW,PRHOA,ZTSTERM1,ZTSTERM2,ZRA,PCDSNOW,PCHSNOW, & 00690 ZQSAT, ZDQSAT, ZRSRA, ZUSTAR2_IC, PRI, & 00691 ZPET_A_COEF_T,ZPEQ_A_COEF_T,ZPET_B_COEF_T,ZPEQ_B_COEF_T ) 00692 ! 00693 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROEBUD",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00694 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00695 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00696 ! 00697 ! Heat transfer: simple diffusion along the thermal gradient 00698 ! 00699 ZSNOWTEMPO1(:) = ZSNOWTEMP(:,1) ! save surface snow temperature before update 00700 ! 00701 CALL SNOWCROSOLVT(PTSTEP,ZSNOWDZMIN,PSNOWDZ,ZSCOND,ZSCAP,PTG, & 00702 PSOILCOND,PD_G,ZRADSINK,ZCT,ZTSTERM1,ZTSTERM2, & 00703 ZPET_A_COEF_T,ZPEQ_A_COEF_T,ZPET_B_COEF_T,ZPEQ_B_COEF_T, & 00704 ZTA_IC,ZQA_IC,PGRNDFLUX, ZSNOWTEMP ,ZSNOWFLUX, & 00705 INLVLS_USE ) 00706 ! 00707 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROSOLVT",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00708 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00709 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00710 ! 00711 !* 8. Surface fluxes 00712 ! -------------- 00713 ! 00714 CALL SNOWCROFLUX(ZSNOWTEMP(:,1),PSNOWDZ(:,1),PEXNS,PEXNA, & 00715 ZUSTAR2_IC, & 00716 PTSTEP,PSNOWALB,PSW_RAD,PEMISNOW,ZLWUPSNOW,PLW_RAD, & 00717 ZTA_IC,ZSFCFRZ,ZQA_IC,PHPSNOW, & 00718 ZSNOWTEMPO1,ZSNOWFLUX,ZCT,ZRADSINK(:,1), & 00719 ZQSAT,ZDQSAT,ZRSRA, & 00720 PRNSNOW,PHSNOW,PGFLUXSNOW,PLES3L,PLEL3L,PEVAP, & 00721 PUSTAR ) 00722 ! 00723 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROFLUX",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00724 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00725 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00726 !* 9. Snow melt 00727 ! --------- 00728 ! 00729 ! First Test to see if snow pack vanishes during this time step: 00730 ! 00731 CALL SNOWCROGONE(PTSTEP,PLEL3L,PLES3L,PSNOWRHO, & 00732 PSNOWHEAT,ZRADSINK,PEVAPCOR,PTHRUFAL,PGRNDFLUX, & 00733 PGFLUXSNOW,PSNOWDZ,PSNOWLIQ,ZSNOWTEMP,ZRADXS, & 00734 PRR,INLVLS_USE) 00735 ! 00736 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROGONE",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00737 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00738 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00739 ! 00740 ! Add radiation not absorbed by snow to soil/vegetation interface flux 00741 ! (for thin snowpacks): 00742 ! 00743 PGRNDFLUX(:) = PGRNDFLUX(:) + ZRADXS(:) 00744 ! 00745 ! Second Test to see if one or several snow layers vanishe during this time 00746 ! step. In such a case, the concerned snow layers are agregated to neighbours 00747 00748 CALL SNOWCROLAYER_GONE(PTSTEP,ZSCAP,ZSNOWTEMP,PSNOWDZ, & 00749 PSNOWRHO,PSNOWLIQ,PSNOWGRAN1,PSNOWGRAN2, & 00750 PSNOWHIST,PSNOWAGE,PLES3L, INLVLS_USE) 00751 ! 00752 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROLAYER_GONE",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00753 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00754 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00755 ! 00756 ! 00757 ! For partial melt: transform excess heat content into snow liquid: 00758 ! 00759 CALL SNOWCROMELT(ZSCAP,ZSNOWTEMP,PSNOWDZ,PSNOWRHO,PSNOWLIQ,INLVLS_USE) 00760 ! 00761 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROMELT",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00762 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00763 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:))! 00764 ! 00765 !* 10. Snow water flow and refreezing 00766 ! ------------------------------ 00767 ! Liquid water vertical transfer and possible snowpack runoff 00768 ! And refreezing/freezing of meltwater/rainfall (ripening of the snow) 00769 ! 00770 CALL SNOWCROREFRZ(PTSTEP,PRR,PSNOWRHO,ZSNOWTEMP,PSNOWDZ,PSNOWLIQ,PTHRUFAL, & 00771 ZSCAP,PLEL3L,INLVLS_USE) 00772 ! 00773 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROREFRZ",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00774 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00775 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00776 ! 00777 !* 11. Snow Evaporation/Sublimation mass updates: 00778 ! ------------------------------------------ 00779 ! 00780 CALL SNOWCROEVAPN(PLES3L,PTSTEP,ZSNOWTEMP(:,1),PSNOWRHO(:,1), & 00781 PSNOWDZ(:,1),PEVAPCOR,PSNOWHMASS) 00782 ! 00783 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROEVAPN",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00784 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00785 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00786 ! 00787 ! If all snow in uppermost layer evaporates/sublimates, re-distribute 00788 ! grid (below could be evoked for vanishingly thin snowpacks): 00789 ! 00790 CALL SNOWCROEVAPGONE(PSNOWHEAT,PSNOWDZ,PSNOWRHO,ZSNOWTEMP,PSNOWLIQ,PSNOWGRAN1,& 00791 PSNOWGRAN2,PSNOWHIST,PSNOWAGE,INLVLS_USE) 00792 ! 00793 IF (GCRODEBUGDETAILSPRINT) CALL SNOWCROPRINTPROFILE("after SNOWCROEVAPGONE",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00794 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),ZSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),PSNOWHEAT(IDEBUG,:),& 00795 PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00796 ! 00797 !* 12. Update surface albedo: 00798 ! ---------------------- 00799 ! Snow clear sky albedo: 00800 ! 00801 CALL SNOWCROALB(TPTIME,OGLACIER,& 00802 PSNOWALB,ZSPECTRALALBEDO,PSNOWDZ(:,1),PSNOWRHO(:,1:2), & 00803 PPERMSNOWFRAC,PSNOWGRAN1(:,1),PSNOWGRAN2(:,1), & 00804 PSNOWAGE(:,1),PSNOWGRAN1(:,2),PSNOWGRAN2(:,2),PSNOWAGE(:,2),& 00805 PPS, PZENITH, INLVLS_USE) 00806 ! 00807 !* 13. Update snow heat content: 00808 ! ------------------------- 00809 ! Update the heat content (variable stored each time step) 00810 ! using current snow temperature and liquid water content: 00811 ! 00812 ! First, make check to make sure heat content not too large 00813 ! (this can result due to signifigant heating of thin snowpacks): 00814 ! add any excess heat to ground flux: 00815 ! 00816 DO JJ=1, SIZE(ZSNOW) 00817 ! active layers 00818 DO JST=1,INLVLS_USE(JJ) 00819 ZWHOLDMAX(JJ,JST) = SNOWCROHOLD(PSNOWRHO(JJ,JST),PSNOWLIQ(JJ,JST),PSNOWDZ(JJ,JST)) 00820 ZLIQHEATXS(JJ) = MAX(0.0, PSNOWLIQ(JJ,JST)*XRHOLW - & 00821 ZWHOLDMAX(JJ,JST)*XRHOLW) *XLMTT/PTSTEP 00822 PSNOWLIQ(JJ,JST) = PSNOWLIQ(JJ,JST) -ZLIQHEATXS(JJ)*PTSTEP/(XRHOLW*XLMTT) 00823 PSNOWLIQ(JJ,JST) = MAX(0.0, PSNOWLIQ(JJ,JST)) 00824 PGRNDFLUX(JJ) = PGRNDFLUX(JJ) + ZLIQHEATXS(JJ) 00825 PSNOWTEMP(JJ,JST) = ZSNOWTEMP(JJ,JST) 00826 ! Heat content using total density 00827 ZSCAP(JJ,JST) = SNOW3LSCAP(PSNOWRHO(JJ,JST)) 00828 PSNOWHEAT(JJ,JST) = PSNOWDZ(JJ,JST)*( ZSCAP(JJ,JST)*(PSNOWTEMP(JJ,JST)-XTT)& 00829 - XLMTT*PSNOWRHO(JJ,JST)) +XLMTT*XRHOLW*PSNOWLIQ(JJ,JST) 00830 ! 00831 PSNOWSWE(JJ,JST) = PSNOWDZ(JJ,JST)*PSNOWRHO(JJ,JST) 00832 ENDDO ! end loop active snow layers 00833 ! 00834 ! unactive layers 00835 DO JST=INLVLS_USE(JJ)+1,SIZE(PSNOWSWE,2) 00836 PSNOWSWE(JJ,JST)=0. 00837 PSNOWHEAT(JJ,JST) =0. 00838 PSNOWRHO(JJ,JST)=999. 00839 PSNOWTEMP(JJ,JST)=0. 00840 PSNOWDZ(JJ,JST)=0. 00841 ENDDO ! end loop unactive snow layers 00842 ! 00843 ENDDO ! end loop grid points 00844 ! 00845 ! print some final diagnostics 00846 ! ! ! IF (INLVLS_USE(I)>0) THEN 00847 ! ! ! write(*,FMT="(I4,2I4,F4.0,A7,F5.2,A10,F7.1,A11,F6.2,A13,F6.2)") & 00848 ! ! ! TPTIME%TDATE%YEAR,TPTIME%TDATE%MONTH,TPTIME%TDATE%DAY,TPTIME%TIME/3600.,& 00849 ! ! ! 'HTN= ',SUM(PSNOWDZ(I,1:INLVLS_USE(I))), 'FLUX Sol=', PGRNDFLUX(I),& 00850 ! ! ! 'Tsurf_sol=',PTG(I)-273.16, 'Tbase_neige=',PSNOWTEMP(I,INLVLS_USE(I))-273.16 00851 ! ! ! ENDIF 00852 00853 IF (GCRODEBUGPRINT) THEN 00854 CALL SNOWCROPRINTDATE() 00855 CALL SNOWCROPRINTPROFILE("CROCUS : end of time step",INLVLS_USE(IDEBUG),LPRINTGRAN,& 00856 PSNOWDZ(IDEBUG,:),PSNOWRHO(IDEBUG,:),PSNOWTEMP(IDEBUG,:),PSNOWLIQ(IDEBUG,:),& 00857 PSNOWHEAT(IDEBUG,:),PSNOWGRAN1(IDEBUG,:),PSNOWGRAN2(IDEBUG,:),& 00858 PSNOWHIST(IDEBUG,:),PSNOWAGE(IDEBUG,:)) 00859 END IF 00860 00861 ! check suspect low temperature 00862 DO JJ=1,SIZE(ZSNOW) 00863 !IF(INLVLS_USE(JJ)>0) write(*,*) 'SOL:',JJ,INLVLS_USE(JJ),PGRNDFLUX(JJ),PTG(JJ),& 00864 ! PSNOWTEMP(jj,INLVLS_USE(JJ)),PSNOWTEMP(jj,1),PZENITH(JJ) 00865 DO JST=1,INLVLS_USE(JJ) 00866 IF (PSNOWTEMP(JJ,JST) < 100.) THEN 00867 write(*,*) 'pb tempe snow :',PSNOWTEMP(JJ,JST) 00868 write(*,FMT='("DATE:",2(I2.2,"/"),I4.4,F3.0)') & 00869 TPTIME%TDATE%DAY,TPTIME%TDATE%MONTH,TPTIME%TDATE%YEAR,TPTIME%TIME/3600. 00870 write(*,*) 'point',JJ,"/",SIZE(ZSNOW) 00871 write(*,*) 'layer',JST 00872 write(*,*) 'pressure',PPS(JJ) 00873 write(*,*) 'slope',ACOS(PDIRCOSZW(JJ))*(180./XPI),"°" 00874 write(*,*) 'XLAT=',PXLAT(JJ),'XLON=',PXLON(JJ) 00875 write(*,*) 'solar radiation=',PSW_RAD(JJ) 00876 write(*,*) 'INLVLS_USE(JJ):',INLVLS_USE(JJ) 00877 write(*,*) PSNOWDZ(JJ,1:INLVLS_USE(JJ)) 00878 write(*,*) PSNOWRHO(JJ,1:INLVLS_USE(JJ)) 00879 write(*,*) PSNOWTEMP(JJ,1:INLVLS_USE(JJ)) 00880 CALL ABOR1_SFX('SNOWCRO: erreur tempe snow') 00881 ENDIF 00882 ENDDO 00883 ENDDO! 00884 00885 00886 !Control and print energy balance 00887 IF (GCRODEBUGPRINTBALANCE) THEN 00888 00889 ZSUMMASS_FIN(IDEBUG)=SUM(PSNOWSWE(IDEBUG,1:INLVLS_USE(JJ))) 00890 ZSUMHEAT_FIN(IDEBUG)=SUM(PSNOWHEAT(IDEBUG,1:INLVLS_USE(JJ))) 00891 00892 CALL SNOWCROPRINTBALANCE(ZSUMMASS_INI(IDEBUG),ZSUMHEAT_INI(IDEBUG),ZSUMMASS_FIN(IDEBUG),& 00893 ZSUMHEAT_FIN(IDEBUG),PSR(IDEBUG),PRR(IDEBUG),PTHRUFAL(IDEBUG),PEVAP(IDEBUG),& 00894 PEVAPCOR(IDEBUG),PGRNDFLUX(IDEBUG),PHSNOW(IDEBUG),PRNSNOW(IDEBUG),PLEL3L(IDEBUG),& 00895 PLES3L(IDEBUG),PHPSNOW(IDEBUG),PSNOWHMASS(IDEBUG),PSNOWDZ(IDEBUG,1),PTSTEP) 00896 00897 ENDIF 00898 00899 IF (LPSTOPBALANCE) THEN 00900 00901 ! bilan pour tous points pour test eventuel sur depassement seuil des residus 00902 DO JJ=1, SIZE(ZSNOW) 00903 ZSUMMASS_FIN(JJ)=SUM(PSNOWSWE(JJ,1:INLVLS_USE(JJ))) 00904 ZSUMHEAT_FIN(JJ)=SUM(PSNOWHEAT(JJ,1:INLVLS_USE(JJ))) 00905 ENDDO ! end loop grid points 00906 00907 CALL SNOWCROSTOPBALANCE(ZSUMMASS_INI(:),ZSUMHEAT_INI(:),ZSUMMASS_FIN(:),& 00908 ZSUMHEAT_FIN(:),PSR(:),PRR(:),PTHRUFAL(:),PEVAP(:),& 00909 PEVAPCOR(:),PGRNDFLUX(:),PHSNOW(:),PRNSNOW(:),PLEL3L(:),& 00910 PLES3L(:),PHPSNOW(:),PSNOWHMASS(:),PSNOWDZ(:,1),PTSTEP) 00911 END IF 00912 00913 ! 00914 IF (LHOOK) CALL DR_HOOK('SNOWCRO',1,ZHOOK_HANDLE) 00915 CONTAINS 00916 ! 00917 ! 00918 ! 00919 !#################################################################### 00920 !#################################################################### 00921 !#################################################################### 00922 SUBROUTINE SNOWCROCOMPACTN(PTSTEP,PSNOWRHO,PSNOWDZ, & 00923 PSNOWTEMP,PSNOW,PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST, & 00924 PSNOWLIQ,INLVLS_USE,PDIRCOSZW) 00925 ! 00926 !! PURPOSE 00927 !! ------- 00928 ! Snow compaction due to overburden and settling. 00929 ! Mass is unchanged: layer thickness is reduced 00930 ! in proportion to density increases. Method 00931 ! directly inherited from Crocus v2.4 and 00932 ! coarsely described in Brun et al., J. Glac 1989 and 1992 00933 ! 00934 ! de/e = -sigma/eta * dt 00935 ! 00936 ! where e is layer thickness, sigma is the vertical stress, dt is the 00937 ! time step and eta is the snow viscosity 00938 ! * sigma is currently calculated taking into account only the overburden 00939 ! (a term representing "metamorphism stress" in fresh snow may be added 00940 ! in the future) 00941 ! * eta is computed as a function of snowtype, density and temperature 00942 ! 00943 ! The local slope is taken into account, through the variable PDIRCOSZW 00944 ! which is directly the cosine of the local slope 00945 ! 00946 ! 00947 ! HISTORY: 00948 ! Basic structure from ISBA-ES model (Boone and Etchevers, 2001) 00949 ! Implementation of Crocus laws : E. Brun, S. Morin, J.-M. Willemet July 2010. 00950 ! Implementation of slope effect on settling : V. Vionnet, S. Morin May 2011 00951 ! 00952 ! 00953 USE MODD_CSTS, ONLY : XTT, XG 00954 USE MODD_SNOW_PAR, ONLY : XRHOSMAX_ES 00955 USE MODD_SNOW_METAMO 00956 ! 00957 IMPLICIT NONE 00958 ! 00959 !* 0.1 declarations of arguments 00960 ! 00961 REAL, INTENT(IN) :: PTSTEP ! Time step UNIT : s 00962 REAL, DIMENSION(:), INTENT(IN) :: PDIRCOSZW ! cosine of local slope 00963 ! 00964 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP ! Snow temperature UNIT : K 00965 ! 00966 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ ! Density UNIT : kg m-3, Layer thickness UNIT : m 00967 ! 00968 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW ! Snowheight UNIT : m 00969 ! 00970 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST, !Snowtype variables 00971 PSNOWLIQ ! Snow liquid water content UNIT ??? 00972 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE ! Number of snow layers used 00973 00974 ! 00975 !* 0.2 declarations of local variables 00976 ! 00977 INTEGER :: JJ,JST ! looping indexes 00978 ! 00979 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHO2, ! work snow density UNIT : kg m-3 00980 ZVISCOSITY, ! Snow viscosity UNIT : N m-2 s (= Pa s) 00981 ZSMASS !, & ! overburden mass for a given layer UNIT : kg m-2 00982 ! ZWSNOWDZ ! mass of each snow layer UNIT : kg m-2 00983 ! 00984 ! Compaction/Settling Coefficients from Crocus v2.4 00985 ! 00986 REAL, PARAMETER :: VVISC1=7.62237E6 ! pre-exponential viscosity factor (UNIT : N m-2 s) 00987 REAL, PARAMETER :: VVISC3=0.023 ! density adjustement in the exponential correction for viscosity (UNIT : m3 kg-1) 00988 REAL, PARAMETER :: VVISC4=.1 ! temperature adjustement in the exponential correction for viscosity (UNIT : K-1) 00989 REAL, PARAMETER :: VVISC5=1. ! factor for viscosity adjustement to grain type - to be checked 00990 REAL, PARAMETER :: VVISC6=60. ! factor for viscosity adjustement to grain type - to be checked 00991 ! (especially this one ; inconsistency with Crocus v2.4) 00992 REAL, PARAMETER :: VVISC7=10. ! factor for viscosity adjustement to grain type - to be checked 00993 REAL, PARAMETER :: VRO11=250. ! normalization term for density dependence of the viscosity calculation (UNIT : kg m-3) 00994 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00995 ! 00996 !------------------------------------------------------------------------------- 00997 ! 00998 ! 1. Cumulative snow mass (kg/m2): 00999 ! -------------------------------- 01000 ! 01001 01002 IF (LHOOK) CALL DR_HOOK('SNOWCROCOMPACTN',0,ZHOOK_HANDLE) 01003 DO JJ=1,SIZE(PSNOW) 01004 ZSMASS(JJ,1) = 0.0 01005 DO JST=2,INLVLS_USE(JJ) 01006 ZSMASS(JJ,JST) = ZSMASS(JJ,JST-1) + PSNOWDZ(JJ,JST-1) * PSNOWRHO(JJ,JST-1) 01007 ENDDO 01008 ZSMASS(JJ,1) = 0.5* PSNOWDZ(JJ,1) * PSNOWRHO(JJ,1) ! overburden of half the mass of the uppermost layer applied to itself 01009 ENDDO 01010 01011 ! 01012 ! 2. Compaction/Settling 01013 ! ---------------------- 01014 01015 DO JJ=1, SIZE(PSNOW) 01016 DO JST=1, INLVLS_USE(JJ) 01017 ! 01018 ! Snow viscosity basic equation (depend on temperature and density only): 01019 ! 01020 ZVISCOSITY(JJ,JST)=VVISC1 *EXP(VVISC3*PSNOWRHO(JJ,JST)+ & 01021 VVISC4*ABS(XTT-PSNOWTEMP(JJ,JST)))*PSNOWRHO(JJ,JST)/VRO11 01022 01023 ! Equations below apply changes to the basic viscosity value, based on snow microstructure properties 01024 01025 IF (PSNOWLIQ(JJ,JST)>0.0) THEN 01026 ZVISCOSITY(JJ,JST)=ZVISCOSITY(JJ,JST)/ & 01027 (VVISC5+VVISC6*PSNOWLIQ(JJ,JST)/PSNOWDZ(JJ,JST)) 01028 ENDIF 01029 01030 IF(PSNOWLIQ(JJ,JST)/PSNOWDZ(JJ,JST) <= 0.01 & 01031 .AND.PSNOWHIST(JJ,JST) >= NVHIS2) THEN 01032 ZVISCOSITY(JJ,JST)=ZVISCOSITY(JJ,JST)*VVISC7 01033 ENDIF 01034 01035 IF ((PSNOWGRAN1(JJ,JST)>=0..AND.PSNOWGRAN1(JJ,JST) < VGRAN6) & 01036 .AND. PSNOWHIST(JJ,JST) == NVHIS1) THEN 01037 ZVISCOSITY(JJ,JST)=ZVISCOSITY(JJ,JST)* & 01038 MIN(4.,EXP(MIN(VDIAM1,PSNOWGRAN2(JJ,JST)-VDIAM4)/VDIAM6)) 01039 ENDIF 01040 01041 ! 01042 ! Calculate new snow snow density: compaction from weight/over-burden 01043 ! 01044 ZSNOWRHO2(JJ,JST) = PSNOWRHO(JJ,JST) + PSNOWRHO(JJ,JST)*PTSTEP* & 01045 ((XG*PDIRCOSZW(JJ)*ZSMASS(JJ,JST)/ZVISCOSITY(JJ,JST))) 01046 ! 01047 ! Calculate new grid thickness in response to density change 01048 ! 01049 PSNOWDZ(JJ,JST) = PSNOWDZ(JJ,JST)*(PSNOWRHO(JJ,JST)/ZSNOWRHO2(JJ,JST)) 01050 ! 01051 ! Update density (kg m-3): 01052 ! 01053 PSNOWRHO(JJ,JST) = ZSNOWRHO2(JJ,JST) 01054 ! 01055 ! 01056 ENDDO ! end loop snow layers 01057 ENDDO ! end loop grid points 01058 ! 01059 ! 01060 ! 3. Update total snow depth: 01061 ! ----------------------------------------------- 01062 ! 01063 DO JJ=1, SIZE(PSNOWDZ,1) 01064 PSNOW(JJ)=SUM(PSNOWDZ(JJ,1:INLVLS_USE(JJ))) 01065 ENDDO 01066 IF (LHOOK) CALL DR_HOOK('SNOWCROCOMPACTN',1,ZHOOK_HANDLE) 01067 01068 ! 01069 !------------------------------------------------------------------------------- 01070 ! 01071 END SUBROUTINE SNOWCROCOMPACTN 01072 01073 01074 !#################################################################### 01075 !#################################################################### 01076 !#################################################################### 01077 SUBROUTINE SNOWCROMETAMO(PSNOWDZ,PSNOWGRAN1, PSNOWGRAN2, & 01078 PSNOWHIST, PSNOWTEMP, PSNOWLIQ, PTSTEP, & 01079 PSNOWSWE,INLVLS_USE) 01080 ! 01081 01082 !**** *METAMO* - METAMORPHOSE DES GRAINS 01083 ! - SNOW METAMORPHISM 01084 ! OBJET. 01085 ! ------ 01086 ! METAMORPHOSE DU MANTEAU NEIGEUX. 01087 ! EVOLUTION DU TYPE DE GRAINS 01088 ! MISE A JOUR DES VARIABLES HISTORIQUES. 01089 ! METAMORPHISM OF THE SNOW GRAINS, 01090 ! HISTORICAL VARIABLES 01091 01092 !** INTERFACE. 01093 ! ---------- 01094 ! FORMALISME ADOPTE POUR LA REPRESENTATION DES GRAINS : 01095 ! FORMALISM FOR THE REPRESENTATION OF GRAINS 01096 ! ----------------------------------------------------- 01097 01098 01099 ! 1 - -1 NEIGE FRAICHE 01100 ! / \ | ------------- 01101 ! / \ | DENDRICITE DECRITE EN TERME 01102 ! / \ | DENDRICITY DE DENDRICITE ET 01103 ! / \ | SPHERICITE 01104 ! 2---------3 - 0 DESCRIBED WITH 01105 ! SPHERICITY AND 01106 ! |---------| DENDRICITY 01107 ! 0 1 01108 ! SPHERICITE 01109 ! SPHERICITY 01110 01111 ! 4---------5 - 01112 ! | | | 01113 ! | | | DIAMETRE (OU TAILLE) 01114 ! | | | DIAMETER (OR SIZE ) 01115 ! | | | 01116 ! | | | NEIGE NON DENDRITIQUE 01117 ! 6---------7 - --------------------- 01118 01119 ! SPHERICITE ET TAILLE 01120 ! SPHERICITY AND SIZE 01121 01122 ! LES VARIABLES DU MODELE : 01123 ! ------------------------- 01124 ! CAS DENDRITIQUE CAS NON DENDRITIQUE 01125 ! 01126 ! SGRAN1(JST) : DENDRICITE SGRAN1(JST) : SPHERICITE 01127 ! SGRAN2(JST) : SPHERICITE SGRAN2(JST) : TAILLE (EN METRE) 01128 ! SIZE 01129 01130 ! 01131 ! CAS DENDRITIQUE/ DENDRITIC CASE 01132 ! ------------------------------- 01133 ! SGRAN1(JST) VARIE DE -VGRAN1 (-99 PAR DEFAUT) (ETOILE) A 0 01134 ! (DENDRICITY) >D OU LA DIVISION PAR -VGRAN1 POUR OBTENIR DES VALEURS 01135 ! ENTRE 1 ET 0 01136 ! VARIES FROM -VGRAN1 (DEFAULT -99) (FRESH SNOW) TO 0 01137 ! DIVISION BY -VGRAN1 TO OBTAIN VALUES BETWEEN 0 AND 1 01138 01139 ! SGRAN2(JST) VARIE DE 0 (CAS COMPLETEMENT ANGULEUX) A VGRAN1 01140 ! (SPHERICITY) (99 PAR DEFAUT) 01141 ! >D OU LA DIVISION PAR VGRAN1 POUR OBTENIR DES VALEURS 01142 ! ENTRE 0 ET 1 01143 ! VARIES FROM 0 (SPHERICITY=0) TO VGRAN1 01144 01145 01146 ! CAS NON DENDRITIQUE / NON DENDRITIC CASE 01147 ! --------------------------------------- 01148 01149 ! SGRAN1(JST) VARIE DE 0 (CAS COMPLETEMENT ANGULEUX) A VGRAN1 01150 ! (SPHERICITY) (99 PAR DEFAUT) (CAS SPHERIQUE) 01151 ! >D OU LA DIVISION PAR VGRAN1 POUR OBTENIR DES VALEURS 01152 ! ENTRE 0 ET 1 01153 ! VARIES FROM 0 TO 99 01154 01155 ! SGRAN2(JST) EST SUPERIEUR A VDIAM1-SPHERICITE (3.E-4 M) ET NE FAIT QUE CROITRE 01156 ! (SIZE) IS GREATER THAN VDIAM1-SPHERICITE (3.E-4 M) ALWAYS INCREASE 01157 01158 01159 ! EXEMPLES : POINTS CARACTERISTIQUES DE LA FIGURE 01160 ! -------- 01161 01162 ! SGRAN1 SGRAN2 DENDRICITE SPHERICITE TAILLE 01163 ! DENDRICITY SPHERICITY SIZE 01164 ! -------------------------------------------------------------- 01165 ! (M) 01166 ! 1 -VGRAN1 VNSPH3 1 0.5 01167 ! 2 0 0 0 0 01168 ! 3 0 VGRAN1 0 1 01169 ! 4 0 VDIAM1 0 4.E-4 01170 ! 5 VGRAN1 VDIAM1-VSPHE1 1 3.E-4 01171 ! 6 0 -- 0 -- 01172 ! 7 VGRAN1 -- 1 -- 01173 01174 ! PAR DEFAUT : VGRAN1 =99 VNSPH3=50 VSPHE1=1. VDIAM1=4.E-4 01175 01176 01177 ! METHODE. 01178 ! -------- 01179 ! EVOLUTION DES TYPES DE GRAINS : SELON LES LOIS DECRITES 01180 ! DANS BRUN ET AL (1992) 01181 ! PLUSIEURS CAS SONT A DISTINGUER 01182 ! 1.2 NEIGE HUMIDE 01183 ! 1.3 METAMORPHOSE NEIGE SECHE 01184 ! 1.3.1 FAIBLE GRADIENT 01185 ! 1.3.2 GRADIENT MOYEN 01186 ! 1.3.3 FORT GRADIENT 01187 ! DANS CHAQUE CAS ON SEPARE NEIGE DENDRITIQUE ET NON DENDRITIQUE 01188 ! LE PASSAGE DENDRITIQUE => NON DENDRITIQUE SE FAIT LORSQUE 01189 ! SGRAN1 DEVIENT > 0 01190 01191 ! TASSEMENT : LOIS DE VISCOSITE ADAPTEE SELON LE TYPE DE GRAINS 01192 01193 ! VARIABLES HISTORIQUES (CAS NON DENDRITIQUE SEULEMENT) 01194 01195 ! MSHIST DEFAUT 01196 ! 0 CAS NORMAL 01197 ! NVHIS1 1 GRAINS ANGULEUX 01198 ! NVHIS2 2 GRAINS AYANT ETE EN PRESENCE D EAU LIQUIDE 01199 ! MAIS N'AYANT PAS EU DE CARATERE ANGULEUX 01200 ! NVHIS3 3 GRAINS AYANT ETE EN PRESENCE D EAU LIQUIDE 01201 ! AYANT EU AUPARAVANT UN CARACTERE ANGULEUX 01202 01203 ! GRAIN METAMORPHISM ACCORDING TO BRUN ET AL (1992) 01204 ! THE DIFFERENT CASES ARE : 01205 ! 1.2 WET SNOW 01206 ! 1.3 DRY SNOW 01207 ! 1.3.1. LOW TEMPERATURE GRADIENT 01208 ! 1.3.2. MODERATE TEMPERATURE GRADIENT 01209 ! 1.3.3. HIGH TEMPERATURE GRADIENT 01210 ! THE CASE OF DENTRITIC OR NON DENDRITIC SNOW IS TREATED SEPARATELY 01211 ! THE LIMIT DENTRITIC ==> NON DENDRITIC IS REACHED WHEN SGRAN1>0 01212 01213 ! SNOW SETTLING : VISCOSITY DEPENDS ON THE GRAIN TYPES 01214 01215 ! HISTORICAL VARIABLES (NON DENDRITIC CASE) 01216 ! MSHIST DEFAUT 01217 ! 0 CAS NORMAL 01218 ! NVHIS1 1 FACETED CRISTAL 01219 ! NVHIS2 2 LIQUID WATER AND NO FACETED CRISTALS BEFORE 01220 ! NVHIS3 3 LIQUID WATER AND FACETED CRISTALS BEFORE 01221 01222 ! EXTERNES. 01223 ! --------- 01224 01225 ! REFERENCES. 01226 ! ----------- 01227 01228 ! AUTEURS. 01229 ! -------- 01230 ! ERIC BRUN ET AL. - JOURNAL OF GLACIOLOGY 1989/1992. 01231 01232 ! MODIFICATIONS. 01233 ! -------------- 01234 ! 08/95: YANNICK DANIELOU - CODAGE A LA NORME DOCTOR. 01235 ! 09/96: ERIC MARTIN - CORRECTION COMMENTAIRES 01236 ! 03/06: JM Willemet - F90 and SI units 01237 ! 08/06: JM Willemet - new formulation for TEL (Mwat/(Mice+Mwat) instead of Mwat/Mice. 01238 ! Threshold on the diameter increasing of the wet grains. 01239 ! 01/07 : JM Willemet - CORRECTION DES COUCHES SATUREES SUBISSANT DU TASSEMENT 01240 ! CORRECTION ON THE SATURATED LAYERS WHICH ARE SETTLED 01241 01242 USE MODD_SNOW_METAMO 01243 USE MODD_CSTS, ONLY : XTT, XPI,XRHOLW 01244 USE MODE_SNOW3L 01245 ! 01246 IMPLICIT NONE 01247 ! 01248 ! 0.1 declarations of arguments 01249 ! 01250 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ, PSNOWTEMP, 01251 PSNOWLIQ, PSNOWSWE 01252 ! 01253 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWGRAN1, PSNOWGRAN2, 01254 PSNOWHIST 01255 ! 01256 REAL, INTENT(IN) :: PTSTEP 01257 ! 01258 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 01259 ! 01260 ! 0.2 declaration of local variables 01261 ! 01262 REAL :: ZGRADT,ZTELM, ZVDENT, 01263 ZDENT, ZSPHE,ZVAP, 01264 ZDANGL 01265 REAL :: diam 01266 INTEGER :: INLVLS 01267 INTEGER JST,JJ !Loop controls 01268 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01269 01270 ! INITIALISATION 01271 ! -------------- 01272 01273 IF (LHOOK) CALL DR_HOOK('SNOWCROMETAMO',0,ZHOOK_HANDLE) 01274 INLVLS=SIZE(PSNOWGRAN1(:,:),2) ! total snow layers 01275 01276 !* 1. METAMORPHOSES DANS LES STRATES. 01277 ! METAMORPHISM 01278 ! ------------------------------- 01279 DO JJ=1,SIZE(PSNOWRHO,1) 01280 DO JST=1,INLVLS_USE(JJ) 01281 01282 ! 1.1 INITIALISATION: GRADIENT DE TEMPERATURE / TEMPERATURE GRADIENT 01283 01284 IF(JST == INLVLS_USE(JJ))THEN 01285 ZGRADT=ABS(PSNOWTEMP(JJ,INLVLS_USE(JJ))-PSNOWTEMP(JJ,INLVLS_USE(JJ)-1))*2. & 01286 /(PSNOWDZ(JJ,INLVLS_USE(JJ))+PSNOWDZ(JJ,INLVLS_USE(JJ)-1)) 01287 ELSEIF( JST == 1)THEN 01288 ZGRADT=ABS(PSNOWTEMP(JJ,2)-PSNOWTEMP(JJ,1))*2./(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 01289 ELSE 01290 ZGRADT=ABS(PSNOWTEMP(JJ,JST+1)-PSNOWTEMP(JJ,JST-1)) & 01291 /(PSNOWDZ(JJ,JST-1)*.5+PSNOWDZ(JJ,JST)+PSNOWDZ(JJ,JST+1)*.5) 01292 ENDIF 01293 01294 IF(PSNOWLIQ(JJ,JST) > UEPSI)THEN 01295 01296 01297 ! 1.2 METAMORPHOSE HUMIDE. 01298 ! WET SNOW METAMORPHISM 01299 01300 ! TENEUR EN EAU LIQUIDE / LIQUID WATER CONTENT 01301 01302 ZTELM=UPOURC*PSNOWLIQ(JJ,JST)*XRHOLW/PSNOWSWE(JJ,JST) 01303 01304 ! VITESSE DE DIMINUTION DE LA DENDRICITE. 01305 ! RATE OF THE DENDRICITY DECREASE 01306 ZVDENT=MAX(VDENT2*ZTELM**NVDENT1,VDENT1*EXP(VVAP1/XTT)) 01307 01308 IF(PSNOWGRAN1(JJ,JST) < -UEPSI)THEN 01309 01310 ! 1.2.1 CAS DENDRITIQUE/DENDRITIC CASE. 01311 01312 01313 ! VARIABLES DESCRIPTIVES DE LA DENDRICITE ET LA SPHERICITE. 01314 ZDENT=-PSNOWGRAN1(JJ,JST)/VGRAN1 01315 ZSPHE=PSNOWGRAN2(JJ,JST)/VGRAN1 01316 ! CALCUL NOUVELLE DENDRICITE ET SPHERICITE. 01317 ZDENT=ZDENT-ZVDENT*PTSTEP 01318 ZSPHE=ZSPHE+ZVDENT*PTSTEP 01319 IF(ZDENT <= UEPSI)THEN 01320 ! EVOLUTION DE SGRAN1 ET SGRAN2 ET TEST PASSAGE 01321 ! DENDRITIQUE > NON DENDRITIQUE. 01322 PSNOWGRAN1(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01323 PSNOWGRAN2(JJ,JST)=VDIAM1-VDIAM5*MIN(ZSPHE,VSPHE1) 01324 ELSE 01325 PSNOWGRAN1(JJ,JST)=-ZDENT*VGRAN1 01326 PSNOWGRAN2(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01327 ENDIF 01328 01329 ELSEIF(PSNOWGRAN1(JJ,JST) < VGRAN1-UEPSI)THEN 01330 01331 ! 1.2.2 CAS NON DENDRITIQUE NON COMPLETEMENT SPHERIQUE. 01332 ! EVOLUTION DE LA SPHERICITE SEULEMENT. 01333 ! NON DENDRITIC AND NOT COMPLETELY SPHERIC CASE 01334 ! EVOLUTION OF SPHERICITY ONLY (NOT SIZE) 01335 01336 ZSPHE=PSNOWGRAN1(JJ,JST)/VGRAN1 01337 ZSPHE=ZSPHE+ZVDENT*PTSTEP 01338 PSNOWGRAN1(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01339 ELSE 01340 01341 ! 1.2.3 CAS NON DENDRITIQUE ET SPHERIQUE/NON DENDRITIC AND SPHERIC. 01342 ! EVOLUTION DE LA TAILLE SEULEMENT/EVOLUTION OF SIZE ONLY 01343 01344 diam=PSNOWGRAN2(JJ,JST) 01345 PSNOWGRAN2(JJ,JST)=2. & 01346 *(3./(4.*XPI)*(4.*XPI/3.*(PSNOWGRAN2(JJ,JST)/2.)**3 & 01347 +(VTAIL1+VTAIL2*ZTELM**NVDENT1)*PTSTEP))**(1./3.) 01348 01349 ENDIF 01350 01351 ELSE 01352 01353 ! 1.3 METAMORPHOSES SECHES/DRY METAMORPHISM. 01354 01355 IF(ZGRADT < VGRAT1)THEN 01356 01357 ! 1.3.1 CALCUL METAMORPHOSE FAIBLE/LOW GRADIENT (0-5 DEG/M). 01358 01359 01360 ZVAP=EXP(VVAP1/PSNOWTEMP(JJ,JST)) 01361 IF(PSNOWGRAN1(JJ,JST) < -UEPSI)THEN 01362 01363 ! 1.3.1.1 CAS DENDRITIQUE / DENDRITIC CASE. 01364 01365 ZDENT=-PSNOWGRAN1(JJ,JST)/VGRAN1 01366 ZSPHE=PSNOWGRAN2(JJ,JST)/VGRAN1 01367 ZDENT=ZDENT-VDENT1*ZVAP*PTSTEP 01368 ZSPHE=ZSPHE+VSPHE2*ZVAP*PTSTEP 01369 IF(ZDENT < UEPSI)THEN 01370 ! EVOLUTION DE SGRAN1 ET SGRAN2 ET TEST PASSAGE 01371 ! DENDRITIQUE > NON DENDRITIQUE. 01372 PSNOWGRAN1(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01373 PSNOWGRAN2(JJ,JST)=VDIAM1-VDIAM5*MIN(VSPHE1,ZSPHE) 01374 ELSE 01375 PSNOWGRAN1(JJ,JST)=-ZDENT*VGRAN1 01376 PSNOWGRAN2(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01377 ENDIF 01378 ELSE 01379 01380 ! 1.3.1.2 CAS NON DENDRITIQUE / NON DENDRITIC CASE. 01381 01382 ZSPHE=PSNOWGRAN1(JJ,JST)/VGRAN1 01383 01384 IF(PSNOWHIST(JJ,JST) /= NVHIS1 .OR. PSNOWGRAN2(JJ,JST) < VDIAM2)THEN 01385 ZSPHE=ZSPHE+VSPHE2*ZVAP*PTSTEP 01386 ELSE 01387 01388 ! CAS HISTORIQUE=2 OU 3 ET GROS GRAINS 01389 ! SPHERICITE LIMITEE 01390 ! CASE HISTORY=2 OR 3 AND BIG GRAINS 01391 ! LIMITED SPHERICITY 01392 01393 ZSPHE=ZSPHE+ & 01394 VSPHE2*ZVAP*PTSTEP*EXP(MIN(0.,VDIAM3-PSNOWGRAN2(JJ,JST))/VDIAM6) 01395 ZSPHE=MIN(VSPHE3,ZSPHE) 01396 ENDIF 01397 PSNOWGRAN1(JJ,JST)=MIN(VGRAN1,ZSPHE*VGRAN1) 01398 ENDIF 01399 ELSEIF (ZGRADT < VGRAT2)THEN 01400 01401 ! 1.3.2 CALCUL METAMORPHOSE GRADIENT MOYEN/MODERATE (5-15). 01402 01403 01404 ZVAP=VDENT1*EXP(VVAP1/PSNOWTEMP(JJ,JST))*(ZGRADT)**VVAP2 01405 IF(PSNOWGRAN1(JJ,JST) < -UEPSI)THEN 01406 01407 ! 1.3.2.1 CAS DENDRITIQUE / DENDRITIC CASE. 01408 01409 ZDENT=-PSNOWGRAN1(JJ,JST)/VGRAN1 01410 ZSPHE=PSNOWGRAN2(JJ,JST)/VGRAN1 01411 ZDENT=ZDENT-ZVAP*PTSTEP 01412 ZSPHE=ZSPHE-ZVAP*PTSTEP 01413 IF(ZDENT < UEPSI)THEN 01414 ! EVOLUTION DE ZSGRAN1 ET ZSGRAN2 ET TEST PASSAGE 01415 ! DENDRITIQUE > NON DENDRITIQUE. 01416 PSNOWGRAN1(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01417 PSNOWGRAN2(JJ,JST)=VDIAM1-VDIAM5*MAX(ZSPHE,0.) 01418 ELSE 01419 PSNOWGRAN1(JJ,JST)=-ZDENT*VGRAN1 01420 PSNOWGRAN2(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01421 ENDIF 01422 ELSE 01423 01424 ! 1.3.2.2 CAS NON DENDRITIQUE / NON DENDRITIC. 01425 01426 ZSPHE=PSNOWGRAN1(JJ,JST)/VGRAN1 01427 ZSPHE=ZSPHE-ZVAP*PTSTEP 01428 PSNOWGRAN1(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01429 ENDIF 01430 ELSE 01431 01432 ! 1.3.3 CALCUL METAMORPHOSE FORT /HIGH GRADIENT 01433 01434 01435 ZVAP=VDENT1*EXP(VVAP1/PSNOWTEMP(JJ,JST))*(ZGRADT)**VVAP2 01436 IF(PSNOWGRAN1(JJ,JST) < -UEPSI)THEN 01437 01438 ! 1.3.3.1 CAS DENDRITIQUE / DENDRITIC CASE. 01439 01440 ZDENT=-PSNOWGRAN1(JJ,JST)/VGRAN1 01441 ZSPHE=PSNOWGRAN2(JJ,JST)/VGRAN1 01442 ZDENT=ZDENT-ZVAP*PTSTEP 01443 ! CAS NON DENDRITIQUE ET ANGULEUX. 01444 ZSPHE=ZSPHE-ZVAP*PTSTEP 01445 IF(ZDENT < UEPSI)THEN 01446 ! EVOLUTION DE SGRAN1 ET SGRAN2 ET TEST PASSAGE 01447 ! DENDRITIQUE > NON DENDRITIQUE. 01448 PSNOWGRAN1(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01449 PSNOWGRAN2(JJ,JST)=VDIAM1-VDIAM5*MAX(ZSPHE,0.) 01450 ELSE 01451 PSNOWGRAN1(JJ,JST)=-ZDENT*VGRAN1 01452 PSNOWGRAN2(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01453 ENDIF 01454 ELSEIF(PSNOWGRAN1(JJ,JST) > 0)THEN 01455 01456 ! 1.3.3.2 CAS NON DENDRITIQUE NON COMPLETEMENT ANGULEUX. 01457 ! NON DENDRITIC AND SPERICITY GT. 0 01458 01459 ZSPHE=PSNOWGRAN1(JJ,JST)/VGRAN1 01460 ZSPHE=ZSPHE-ZVAP*PTSTEP 01461 diam=PSNOWGRAN1(JJ,JST) 01462 PSNOWGRAN1(JJ,JST)=MAX(0.,ZSPHE*VGRAN1) 01463 ELSE 01464 01465 ! 1.3.3.3 CAS NON DENDRITIQUE ET ANGULEUX/DENDRITIC AND SPERICITY=0. 01466 ZDANGL = SNOW3L_MARBOUTY(PSNOWRHO(JJ,JST),PSNOWTEMP(JJ,JST),ZGRADT) 01467 diam=PSNOWGRAN2(JJ,JST) 01468 PSNOWGRAN2(JJ,JST) = PSNOWGRAN2(JJ,JST) + ZDANGL*VFI*PTSTEP 01469 ENDIF 01470 ENDIF 01471 ENDIF 01472 ENDDO 01473 ENDDO 01474 01475 !* 2. MISE A JOUR VARIABLES HISTORIQUES (CAS NON DENDRITIQUE). 01476 ! UPDATE OF THE HISTORICAL VARIABLES 01477 ! -------------------------------------------------------- 01478 DO JJ=1,SIZE(PSNOWRHO,1) 01479 DO JST=1,INLVLS_USE(JJ) 01480 IF(PSNOWGRAN1(JJ,JST) >= 0.)THEN 01481 IF(PSNOWGRAN1(JJ,JST) < VSPHE4 .AND. PSNOWHIST(JJ,JST) == 0)THEN 01482 PSNOWHIST(JJ,JST)=NVHIS1 01483 ELSEIF(VGRAN1-PSNOWGRAN1(JJ,JST) < VSPHE4 .AND. & 01484 PSNOWLIQ(JJ,JST)/(PSNOWDZ(JJ,JST)) > VTELV1)THEN 01485 IF(PSNOWHIST(JJ,JST) == 0) PSNOWHIST(JJ,JST)=NVHIS2 01486 IF(PSNOWHIST(JJ,JST) == NVHIS1) PSNOWHIST(JJ,JST)=NVHIS3 01487 ELSEIF(PSNOWTEMP(JJ,JST) < XTT)THEN 01488 IF(PSNOWHIST(JJ,JST) == NVHIS2) PSNOWHIST(JJ,JST)=NVHIS4 01489 IF(PSNOWHIST(JJ,JST) == NVHIS3) PSNOWHIST(JJ,JST)=NVHIS5 01490 ENDIF 01491 ENDIF 01492 ENDDO 01493 ENDDO 01494 IF (LHOOK) CALL DR_HOOK('SNOWCROMETAMO',1,ZHOOK_HANDLE) 01495 01496 END SUBROUTINE SNOWCROMETAMO 01497 01498 !#################################################################### 01499 !#################################################################### 01500 !#################################################################### 01501 ! 01502 !#################################################################### 01503 !#################################################################### 01504 !#################################################################### 01505 ! 01506 SUBROUTINE SNOWCROALB(TPTIME,OGLACIER,& 01507 PALBEDOSC,PSPECTRALALBEDO,PSNOWDZ,& 01508 PSNOWRHO,PPERMSNOWFRAC, & 01509 PSNOWGRAN1_TOP,PSNOWGRAN2_TOP,PSNOWAGE_TOP, & 01510 PSNOWGRAN1_BOT,PSNOWGRAN2_BOT,PSNOWAGE_BOT, & 01511 PPS, PZENITH, INLVLS_USE) 01512 ! 01513 !! PURPOSE 01514 !! ------- 01515 ! Calculate the snow surface albedo. Use the method of original 01516 ! Crocus which considers a specified spectral distribution of solar 01517 ! solar radiation (to be replaced by an input forcing when available) 01518 ! In addition to original crocus, the top 2 surface snow layers are 01519 ! considered in the calculation, using an arbitrary weighting, in order 01520 ! to avoid time discontinuities due to layers agregation 01521 ! Ageing depends on the presence of permanent snow cover 01522 ! 01523 USE MODD_CSTS, ONLY : XDAY 01524 USE MODD_SNOW_PAR, ONLY : XWCRN, XANSMAX, XANSMIN, XANS_TODRY, & 01525 XSNOWDMIN, XANS_T, XAGLAMIN, XAGLAMAX, & 01526 XALBICE1,XALBICE2,XALBICE3,& 01527 XRHOTHRESHOLD_ICE,XZ0ICEZ0SNOW,& 01528 XD1,XD2,XD3,XX,XVALB2,XVALB3,& 01529 XVALB4,XVALB5,XVALB6,XVALB7,XVALB8,& 01530 XVALB9,XVALB10,XVALB11,XVDIOP1, & 01531 XVRPRE1,XVRPRE2,XVAGING_NOGLACIER,& 01532 XVAGING_GLACIER,XVPRES1,& 01533 XVSPEC1,XVSPEC2,XVSPEC3,& 01534 XVW1,XVW2,XVD1,XVD2 01535 01536 USE MODD_TYPE_DATE_SURF, ONLY : DATE_TIME 01537 ! 01538 ! 01539 USE MODE_SNOW3L 01540 ! 01541 IMPLICIT NONE 01542 ! 01543 !* 0.1 declarations of arguments 01544 ! 01545 TYPE(DATE_TIME), INTENT(IN) :: TPTIME ! current date and time 01546 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 01547 ! initialise WGI=WSAT, 01548 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 01549 ! False = No specific treatment 01550 REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ,PPERMSNOWFRAC 01551 ! 01552 REAL,DIMENSION(:,:),INTENT(IN) :: PSNOWRHO ! For now only the 2 first layers are required 01553 ! 01554 REAL, DIMENSION(:), INTENT(INOUT) :: PALBEDOSC 01555 ! 01556 REAL, DIMENSION(:,:), INTENT(OUT) :: PSPECTRALALBEDO ! Albedo in the different spectral bands 01557 ! 01558 REAL, DIMENSION(:), INTENT(IN) :: PSNOWGRAN1_TOP,PSNOWGRAN2_TOP,PSNOWAGE_TOP, 01559 PSNOWGRAN1_BOT,PSNOWGRAN2_BOT,PSNOWAGE_BOT, PPS 01560 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 01561 ! 01562 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle for future use 01563 ! 01564 !* 0.2 declarations of local variables 01565 ! 01566 INTEGER :: JJ ! looping indexes 01567 ! 01568 REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZANSMIN,ZANSMAX 01569 ! 01570 REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZDIAM_TOP,ZDIAM_BOT 01571 REAL, DIMENSION(SIZE(PSNOWRHO),3) :: ZALB,ZALB_TOP,ZALB_BOT 01572 ! 01573 REAL, DIMENSION(SIZE(PALBEDOSC)) :: NVAGE1 01574 REAL :: ZAGE_NOW 01575 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01576 01577 !------------------------------------------------------------------------------- 01578 ! 01579 ! 0. Initialize: 01580 ! ------------------ 01581 ! 01582 01583 01584 ! PRINT*,XVAGING_NOGLACIER,XVAGING_GLACIER 01585 01586 01587 IF (LHOOK) CALL DR_HOOK('SNOWCROALB',0,ZHOOK_HANDLE) 01588 IF(OGLACIER)THEN 01589 ZANSMIN(:) = XAGLAMIN * PPERMSNOWFRAC(:) + XANSMIN * (1.0-PPERMSNOWFRAC(:)) 01590 ZANSMAX(:) = XAGLAMAX * PPERMSNOWFRAC(:) + XANSMAX * (1.0-PPERMSNOWFRAC(:)) 01591 NVAGE1(:) = XVAGING_GLACIER * PPERMSNOWFRAC(:) & 01592 + XVAGING_NOGLACIER * (1.0-PPERMSNOWFRAC(:)) 01593 ELSE 01594 ZANSMIN(:) = XANSMIN 01595 ZANSMAX(:) = XANSMAX 01596 NVAGE1(:) = XVAGING_NOGLACIER 01597 ENDIF 01598 01599 ! ! ! ! ! ! date computation for ageing effects 01600 ! ! ! ! ! CALL GREGODSTRATI(TPTIME%TDATE%YEAR,TPTIME%TDATE%MONTH,TPTIME%TDATE%DAY, & 01601 ! ! ! ! ! TPTIME%TIME,ZAGE_NOW) 01602 01603 ! coherence control 01604 ! to remove when initialization routines will be updated 01605 IF (MINVAL(PSNOWAGE_BOT)<0) THEN 01606 CALL ABOR1_SFX('FATAL ERROR in SNOWCRO: Snow layer age inconsistent : check initialization routine. !') 01607 END IF 01608 01609 ! ! ! ! ! ! should be moved with other time controls to not compute MAXVAL(PSNOWAGE_TOP) at each time step 01610 ! ! ! ! ! IF ((ZAGE_NOW - MAXVAL(PSNOWAGE_TOP))<-0.001) THEN 01611 ! ! ! ! ! WRITE(*,*),"ZAGE_NOW=",ZAGE_NOW 01612 ! ! ! ! ! WRITE(*,*),"MAXVAL(PSNOWAGE_TOP)=",MAXVAL(PSNOWAGE_TOP) 01613 ! ! ! ! ! CALL ABOR1_SFX(& 01614 ! ! ! ! ! 'FATAL ERROR in SNOWCRO: Snow layer date inconsistent with the current day !') 01615 ! ! ! ! ! END IF 01616 01617 DO JJ=1, SIZE(PALBEDOSC) 01618 IF (INLVLS_USE(JJ)==0) THEN 01619 ! case with no snow on the ground 01620 PALBEDOSC(JJ) = ZANSMIN(JJ) 01621 ELSE 01622 ! case with snow on the ground 01623 ! top layer 01624 01625 IF (PSNOWRHO(JJ,1) < XRHOTHRESHOLD_ICE) THEN 01626 ! Normal case (snow) 01627 01628 IF (PSNOWGRAN1_TOP(JJ)<0.) THEN 01629 ZDIAM_TOP(JJ) = -PSNOWGRAN1_TOP(JJ)*XD1/XX+(1.+PSNOWGRAN1_TOP(JJ)/XX)* & 01630 (PSNOWGRAN2_TOP(JJ)*XD2/XX+(1.-PSNOWGRAN2_TOP(JJ)/XX)*XD3) 01631 ZDIAM_TOP(JJ)=ZDIAM_TOP(JJ)/10000. 01632 ELSE 01633 ZDIAM_TOP(JJ)=MAX(0.0003,0.5*(1.+PSNOWGRAN1_TOP(JJ)/XX)*PSNOWGRAN2_TOP(JJ)) 01634 ENDIF 01635 ZALB_TOP(JJ,1)=MIN(XVALB2-XVALB3*SQRT(ZDIAM_TOP(JJ)),XVALB4) 01636 ZALB_TOP(JJ,2)=MAX(0.3,XVALB5-XVALB6*SQRT(ZDIAM_TOP(JJ))) 01637 ZDIAM_TOP(JJ)=MIN(ZDIAM_TOP(JJ),XVDIOP1) 01638 ZALB_TOP(JJ,3)=MAX(0.,XVALB7*ZDIAM_TOP(JJ)-XVALB8*SQRT(ZDIAM_TOP(JJ))+XVALB9) 01639 ! 01640 ! AGE CORRECTION ONLY FOR VISIBLE BAND 01641 ! 01642 ! ! ! ! ! ZALB_TOP(JJ,1)=MAX(XVALB11,ZALB_TOP(JJ,1)-MIN(MAX(PPS(JJ)/XVPRES1,XVRPRE1), & 01643 ! ! ! ! ! XVRPRE2)*XVALB10*(ZAGE_NOW-PSNOWAGE_TOP(JJ))/NVAGE1(JJ)) 01644 ZALB_TOP(JJ,1) = MAX(XVALB11,ZALB_TOP(JJ,1)-MIN(MAX(PPS(JJ)/XVPRES1,XVRPRE1), & 01645 XVRPRE2)*XVALB10*PSNOWAGE_TOP(JJ)/NVAGE1(JJ)) 01646 01647 01648 ELSE 01649 01650 ! Prescribed spectral albedo for surface ice 01651 ZALB_TOP(JJ,1)=XALBICE1 01652 ZALB_TOP(JJ,2)=XALBICE2 01653 ZALB_TOP(JJ,3)=XALBICE3 01654 01655 END IF 01656 01657 ! IF (INLVLS_USE(JJ)>=1) THEN 01658 IF (INLVLS_USE(JJ)>=2) THEN !modif ML 01659 ! second surface layer when it exists 01660 01661 IF (PSNOWRHO(JJ,2) < XRHOTHRESHOLD_ICE) THEN 01662 ! Normal case (snow) 01663 01664 IF(PSNOWGRAN1_BOT(JJ)<0.) THEN 01665 ZDIAM_BOT(JJ) = -PSNOWGRAN1_BOT(JJ)*XD1/XX+(1.+PSNOWGRAN1_BOT(JJ)/XX)* & 01666 (PSNOWGRAN2_BOT(JJ)*XD2/XX+(1.-PSNOWGRAN2_BOT(JJ)/XX)*XD3) 01667 ZDIAM_BOT(JJ)=ZDIAM_BOT(JJ)/10000. 01668 ELSE 01669 ZDIAM_BOT(JJ)=MAX(0.0003,0.5*(1.+PSNOWGRAN1_BOT(JJ)/XX)*PSNOWGRAN2_BOT(JJ)) 01670 ENDIF 01671 01672 ZALB_BOT(JJ,1)=MIN(XVALB2-XVALB3*SQRT(ZDIAM_BOT(JJ)),XVALB4) 01673 ZALB_BOT(JJ,2)=MAX(0.3,XVALB5-XVALB6*SQRT(ZDIAM_BOT(JJ))) 01674 ZDIAM_BOT(JJ)=MIN(ZDIAM_BOT(JJ),XVDIOP1) 01675 ZALB_BOT(JJ,3)=MAX(0.,XVALB7*ZDIAM_BOT(JJ)-XVALB8*SQRT(ZDIAM_BOT(JJ))+XVALB9) 01676 ! 01677 ! AGE CORRECTION ONLY FOR VISIBLE BAND 01678 01679 ! ! ! ! ! ZALB_BOT(JJ,1)=MAX(XVALB11,ZALB_BOT(JJ,1)-MIN(MAX(PPS(JJ)/XVPRES1,XVRPRE1), & 01680 ! ! ! ! ! XVRPRE2)*XVALB10*MIN(365.,ZAGE_NOW-PSNOWAGE_BOT(JJ))/NVAGE1(JJ)) 01681 01682 ZALB_BOT(JJ,1) = MAX(XVALB11,ZALB_BOT(JJ,1)-MIN(MAX(PPS(JJ)/XVPRES1,XVRPRE1), & 01683 XVRPRE2)*XVALB10*MIN(365.,PSNOWAGE_BOT(JJ))/NVAGE1(JJ)) 01684 ELSE 01685 01686 ! Prescribed spectral albedo for surface ice 01687 ZALB_BOT(JJ,1)=XALBICE1 01688 ZALB_BOT(JJ,2)=XALBICE2 01689 ZALB_BOT(JJ,3)=XALBICE3 01690 01691 ENDIF 01692 01693 ELSE 01694 01695 ! when it does not exist, the second surface layer gets top layer albedo 01696 ZALB_BOT(JJ,1)= ZALB_TOP(JJ,1) 01697 ZALB_BOT(JJ,2)= ZALB_TOP(JJ,2) 01698 ZALB_BOT(JJ,3)= ZALB_TOP(JJ,3) 01699 01700 ENDIF 01701 ! 01702 ! computation of spectral albedo over 3 bands taking into account the respective 01703 ! depths of top layers 01704 ! 01705 PSPECTRALALBEDO(JJ,1)=(XVW1*min(1.,PSNOWDZ(JJ)/XVD1)+ & 01706 XVW2*min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))) *ZALB_TOP(JJ,1) + & 01707 (XVW1*(1.-min(1.,PSNOWDZ(JJ)/XVD1))+ & 01708 XVW2*(1.- min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))))*ZALB_BOT(JJ,1) 01709 PSPECTRALALBEDO(JJ,2)=(XVW1*min(1.,PSNOWDZ(JJ)/XVD1)+ & 01710 XVW2*min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))) *ZALB_TOP(JJ,2) + & 01711 (XVW1*(1.-min(1.,PSNOWDZ(JJ)/XVD1))+ & 01712 XVW2*(1.- min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))))*ZALB_BOT(JJ,2) 01713 PSPECTRALALBEDO(JJ,3)=(XVW1*min(1.,PSNOWDZ(JJ)/XVD1)+ & 01714 XVW2*min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))) *ZALB_TOP(JJ,3) + & 01715 (XVW1*(1.-min(1.,PSNOWDZ(JJ)/XVD1))+ & 01716 XVW2*(1.- min(1.,max(0.,(PSNOWDZ(JJ)-XVD1)/XVD2))))*ZALB_BOT(JJ,3) 01717 ! 01718 ! arbitrarily specified spectral distribution 01719 ! to be changed when solar radiation distribution is an input variable 01720 PALBEDOSC(JJ)=XVSPEC1*PSPECTRALALBEDO(JJ,1)+XVSPEC2*PSPECTRALALBEDO(JJ,2)+XVSPEC3*PSPECTRALALBEDO(JJ,3) 01721 ! 01722 01723 ENDIF ! end case with snow on the ground 01724 ENDDO ! end loop grid points 01725 IF (LHOOK) CALL DR_HOOK('SNOWCROALB',1,ZHOOK_HANDLE) 01726 !------------------------------------------------------------------------------- 01727 ! 01728 END SUBROUTINE SNOWCROALB 01729 !#################################################################### 01730 !#################################################################### 01731 !#################################################################### 01732 SUBROUTINE SNOWCRORAD(TPTIME, OGLACIER, & 01733 PSW_RAD, PSNOWALB, PSNOWDZ, & 01734 PSNOWRHO, PALB, PRADSINK, PRADXS, & 01735 PSNOWGRAN1, PSNOWGRAN2, PSNOWAGE,PPS, & 01736 PZENITH, PPERMSNOWFRAC,INLVLS_USE ) 01737 ! 01738 !! PURPOSE 01739 !! ------- 01740 ! Calculate the transmission of shortwave (solar) radiation 01741 ! through the snowpack (using a form of Beer's Law: exponential 01742 ! decay of radiation with increasing snow depth). 01743 ! Needs a first calculation of the albedo to stay coherent with 01744 ! ISBA-ES ==> make sure to keep SNOWCRORAD coherent with SNOWCROALB 01745 ! 01746 USE MODD_CSTS, ONLY : XDAY 01747 USE MODD_SNOW_PAR, ONLY : XWCRN, XANSMAX, XANSMIN, XANS_TODRY, & 01748 XSNOWDMIN, XANS_T, XAGLAMIN, XAGLAMAX ,& 01749 XD1,XD2,XD3,XX,XVSPEC1,XVSPEC2,XVSPEC3 01750 USE MODD_TYPE_DATE_SURF, ONLY: DATE_TIME 01751 ! 01752 IMPLICIT NONE 01753 ! 01754 !* 0.1 declarations of arguments 01755 ! 01756 TYPE(DATE_TIME), INTENT(IN) :: TPTIME ! current date and time 01757 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 01758 ! initialise WGI=WSAT, 01759 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 01760 ! False = No specific treatment 01761 ! 01762 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD, PSNOWALB, PALB,PPERMSNOWFRAC 01763 ! 01764 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWDZ 01765 ! 01766 REAL, DIMENSION(:), INTENT(OUT) :: PRADXS 01767 ! 01768 REAL, DIMENSION(:,:), INTENT(OUT) :: PRADSINK 01769 ! 01770 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWGRAN1, PSNOWGRAN2, PSNOWAGE 01771 REAL, DIMENSION(:), INTENT(IN) :: PPS 01772 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 01773 REAL, DIMENSION(:), INTENT(IN) :: PZENITH 01774 !* 0.2 declarations of local variables 01775 ! 01776 INTEGER :: JJ,JST ! looping indexes 01777 ! 01778 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZRADTOT, ZANSMIN, ZANSMAX 01779 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZALB_NEW 01780 REAL :: ZPROJLAT 01781 REAL, DIMENSION(SIZE(PSNOWRHO,2)) :: ZDIAM 01782 REAL, DIMENSION(SIZE(PSNOWRHO,2),3) :: ZBETA 01783 REAL, DIMENSION(SIZE(PSNOWRHO,1),3) :: ZALB !albedo 3 bands 01784 REAL, DIMENSION(3) :: ZOPTICALPATH 01785 01786 ! calibration coefficients for exctinction computation 01787 REAL, PARAMETER :: PPVBETA1=1.92E-3,PPVBETA2=40.,PPVBETA3=1.098E-2, 01788 PPVBETA4=100., PPVBETA5=2000. 01789 01790 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01791 !------------------------------------------------------------------------------- 01792 ! 01793 ! 0. Initialization: 01794 ! ------------------ 01795 ! 01796 IF (LHOOK) CALL DR_HOOK('SNOWCRORAD',0,ZHOOK_HANDLE) 01797 PRADSINK(:,:)=0. 01798 ! 01799 IF(OGLACIER)THEN 01800 ZANSMIN(:) = XAGLAMIN * PPERMSNOWFRAC(:) + XANSMIN * (1.0-PPERMSNOWFRAC(:)) 01801 ZANSMAX(:) = XAGLAMAX * PPERMSNOWFRAC(:) + XANSMAX * (1.0-PPERMSNOWFRAC(:)) 01802 ELSE 01803 ZANSMIN(:) = XANSMIN 01804 ZANSMAX(:) = XANSMAX 01805 ENDIF 01806 ! 01807 ! 01808 ! 1. Computation of the new albedo (see SNOWCROALB): 01809 ! ----------------------------------- 01810 ! 01811 ! 01812 CALL SNOWCROALB(TPTIME,OGLACIER,& 01813 ZALB_NEW,ZALB,PSNOWDZ(:,1),PSNOWRHO(:,1:2), & 01814 PPERMSNOWFRAC,PSNOWGRAN1(:,1),PSNOWGRAN2(:,1), & 01815 PSNOWAGE(:,1),PSNOWGRAN1(:,2),PSNOWGRAN2(:,2),PSNOWAGE(:,2),& 01816 PPS, PZENITH, INLVLS_USE) 01817 01818 DO JJ=1, SIZE(PSW_RAD) 01819 ! Coefficient for taking into account the increase of path length of rays 01820 ! in snow due to zenithal angle 01821 ZPROJLAT=1./MAX(UEPSI,COS(PZENITH(JJ))) 01822 ! 01823 DO JST=1, INLVLS_USE(JJ) 01824 IF (PSNOWGRAN1(JJ,JST)<0.) THEN 01825 ZDIAM(JST)=-PSNOWGRAN1(JJ,JST)*XD1/XX+(1.+PSNOWGRAN1(JJ,JST)/XX)* & 01826 (PSNOWGRAN2(JJ,JST)*XD2/XX+(1.-PSNOWGRAN2(JJ,JST)/XX)*XD3) 01827 ZDIAM(JST)=ZDIAM(JST)/10000. 01828 ELSE 01829 ZDIAM(JST)=MAX(0.0003,0.5*(1.+PSNOWGRAN1(JJ,JST)/XX)*PSNOWGRAN2(JJ,JST)) 01830 ENDIF 01831 ENDDO ! end loop snow layers 01832 01833 ! 01834 ! 2. Extinction of net shortwave radiation 01835 ! ---------------------------------------- 01836 ! First calculates extinction coefficients fn of grain size and density 01837 ! then calculates exctinction in the layer and increases optical path length 01838 ! 01839 ! Initialize optical depth 01840 ZOPTICALPATH(1)=0. 01841 ZOPTICALPATH(2)=0. 01842 ZOPTICALPATH(3)=0. 01843 DO JST=1, INLVLS_USE(JJ) 01844 ZBETA(JST,1)=MAX(PPVBETA1*PSNOWRHO(JJ,JST)/SQRT(ZDIAM(JST)),PPVBETA2) 01845 ZBETA(JST,2)=MAX(PPVBETA3*PSNOWRHO(JJ,JST)/SQRT(ZDIAM(JST)),PPVBETA4) 01846 ZBETA(JST,3)=PPVBETA5 01847 ZOPTICALPATH(1)=ZOPTICALPATH(1)+ ZBETA(JST,1)*PSNOWDZ(JJ,JST) 01848 ZOPTICALPATH(2)=ZOPTICALPATH(2)+ ZBETA(JST,2)*PSNOWDZ(JJ,JST) 01849 ZOPTICALPATH(3)=ZOPTICALPATH(3)+ ZBETA(JST,3)*PSNOWDZ(JJ,JST) 01850 PRADSINK(JJ,JST) = -PSW_RAD(JJ)*(1.-PSNOWALB(JJ))/(1.-ZALB_NEW(JJ))* & 01851 ( XVSPEC1*(1.-ZALB(JJ,1))*EXP(-ZOPTICALPATH(1)*ZPROJLAT) & 01852 + XVSPEC2*(1.-ZALB(JJ,2))*EXP(-ZOPTICALPATH(2)*ZPROJLAT) & 01853 + XVSPEC3*(1.-ZALB(JJ,3))*EXP(-ZOPTICALPATH(3)*ZPROJLAT)) 01854 ENDDO ! end loop snow layers 01855 ! 01856 ! For thin snow packs, radiation might reach base of 01857 ! snowpack and the reflected energy can be absorbed by the bottom of snow layer: 01858 ! THIS PROCESS IS NOT SIMULATED 01859 ! 01860 ! 01861 ! 4. Excess radiation not absorbed by snowpack (W/m2)JJ 01862 ! ---------------------------------------------------- 01863 ! 01864 PRADXS(JJ) = -PRADSINK(JJ,INLVLS_USE(JJ)) 01865 01866 ENDDO !end loop grid points 01867 IF (LHOOK) CALL DR_HOOK('SNOWCRORAD',1,ZHOOK_HANDLE) 01868 ! 01869 !------------------------------------------------------------------------------- 01870 ! 01871 END SUBROUTINE SNOWCRORAD 01872 !#################################################################### 01873 !#################################################################### 01874 !#################################################################### 01875 SUBROUTINE SNOWCROTHRM(PSNOWRHO,PSCOND,PSNOWTEMP,PPS,PSNOWLIQ, & 01876 OCOND_GRAIN,OCOND_YEN) 01877 ! 01878 !! PURPOSE 01879 !! ------- 01880 ! Calculate snow thermal conductivity from 01881 ! Sun et al. 1999, J. of Geophys. Res., 104, 19587-19579 01882 ! (vapor) and Anderson, 1976, NOAA Tech. Rep. NWS 19 (snow). 01883 ! 01884 ! Upon activation of flag OCOND_YEN, use the Yen (1981) formula for thermal conductivity 01885 ! This formula was originally used in Crocus. 01886 ! 01887 ! 01888 USE MODD_CSTS,ONLY : XP00, XCONDI, XRHOLW 01889 ! 01890 IMPLICIT NONE 01891 ! 01892 !* 0.1 declarations of arguments 01893 ! 01894 REAL, DIMENSION(:), INTENT(IN) :: PPS 01895 ! 01896 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP, PSNOWRHO, PSNOWLIQ 01897 ! 01898 REAL, DIMENSION(:,:), INTENT(OUT) :: PSCOND 01899 ! 01900 LOGICAL, INTENT(IN) :: OCOND_GRAIN, OCOND_YEN 01901 ! 01902 !* 0.2 declarations of local variables 01903 ! 01904 INTEGER JJ, JST ! looping indexes 01905 ! 01906 ! ISBA-ES Thermal conductivity coefficients from Anderson (1976): 01907 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 01908 ! 01909 REAL, PARAMETER :: ZSNOWTHRMCOND1 = 0.02 ! [W/(m K)] 01910 REAL, PARAMETER :: ZSNOWTHRMCOND2 = 2.5E-6 ! [W m5/(kg2 K)] 01911 ! 01912 ! ISBA-ES Thermal conductivity: Implicit vapor diffn effects 01913 ! (sig only for new snow OR high altitudes) 01914 ! from Sun et al. (1999): based on data from Jordan (1991) 01915 ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 01916 ! 01917 REAL, PARAMETER :: ZSNOWTHRMCOND_AVAP = -0.06023 ! [W/(m K)] 01918 REAL, PARAMETER :: ZSNOWTHRMCOND_BVAP = -2.5425 ! (W/m) 01919 REAL, PARAMETER :: ZSNOWTHRMCOND_CVAP = -289.99 ! (K) 01920 !------------------------------------------------------------------------------- 01921 ! Crocus thermal conducitivity coefficient from Yen (1981) 01922 REAL, PARAMETER :: VRKZ6 = 1.88 01923 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01924 ! 01925 ! 1. Snow thermal conductivity 01926 ! ---------------------------- 01927 ! 01928 01929 IF (LHOOK) CALL DR_HOOK('SNOWCROTHRM',0,ZHOOK_HANDLE) 01930 DO JJ=1,SIZE(PSNOWRHO(:,:),1) 01931 DO JST=1,SIZE(PSNOWRHO(:,:),2) 01932 ! 01933 IF (OCOND_YEN) THEN 01934 PSCOND(JJ,JST) = XCONDI * EXP(VRKZ6 * LOG(PSNOWRHO(JJ,JST)/XRHOLW)) 01935 ELSE 01936 PSCOND(JJ,JST) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*PSNOWRHO(JJ,JST)* & 01937 PSNOWRHO(JJ,JST)) + & 01938 MAX(0.0,(ZSNOWTHRMCOND_AVAP+ & 01939 (ZSNOWTHRMCOND_BVAP/(PSNOWTEMP(JJ,JST)+ & 01940 ZSNOWTHRMCOND_CVAP)))*(XP00/PPS(JJ))) 01941 ENDIF 01942 ! 01943 ! Snow thermal conductivity is set to be above 0.04 W m-1 K-1 01944 IF(OCOND_GRAIN) THEN 01945 PSCOND(JJ,JST) = MAX(0.04,PSCOND(JJ,JST)) 01946 ! Snow thermal conductivity is annihilated in presence of liquid water 01947 IF(PSNOWLIQ(JJ,JST)>UEPSI) PSCOND(JJ,JST) =0.01*PSCOND(JJ,JST) 01948 ENDIF 01949 ! 01950 ENDDO ! end loop JST 01951 ENDDO ! end loop JST 01952 IF (LHOOK) CALL DR_HOOK('SNOWCROTHRM',1,ZHOOK_HANDLE) 01953 ! 01954 !------------------------------------------------------------------------------- 01955 ! 01956 END SUBROUTINE SNOWCROTHRM 01957 !#################################################################### 01958 !#################################################################### 01959 !#################################################################### 01960 SUBROUTINE SNOWCROEBUD(HSNOWRES, HIMPLICIT_WIND, & 01961 PPEW_A_COEF, PPEW_B_COEF, & 01962 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, & 01963 PSNOWDZMIN, & 01964 PZREF,PTS,PSNOWRHO,PSNOWLIQ,PSCAP,PSCOND1,PSCOND2, & 01965 PUREF,PEXNS,PEXNA,PDIRCOSZW,PVMOD, & 01966 PLW_RAD,PSW_RAD,PTA,PQA,PPS,PTSTEP, & 01967 PSNOWDZ1,PSNOWDZ2,PALBT,PZ0,PZ0EFF,PZ0H, & 01968 PSFCFRZ,PRADSINK,PHPSNOW, & 01969 PCT,PEMIST,PRHOA,PTSTERM1,PTSTERM2,PRA,PCDSNOW,PCHSNOW, & 01970 PQSAT,PDQSAT,PRSRA,PUSTAR2_IC, PRI, & 01971 PPET_A_COEF_T,PPEQ_A_COEF_T,PPET_B_COEF_T,PPEQ_B_COEF_T ) 01972 ! 01973 !! PURPOSE 01974 !! ------- 01975 ! Calculate surface energy budget linearization (terms) and turbulent 01976 ! exchange coefficients/resistance between surface and atmosphere. 01977 ! (Noilhan and Planton 1989; Giordani 1993; Noilhan and Mahfouf 1996) 01978 ! 01979 !! MODIFICATIONS 01980 !! ------------- 01981 !! Original A. Boone 01982 !! Modified by E. Brun (24/09/2012) : 01983 !! * Correction coupling coefficient for specific humidity in SNOWCROEBUD 01984 !! * PSFCFRZ(:) = 1.0 for systematic solid/vapor latent fluxes in SNOWCROEBUD 01985 !! Modified by B. Decharme 09/12 new wind implicitation 01986 ! 01987 USE MODD_SURF_PAR, ONLY : XUNDEF 01988 USE MODD_CSTS, ONLY : XCPD, XRHOLW, XSTEFAN, XLVTT, XLSTT, XRHOLW 01989 USE MODD_SNOW_PAR, ONLY : X_RI_MAX, XEMISSN 01990 ! 01991 USE MODE_THERMOS 01992 ! 01993 USE MODI_SURFACE_RI 01994 USE MODI_SURFACE_AERO_COND 01995 USE MODI_SURFACE_CD 01996 ! 01997 ! 01998 IMPLICIT NONE 01999 ! 02000 !* 0.1 declarations of arguments 02001 ! 02002 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN 02003 ! 02004 CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES ! type of sfc resistance 02005 ! DEFAULT=Louis (1979), standard ISBA 02006 ! method. Option to limit Ri number 02007 ! for very stable conditions 02008 ! 02009 CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option 02010 ! ! 'OLD' = direct 02011 ! ! 'NEW' = Taylor serie, order 1 02012 ! 02013 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, 02014 PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, 02015 PPEQ_B_COEF 02016 ! PPEW_A_COEF = wind coefficient (m2s/kg) 02017 ! PPEW_B_COEF = wind coefficient (m/s) 02018 ! PPET_A_COEF = A-air temperature coefficient 02019 ! PPET_B_COEF = B-air temperature coefficient 02020 ! PPEQ_A_COEF = A-air specific humidity coefficient 02021 ! PPEQ_B_COEF = B-air specific humidity coefficient 02022 ! 02023 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PTS, PSNOWDZ1, PSNOWDZ2, 02024 PRADSINK, PSNOWRHO, PSNOWLIQ, PSCAP, 02025 PSCOND1, PSCOND2, 02026 PZ0, PHPSNOW, 02027 PALBT, PZ0EFF, PZ0H 02028 ! 02029 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD, PLW_RAD, PTA, PQA, PPS, PRHOA 02030 ! 02031 REAL, DIMENSION(:), INTENT(IN) :: PUREF, PEXNS, PEXNA, PDIRCOSZW, PVMOD 02032 ! 02033 REAL, DIMENSION(:), INTENT(OUT) :: PTSTERM1, PTSTERM2, PEMIST, PRA, 02034 PCT, PSFCFRZ, PCDSNOW, PCHSNOW, 02035 PQSAT, PDQSAT, PRSRA 02036 ! 02037 REAL, DIMENSION(:), INTENT(OUT) :: PUSTAR2_IC, 02038 PPET_A_COEF_T, PPEQ_A_COEF_T, 02039 PPET_B_COEF_T, PPEQ_B_COEF_T 02040 ! 02041 REAL, DIMENSION(:), INTENT(OUT) :: PRI 02042 ! 02043 ! 02044 !* 0.2 declarations of local variables 02045 ! 02046 REAL, DIMENSION(SIZE(PTS)) :: ZAC, ZRI, 02047 ZSCONDA, ZA, ZB, ZC, 02048 ZCDN, ZSNOWDZM1, ZSNOWDZM2, 02049 ZVMOD, ZUSTAR2, ZTS3, ZLVT, 02050 Z_CCOEF 02051 REAL, DIMENSION(SIZE(PTS)) :: ZSNOWEVAPX ! 02052 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02053 ! 02054 INTEGER :: JJ ! looping indexes 02055 !------------------------------------------------------------------------------- 02056 ! 02057 ! 1. New saturated specific humidity and derrivative: 02058 ! --------------------------------------------------- 02059 ! 02060 IF (LHOOK) CALL DR_HOOK('SNOWCROEBUD',0,ZHOOK_HANDLE) 02061 ! 02062 ZRI (:) = XUNDEF 02063 ! 02064 PQSAT (:) = QSATI(PTS(:),PPS(:)) 02065 PDQSAT(:) = DQSATI(PTS(:),PPS(:),PQSAT(:)) 02066 ! 02067 ! 02068 ! 2. Surface properties: 02069 ! ---------------------- 02070 ! It might be of interest to use snow-specific roughness 02071 ! or a temperature dependence on emissivity: 02072 ! but for now, use ISBA defaults. 02073 ! 02074 PEMIST(:) = XEMISSN 02075 ! 02076 ! 2. Computation of resistance and drag coefficient 02077 ! ------------------------------------------------- 02078 ! 02079 CALL SURFACE_RI(PTS, PQSAT, PEXNS, PEXNA, PTA, PQA, & 02080 PZREF, PUREF, PDIRCOSZW, PVMOD, ZRI ) 02081 ! 02082 ! Simple adaptation of method by Martin and Lejeune (1998) 02083 ! to apply a lower limit to turbulent transfer coef 02084 ! by defining a maximum Richarson number for stable 02085 ! conditions: 02086 ! 02087 IF(HSNOWRES=='RIL') ZRI(:) = MIN(X_RI_MAX, ZRI(:)) 02088 ! 02089 PRI(:)=ZRI(:) 02090 ! 02091 ! Surface aerodynamic resistance for heat transfers 02092 ! 02093 CALL SURFACE_AERO_COND(ZRI, PZREF, PUREF, PVMOD, PZ0, PZ0H, ZAC, PRA, PCHSNOW) 02094 ! 02095 PRSRA(:) = PRHOA(:) / PRA(:) 02096 ! 02097 ! For atmospheric model coupling: 02098 ! 02099 CALL SURFACE_CD(ZRI, PZREF, PUREF, PZ0EFF, PZ0H, PCDSNOW, ZCDN) 02100 ! 02101 ! 02102 ! Modify flux-form implicit coupling coefficients: 02103 ! - wind components: 02104 ! 02105 IF(HIMPLICIT_WIND=='OLD')THEN 02106 ! old implicitation 02107 ZUSTAR2(:) = (PCDSNOW(:)*PVMOD(:)*PPEW_B_COEF(:))/ & 02108 (1.0-PRHOA(:)*PCDSNOW(:)*PVMOD(:)*PPEW_A_COEF(:)) 02109 ELSE 02110 ! new implicitation 02111 ZUSTAR2(:) = (PCDSNOW(:)*PVMOD(:)*(2.*PPEW_B_COEF(:)-PVMOD(:))) & 02112 / (1.0-2.0*PRHOA(:)*PCDSNOW(:)*PVMOD(:)*PPEW_A_COEF(:)) 02113 ENDIF 02114 ! 02115 ZVMOD(:) = PRHOA(:)*PPEW_A_COEF(:)*ZUSTAR2(:) + PPEW_B_COEF(:) 02116 ZVMOD(:) = MAX(ZVMOD(:),0.) 02117 ! 02118 WHERE(PPEW_A_COEF(:)/= 0.) 02119 ZUSTAR2(:) = MAX( ( ZVMOD(:) - PPEW_B_COEF(:) ) / (PRHOA(:)*PPEW_A_COEF(:)), 0.) 02120 ENDWHERE 02121 ! 02122 ! implicit wind friction 02123 ZUSTAR2(:) = MAX(ZUSTAR2(:),0.) 02124 ! 02125 PUSTAR2_IC(:) = ZUSTAR2(:) 02126 ! 02127 ! 3. Calculate linearized surface energy budget components: 02128 ! --------------------------------------------------------- 02129 ! To prevent numerical difficulties for very thin snow 02130 ! layers, limit the grid "thinness": this is important as 02131 ! layers become vanishing thin: 02132 ! 02133 ZSNOWDZM1(:) = MAX(PSNOWDZ1(:), PSNOWDZMIN) 02134 ZSNOWDZM2(:) = MAX(PSNOWDZ2(:), PSNOWDZMIN) 02135 ! 02136 ! Surface thermal inertia: 02137 ! 02138 PCT(:) = 1.0/(PSCAP(:)*ZSNOWDZM1(:)) 02139 ! 02140 ! Fraction of surface frozen (sublimation) with the remaining 02141 ! fraction being liquid (evaporation): 02142 ! 02143 PSFCFRZ(:) = 1.0 02144 ! 02145 ! Thermal conductivity between uppermost and lower snow layers: 02146 ! 02147 ZSCONDA(:) = (ZSNOWDZM1(:)*PSCOND1(:) + ZSNOWDZM2(:)*PSCOND2(:))/ & 02148 (ZSNOWDZM1(:)+ZSNOWDZM2(:)) 02149 ! 02150 ! Transform implicit coupling coefficients: 02151 ! Note, surface humidity is 100% over snow. 02152 ! 02153 ! - specific humidity: 02154 ! 02155 Z_CCOEF(:) = 1.0 - PPEQ_A_COEF(:)*PRSRA(:) 02156 ! 02157 PPEQ_A_COEF_T(:) = - PPEQ_A_COEF(:)*PRSRA(:)*PDQSAT(:)/Z_CCOEF(:) 02158 ! 02159 PPEQ_B_COEF_T(:) = ( PPEQ_B_COEF(:) - PPEQ_A_COEF(:)*PRSRA(:)*(PQSAT(:) - & 02160 PDQSAT(:)*PTS(:)) )/Z_CCOEF(:) 02161 ! 02162 ! - air temperature: 02163 ! (assumes A and B correspond to potential T): 02164 ! 02165 Z_CCOEF(:) = (1.0 - PPET_A_COEF(:)*PRSRA(:))/PEXNA(:) 02166 ! 02167 PPET_A_COEF_T(:) = - PPET_A_COEF(:)*PRSRA(:)/(PEXNS(:)*Z_CCOEF(:)) 02168 ! 02169 PPET_B_COEF_T(:) = PPET_B_COEF(:)/Z_CCOEF(:) 02170 ! 02171 ! 02172 ! Energy budget solution terms: 02173 ! 02174 ZTS3(:) = PEMIST(:) * XSTEFAN * PTS(:)**3 02175 ZLVT(:) = (1.-PSFCFRZ(:))*XLVTT + PSFCFRZ(:)*XLSTT 02176 ! 02177 ZA(:) = 1. / PTSTEP + PCT(:) * (4. * ZTS3(:) + & 02178 PRSRA(:) * ZLVT(:) * (PDQSAT(:) - PPEQ_A_COEF_T(:)) & 02179 + PRSRA(:) * XCPD * ( (1./PEXNS(:))-(PPET_A_COEF_T(:)/PEXNA(:)) ) & 02180 + (2*ZSCONDA(:)/(ZSNOWDZM2(:)+ZSNOWDZM1(:))) ) 02181 ! 02182 ZB(:) = 1. / PTSTEP + PCT(:) * (3. * ZTS3(:) + & 02183 PRSRA(:) * PDQSAT(:) * ZLVT(:) ) 02184 ! 02185 ZC(:) = PCT(:) * (PRSRA(:) * XCPD * PPET_B_COEF_T(:)/PEXNA(:) + PSW_RAD(:) * & 02186 (1. - PALBT(:)) + PEMIST(:)*PLW_RAD(:) - PRSRA(:) * & 02187 ZLVT(:) * (PQSAT(:)-PPEQ_B_COEF_T(:)) & 02188 + PHPSNOW(:) + PRADSINK(:) ) 02189 ! 02190 ! 02191 ! Coefficients needed for implicit solution 02192 ! of linearized surface energy budget: 02193 ! 02194 PTSTERM2(:) = 2*ZSCONDA(:)*PCT(:)/(ZA(:)*(ZSNOWDZM2(:)+ZSNOWDZM1(:))) 02195 ! 02196 PTSTERM1(:) = (PTS(:)*ZB(:) + ZC(:))/ZA(:) 02197 IF (LHOOK) CALL DR_HOOK('SNOWCROEBUD',1,ZHOOK_HANDLE) 02198 ! 02199 !------------------------------------------------------------------------------- 02200 ! 02201 END SUBROUTINE SNOWCROEBUD 02202 !#################################################################### 02203 !#################################################################### 02204 !#################################################################### 02205 SUBROUTINE SNOWCROSOLVT(PTSTEP,PSNOWDZMIN, & 02206 PSNOWDZ,PSCOND,PSCAP,PTG, & 02207 PSOILCOND,PD_G, & 02208 PRADSINK,PCT,PTERM1,PTERM2, & 02209 PPET_A_COEF_T,PPEQ_A_COEF_T, & 02210 PPET_B_COEF_T,PPEQ_B_COEF_T, & 02211 PTA_IC, PQA_IC, & 02212 PGBAS,PSNOWTEMP,PSNOWFLUX, & 02213 INLVLS_USE ) 02214 ! 02215 !! PURPOSE 02216 !! ------- 02217 ! This subroutine solves the 1-d diffusion of 'ZSNOWTEMP' using a 02218 ! layer averaged set of equations which are time differenced 02219 ! using the backward-difference scheme (implicit). 02220 ! The eqs are solved rapidly by taking advantage of the 02221 ! fact that the matrix is tridiagonal. This is a very 02222 ! general routine and can be used for the 1-d diffusion of any 02223 ! quantity as long as the diffusity is not a function of the 02224 ! quantity being diffused. Aaron Boone 8-98. Soln to the eq: 02225 ! 02226 ! c dQ d k dQ dS 02227 ! -- = -- -- - -- 02228 ! dt dx dx dx 02229 ! 02230 ! where k = k(x) (thermal conductivity), c = c(x) (heat capacity) 02231 ! as an eg. for temperature/heat eq. S is a sink (-source) term. 02232 ! Diffusivity is k/c 02233 ! 02234 !! MODIFICATIONS 02235 !! ------------- 02236 !! Original A. Boone 02237 !! 05/2011: Brun Special treatment to tackle the variable number 02238 !! of snow layers 02239 ! 02240 USE MODD_CSTS,ONLY : XTT 02241 ! 02242 IMPLICIT NONE 02243 ! 02244 !* 0.1 declarations of arguments 02245 ! 02246 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN 02247 ! 02248 REAL, DIMENSION(:), INTENT(IN) :: PTG, PSOILCOND, PD_G, 02249 PCT, PTERM1, PTERM2 02250 02251 ! 02252 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ, PSCOND, PSCAP, 02253 PRADSINK 02254 ! 02255 REAL, DIMENSION(:), INTENT(IN) :: PPET_A_COEF_T, PPEQ_A_COEF_T, 02256 PPET_B_COEF_T, PPEQ_B_COEF_T 02257 ! 02258 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP 02259 ! 02260 REAL, DIMENSION(:), INTENT(OUT) :: PGBAS, PSNOWFLUX, PTA_IC, PQA_IC 02261 ! 02262 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 02263 ! 02264 !* 0.2 declarations of local variables 02265 ! 02266 ! 02267 INTEGER :: JJ ! looping indexes 02268 ! 02269 INTEGER :: INLVLS 02270 ! 02271 REAL, DIMENSION(SIZE(PTG)) :: ZGBAS, ZSNOWTEMP_DELTA 02272 ! 02273 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSNOWTEMP, ZDTERM, ZCTERM, 02274 ZFRCV, ZAMTRX, ZBMTRX, 02275 ZCMTRX 02276 ! 02277 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZWORK1, ZWORK2, ZDZDIF, 02278 ZSNOWDZM 02279 ! 02280 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)-1) :: ZSNOWTEMP_M, 02281 ZFRCV_M, ZAMTRX_M, 02282 ZBMTRX_M, ZCMTRX_M 02283 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02284 !------------------------------------------------------------------------------- 02285 ! 02286 ! 0. Initialize: 02287 ! ------------------ 02288 ! 02289 IF (LHOOK) CALL DR_HOOK('SNOWCROSOLVT',0,ZHOOK_HANDLE) 02290 ZSNOWTEMP(:,:) = PSNOWTEMP(:,:) 02291 INLVLS = SIZE(PSNOWDZ(:,:),2) 02292 ! 02293 ! 02294 ! 1. Calculate tri-diagnoal matrix coefficients: 02295 ! ---------------------------------------------- 02296 ! For heat transfer, assume a minimum grid 02297 ! thickness (to prevent numerical 02298 ! problems for very thin snow cover): 02299 ! 02300 DO JJ=1,SIZE(PTG) 02301 ZSNOWDZM(JJ,1:INLVLS_USE(JJ)) = MAX(PSNOWDZ(JJ,1:INLVLS_USE(JJ)), PSNOWDZMIN) 02302 ZDZDIF(JJ,1:INLVLS_USE(JJ)-1) = ZSNOWDZM(JJ,1:INLVLS_USE(JJ)-1) + & 02303 ZSNOWDZM(JJ,2:INLVLS_USE(JJ)) 02304 ZDZDIF(JJ,INLVLS_USE(JJ)) = ZSNOWDZM(JJ,INLVLS_USE(JJ)) + PD_G(JJ) 02305 ZWORK1(JJ,1:INLVLS_USE(JJ)) = ZSNOWDZM(JJ,1:INLVLS_USE(JJ))* & 02306 PSCOND(JJ,1:INLVLS_USE(JJ)) 02307 ZWORK2(JJ,1:INLVLS_USE(JJ)-1) = ZSNOWDZM(JJ,2:INLVLS_USE(JJ))*PSCOND(JJ,2:INLVLS_USE(JJ)) 02308 ZWORK2(JJ,INLVLS_USE(JJ)) = PD_G(JJ)*PSOILCOND(JJ) 02309 ! 02310 ZDTERM(JJ,1:INLVLS_USE(JJ)) = 2.0*(ZWORK1(JJ,1:INLVLS_USE(JJ))+ & 02311 ZWORK2(JJ,1:INLVLS_USE(JJ)))/ & 02312 (ZDZDIF(JJ,1:INLVLS_USE(JJ))* & 02313 ZDZDIF(JJ,1:INLVLS_USE(JJ))) 02314 ! 02315 ZCTERM(JJ,1:INLVLS_USE(JJ)) = PSCAP(JJ,1:INLVLS_USE(JJ))* & 02316 ZSNOWDZM(JJ,1:INLVLS_USE(JJ))/PTSTEP 02317 ENDDO 02318 ! 02319 ! 2. Set up tri-diagonal matrix 02320 ! ----------------------------- 02321 ! 02322 ! Upper BC 02323 ! 02324 ZAMTRX(:,1) = 0.0 02325 ZBMTRX(:,1) = 1./(PCT(:)*PTSTEP) 02326 ZCMTRX(:,1) = -PTERM2(:)*ZBMTRX(:,1) 02327 ZFRCV(:,1) = PTERM1(:)*ZBMTRX(:,1) 02328 ! 02329 ! 02330 ! 02331 DO JJ=1,SIZE(PTG) 02332 ! Interior Grid 02333 ZAMTRX(JJ,2:INLVLS_USE(JJ)-1) = -ZDTERM(JJ,1:INLVLS_USE(JJ)-2) 02334 ZBMTRX(JJ,2:INLVLS_USE(JJ)-1) = ZCTERM(JJ,2:INLVLS_USE(JJ)-1) + & 02335 ZDTERM(JJ,1:INLVLS_USE(JJ)-2) + & 02336 ZDTERM(JJ,2:INLVLS_USE(JJ)-1) 02337 ZCMTRX(JJ,2:INLVLS_USE(JJ)-1) = -ZDTERM(JJ,2:INLVLS_USE(JJ)-1) 02338 ZFRCV(JJ,2:INLVLS_USE(JJ)-1) = ZCTERM(JJ,2:INLVLS_USE(JJ)-1)* & 02339 PSNOWTEMP(JJ,2:INLVLS_USE(JJ)-1) - & 02340 (PRADSINK(JJ,1:INLVLS_USE(JJ)-2)- & 02341 PRADSINK(JJ,2:INLVLS_USE(JJ)-1)) 02342 !Lower BC 02343 ZAMTRX(JJ,INLVLS_USE(JJ)) = -ZDTERM(JJ,INLVLS_USE(JJ)-1) 02344 ZBMTRX(JJ,INLVLS_USE(JJ)) = ZCTERM(JJ,INLVLS_USE(JJ)) + & 02345 ZDTERM(JJ,INLVLS_USE(JJ)-1) +ZDTERM(JJ,INLVLS_USE(JJ)) 02346 ZCMTRX(JJ,INLVLS_USE(JJ)) = 0.0 02347 ZFRCV(JJ,INLVLS_USE(JJ)) = ZCTERM(JJ,INLVLS_USE(JJ))* & 02348 PSNOWTEMP(JJ,INLVLS_USE(JJ)) + ZDTERM(JJ,INLVLS_USE(JJ))*PTG(JJ) & 02349 - (PRADSINK(JJ,INLVLS_USE(JJ)-1)-PRADSINK(JJ,INLVLS_USE(JJ))) 02350 ENDDO 02351 ! 02352 ! - - ------------------------------------------------- 02353 ! 02354 ! 4. Compute solution vector 02355 ! -------------------------- 02356 ! 02357 CALL TRIDIAG_GROUND_SNOWCRO(ZAMTRX,ZBMTRX,ZCMTRX,ZFRCV,ZSNOWTEMP, & 02358 INLVLS_USE,0) 02359 ! 02360 ! Heat flux between surface and 2nd snow layers: (W/m2) 02361 ! 02362 PSNOWFLUX(:) = ZDTERM(:,1)*(ZSNOWTEMP(:,1) - ZSNOWTEMP(:,2)) 02363 ! 02364 ! 02365 ! 5. Snow melt case 02366 ! ----------------- 02367 ! If melting in uppermost layer, assume surface layer 02368 ! temperature at freezing point and re-evaluate lower 02369 ! snowpack temperatures. This is done as it is most likely 02370 ! most signigant melting will occur within a time step in surface layer. 02371 ! Surface energy budget (and fluxes) will 02372 ! be re-calculated (outside of this routine): 02373 ! 02374 ZAMTRX_M(:,1) = 0.0 02375 ZBMTRX_M(:,1) = ZCTERM(:,2) + ZDTERM(:,1) + ZDTERM(:,2) 02376 ZCMTRX_M(:,1) = -ZDTERM(:,2) 02377 ZFRCV_M(:,1) = ZCTERM(:,2)*PSNOWTEMP(:,2) + XTT*ZDTERM(:,1) - & 02378 (PRADSINK(:,1)-PRADSINK(:,2)) 02379 ! 02380 DO JJ=1,SIZE(PTG) 02381 ZAMTRX_M(JJ,2:INLVLS_USE(JJ)-1) = ZAMTRX(JJ,3:INLVLS_USE(JJ)) 02382 ZBMTRX_M(JJ,2:INLVLS_USE(JJ)-1) = ZBMTRX(JJ,3:INLVLS_USE(JJ)) 02383 ZCMTRX_M(JJ,2:INLVLS_USE(JJ)-1) = ZCMTRX(JJ,3:INLVLS_USE(JJ)) 02384 ZFRCV_M(JJ,2:INLVLS_USE(JJ)-1) = ZFRCV(JJ,3:INLVLS_USE(JJ)) 02385 ZSNOWTEMP_M(JJ,2:INLVLS_USE(JJ)-1) = PSNOWTEMP(JJ,3:INLVLS_USE(JJ)) 02386 ENDDO 02387 ! 02388 CALL TRIDIAG_GROUND_SNOWCRO(ZAMTRX_M,ZBMTRX_M,ZCMTRX_M,ZFRCV_M,ZSNOWTEMP_M, & 02389 INLVLS_USE,1) 02390 ! 02391 ! If melting for 2 consecuative time steps, then replace current T-profile 02392 ! with one assuming T=Tf in surface layer: 02393 ! 02394 ZSNOWTEMP_DELTA(:) = 0.0 02395 ! 02396 WHERE(ZSNOWTEMP(:,1) > XTT .AND. PSNOWTEMP(:,1) >= XTT) 02397 PSNOWFLUX(:) = ZDTERM(:,1)*(XTT - ZSNOWTEMP_M(:,1)) 02398 ZSNOWTEMP_DELTA(:) = 1.0 02399 END WHERE 02400 ! 02401 DO JJ=1,SIZE(PTG) 02402 ZSNOWTEMP(JJ,2:INLVLS_USE(JJ)) = ZSNOWTEMP_DELTA(JJ)* & 02403 ZSNOWTEMP_M(JJ,1:INLVLS_USE(JJ)-1) + & 02404 (1.0-ZSNOWTEMP_DELTA(JJ))*ZSNOWTEMP(JJ,2:INLVLS_USE(JJ)) 02405 ENDDO 02406 ! 02407 ! 02408 ! 6. Lower boundary flux: 02409 ! ----------------------- 02410 ! NOTE: evaluate this term assuming the snow layer 02411 ! can't exceed the freezing point as this adjustment 02412 ! is made in melting routine. Then must adjust temperature 02413 ! to conserve energy: 02414 ! 02415 DO JJ=1, SIZE(PTG) 02416 ZGBAS(JJ) = ZDTERM(JJ,INLVLS_USE(JJ))* & 02417 (ZSNOWTEMP(JJ,INLVLS_USE(JJ)) -PTG(JJ)) 02418 PGBAS(JJ) = ZDTERM(JJ,INLVLS_USE(JJ))* & 02419 (MIN(XTT,ZSNOWTEMP(JJ,INLVLS_USE(JJ)))-PTG(JJ)) 02420 ZSNOWTEMP(JJ,INLVLS_USE(JJ)) = ZSNOWTEMP(JJ,INLVLS_USE(JJ)) + & 02421 (ZGBAS(JJ)-PGBAS(JJ))/ZCTERM(JJ,INLVLS_USE(JJ)) 02422 ENDDO 02423 ! 02424 02425 ! 02426 ! 7. Update temperatute profile in time: 02427 ! -------------------------------------- 02428 ! 02429 DO JJ=1, SIZE(PTG) 02430 PSNOWTEMP(JJ,1:INLVLS_USE(JJ)) = ZSNOWTEMP(JJ,1:INLVLS_USE(JJ)) 02431 ENDDO 02432 ! 02433 ! 02434 ! 8. Compute new (implicit) air T and specific humidity 02435 ! ----------------------------------------------------- 02436 ! 02437 PTA_IC(:) = PPET_B_COEF_T(:) + PPET_A_COEF_T(:)* PSNOWTEMP(:,1) 02438 ! 02439 PQA_IC(:) = PPEQ_B_COEF_T(:) + PPEQ_A_COEF_T(:)* PSNOWTEMP(:,1) 02440 IF (LHOOK) CALL DR_HOOK('SNOWCROSOLVT',1,ZHOOK_HANDLE) 02441 ! 02442 ! 02443 END SUBROUTINE SNOWCROSOLVT 02444 !#################################################################### 02445 !#################################################################### 02446 !#################################################################### 02447 SUBROUTINE SNOWCROMELT(PSCAP,PSNOWTEMP,PSNOWDZ, & 02448 PSNOWRHO,PSNOWLIQ,INLVLS_USE ) 02449 ! 02450 ! 02451 !! PURPOSE 02452 !! ------- 02453 ! Calculate snow melt (resulting from surface fluxes, ground fluxes, 02454 ! or internal shortwave radiation absorbtion). It is used to 02455 ! augment liquid water content, maintain temperatures 02456 ! at or below freezing, and possibly reduce the mass 02457 ! or compact the layer(s). 02458 ! 02459 ! 02460 USE MODD_CSTS,ONLY : XTT, XLMTT, XRHOLW, XRHOLI 02461 ! 02462 USE MODE_SNOW3L 02463 ! 02464 IMPLICIT NONE 02465 ! 02466 !* 0.1 declarations of arguments 02467 ! 02468 REAL, DIMENSION(:,:), INTENT(IN) :: PSCAP 02469 ! 02470 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWRHO, 02471 PSNOWLIQ 02472 ! 02473 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 02474 ! 02475 !* 0.2 declarations of local variables 02476 ! 02477 INTEGER :: JJ, JST ! looping indexes 02478 ! 02479 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZPHASE, ZCMPRSFACT, 02480 ZSNOWLWE, 02481 ZSNOWMELT, ZSNOWTEMP 02482 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02483 !------------------------------------------------------------------------------- 02484 ! 02485 ! 0. Initialize: 02486 ! --------------------------- 02487 ! 02488 IF (LHOOK) CALL DR_HOOK('SNOWCROMELT',0,ZHOOK_HANDLE) 02489 DO JJ=1,SIZE(PSNOWDZ,1) 02490 DO JST=1,INLVLS_USE(JJ) 02491 ZPHASE(JJ,JST) = 0.0 02492 ZCMPRSFACT(JJ,JST) = 0.0 02493 ZSNOWLWE(JJ,JST) = 0.0 02494 ZSNOWMELT(JJ,JST) = 0.0 02495 ZSNOWTEMP(JJ,JST) = 0.0 02496 ENDDO 02497 ENDDO 02498 ! 02499 ! 1. Determine amount of melt in each layer: 02500 ! ------------------------------------------ 02501 ! 02502 DO JJ=1,SIZE(PSNOWDZ,1) 02503 DO JST=1, INLVLS_USE(JJ) 02504 ! 02505 ! Total Liquid equivalent water content of snow (m): 02506 ! 02507 ZSNOWLWE(JJ,JST) = PSNOWRHO(JJ,JST)*PSNOWDZ(JJ,JST)/XRHOLW 02508 ! 02509 ! Melt snow if excess energy and snow available: 02510 ! Phase change (J/m2) 02511 ! 02512 ZPHASE(JJ,JST) = MIN(PSCAP(JJ,JST)*MAX(0.0, PSNOWTEMP(JJ,JST) - XTT)* & 02513 PSNOWDZ(JJ,JST), & 02514 MAX(0.0,ZSNOWLWE(JJ,JST)-PSNOWLIQ(JJ,JST))*XLMTT*XRHOLW) 02515 ! 02516 ! Update snow liq water content and temperature if melting: 02517 ! liquid water available for next layer from melting of snow 02518 ! which is assumed to be leaving the current layer (m): 02519 ! 02520 ZSNOWMELT(JJ,JST) = ZPHASE(JJ,JST)/(XLMTT*XRHOLW) 02521 ! 02522 ! Cool off snow layer temperature due to melt: 02523 ! 02524 ZSNOWTEMP(JJ,JST) = PSNOWTEMP(JJ,JST) - ZPHASE(JJ,JST)/(PSCAP(JJ,JST)*PSNOWDZ(JJ,JST)) 02525 ! 02526 ! Difference with ISBA_ES: ZMELTXS should never be different of 0. 02527 ! because of the introduction of the tests in LLAYERGONE 02528 ! 02529 PSNOWTEMP(JJ,JST) = ZSNOWTEMP(JJ,JST) 02530 ! The control below should be suppressed after further tests 02531 IF (PSNOWTEMP(JJ,JST)-XTT > UEPSI) THEN 02532 WRITE(*,*) 'pb dans MELT PSNOWTEMP(JJ,JST) >XTT:', JJ,JST, PSNOWTEMP(JJ,JST) 02533 CALL ABOR1_SFX('SNOWCRO: pb dans MELT') 02534 ENDIF 02535 ! 02536 ! Loss of snowpack depth: (m) and liquid equiv (m): 02537 ! Compression factor for melt loss: this decreases 02538 ! layer thickness and increases density thereby leaving 02539 ! total SWE constant. 02540 ! 02541 ! Difference with ISBA_ES: All melt is considered to decrease the depth 02542 ! without consideration to the irreducible content 02543 ! 02544 ZCMPRSFACT(JJ,JST) = (ZSNOWLWE(JJ,JST)-(PSNOWLIQ(JJ,JST)+ZSNOWMELT(JJ,JST)))& 02545 / (ZSNOWLWE(JJ,JST)-PSNOWLIQ(JJ,JST)) 02546 PSNOWDZ(JJ,JST) = PSNOWDZ(JJ,JST)*ZCMPRSFACT(JJ,JST) 02547 PSNOWRHO(JJ,JST) = ZSNOWLWE(JJ,JST)*XRHOLW/PSNOWDZ(JJ,JST) 02548 ! 02549 ! 2. Add snow melt to current snow liquid water content: 02550 ! ------------------------------------------------------ 02551 ! 02552 PSNOWLIQ(JJ,JST) = PSNOWLIQ(JJ,JST) + ZSNOWMELT(JJ,JST) 02553 ! 02554 ENDDO ! loop JST active snow layers 02555 ENDDO ! loop JJ grid points 02556 IF (LHOOK) CALL DR_HOOK('SNOWCROMELT',1,ZHOOK_HANDLE) 02557 ! 02558 ! 02559 END SUBROUTINE SNOWCROMELT 02560 !#################################################################### 02561 !#################################################################### 02562 !#################################################################### 02563 SUBROUTINE SNOWCROREFRZ(PTSTEP,PRR, & 02564 PSNOWRHO,PSNOWTEMP,PSNOWDZ,PSNOWLIQ, & 02565 PTHRUFAL, PSCAP, PLEL3L,INLVLS_USE) 02566 ! 02567 ! 02568 !! PURPOSE 02569 !! ------- 02570 ! Calculate any freezing/refreezing of liquid water in the snowpack. 02571 ! Also, calculate liquid water transmission and snow runoff. 02572 ! Refreezing causes densification of a layer. 02573 ! 02574 ! 02575 USE MODD_CSTS, ONLY : XTT, XLMTT, XRHOLW, XCI,XRHOLI 02576 USE MODD_SNOW_PAR, ONLY : XSNOWDMIN 02577 ! 02578 USE MODE_SNOW3L 02579 ! 02580 IMPLICIT NONE 02581 ! 02582 !* 0.1 declarations of arguments 02583 ! 02584 REAL, INTENT(IN) :: PTSTEP 02585 ! 02586 REAL, DIMENSION(:), INTENT(IN) :: PRR 02587 ! 02588 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWLIQ, PSNOWRHO 02589 ! 02590 REAL, DIMENSION(:), INTENT(INOUT) :: PTHRUFAL 02591 ! 02592 ! 02593 ! modifs_EB layers 02594 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 02595 REAL, DIMENSION(:,:), INTENT(IN) :: PSCAP 02596 REAL, DIMENSION(:), INTENT(IN) :: PLEL3L 02597 02598 ! 02599 !* 0.2 declarations of local variables 02600 ! 02601 INTEGER :: JJ, JST ! looping indexes 02602 ! 02603 INTEGER :: INLVLS ! maximum snow layers number 02604 ! 02605 ! 02606 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZPHASE, 02607 ZSNOWLIQ, ZSNOWRHO, 02608 ZWHOLDMAX, ZSNOWDZ, 02609 ZSNOWTEMP 02610 ! 02611 REAL, DIMENSION(SIZE(PSNOWRHO,1),0:SIZE(PSNOWRHO,2)) :: ZFLOWLIQ 02612 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02613 ! 02614 !------------------------------------------------------------------------------- 02615 ! 02616 ! 0. Initialize: 02617 ! -------------- 02618 ! 02619 IF (LHOOK) CALL DR_HOOK('SNOWCROREFRZ',0,ZHOOK_HANDLE) 02620 INLVLS = SIZE(PSNOWDZ,2) 02621 DO JJ=1,SIZE(PSNOWDZ,1) 02622 DO JST=1,INLVLS_USE(JJ) 02623 ZSNOWRHO(JJ,JST) = PSNOWRHO(JJ,JST) 02624 ZSNOWTEMP(JJ,JST) = PSNOWTEMP(JJ,JST) 02625 ZWHOLDMAX(JJ,JST) = SNOWCROHOLD(PSNOWRHO(JJ,JST),PSNOWLIQ(JJ,JST),PSNOWDZ(JJ,JST)) 02626 ENDDO 02627 ENDDO 02628 ! 02629 DO JJ=1,SIZE(PSNOWDZ,1) ! loop JJ grid points 02630 ! 02631 ! 1. Increases Liquid Water of top layer from rain 02632 ! --------------------------------------------- 02633 ! 02634 ! Rainfall (m) initialises the liquid flow whih feeds the top layer 02635 ! and evaporation/condensation are taken into account 02636 ! 02637 IF(INLVLS_USE(JJ) > 0.) THEN 02638 ZFLOWLIQ(JJ,0)= PRR(JJ)*PTSTEP/XRHOLW 02639 ZFLOWLIQ(JJ,0)= MAX(0.,ZFLOWLIQ(JJ,0)- PLEL3L(JJ)*PTSTEP/(XLVTT*XRHOLW)) 02640 ELSE 02641 ZFLOWLIQ(JJ,0)=0 02642 ENDIF 02643 ! 02644 ! 02645 DO JST=1,INLVLS_USE(JJ) ! loop JST active snow layers 02646 02647 02648 ! 2. Increases Liquid Water from the upper layers flow (or rain for top layer) 02649 ! ----------------------------- 02650 ! 02651 PSNOWLIQ(JJ,JST) = PSNOWLIQ(JJ,JST) + ZFLOWLIQ(JJ,JST-1) 02652 ! 02653 ! 3. Freezes liquid water in any cold layers 02654 ! --------------------------------------- 02655 ! 02656 ! Calculate the maximum possible refreezing 02657 ! 02658 ZPHASE(JJ,JST) = MIN(PSCAP(JJ,JST)* & 02659 MAX(0.0, XTT - ZSNOWTEMP(JJ,JST))*PSNOWDZ(JJ,JST), & 02660 PSNOWLIQ(JJ,JST)*XLMTT*XRHOLW) 02661 ! 02662 ! Reduce liquid content if freezing occurs: 02663 ! 02664 ZSNOWLIQ(JJ,JST) = PSNOWLIQ(JJ,JST) - ZPHASE(JJ,JST)/(XLMTT*XRHOLW) 02665 ! 02666 ! Warm layer and reduce liquid if freezing occurs: 02667 ! 02668 ZSNOWDZ(JJ,JST) = MAX(XSNOWDMIN/INLVLS, PSNOWDZ(JJ,JST)) 02669 ! 02670 ! 02671 ! Difference with ISBA-ES: a possible cooling of current refreezing water 02672 ! is taken into account to calculate temperature change 02673 ! 02674 PSNOWTEMP(JJ,JST)=XTT+(ZSNOWTEMP(JJ,JST)-XTT)* & 02675 (ZSNOWRHO(JJ,JST)*ZSNOWDZ(JJ,JST)-(PSNOWLIQ(JJ,JST)-ZFLOWLIQ(JJ,JST-1))*XRHOLW)/ & 02676 (ZSNOWRHO(JJ,JST)*ZSNOWDZ(JJ,JST)-(ZSNOWLIQ(JJ,JST)-ZFLOWLIQ(JJ,JST-1))*XRHOLW) & 02677 + ZPHASE(JJ,JST)/(XCI*(ZSNOWRHO(JJ,JST)*ZSNOWDZ(JJ,JST)- & 02678 (ZSNOWLIQ(JJ,JST)-ZFLOWLIQ(JJ,JST-1))*XRHOLW)) 02679 02680 ! 02681 ! 4. Calculate flow from the excess of holding capacity 02682 ! -------------------------------------------------------------- 02683 ! Any water in excess of the maximum holding space for liquid water 02684 ! amount is drained into next layer down. 02685 ! 02686 ZFLOWLIQ(JJ,JST) = MAX(0.,ZSNOWLIQ(JJ,JST)-ZWHOLDMAX(JJ,JST)) 02687 ! 02688 ZSNOWLIQ(JJ,JST) = ZSNOWLIQ(JJ,JST) - ZFLOWLIQ(JJ,JST) 02689 ! 02690 ! 5. Density is adjusted to conserve the mass 02691 ! -------------------------------------------------------------- 02692 ZSNOWRHO(JJ,JST) = (ZSNOWRHO(JJ,JST)*PSNOWDZ(JJ,JST) & 02693 -(ZFLOWLIQ(JJ,JST)-ZFLOWLIQ(JJ,JST-1))*XRHOLW)/ZSNOWDZ(JJ,JST) 02694 ! keeps snow density below ice density 02695 IF(ZSNOWRHO(JJ,JST)> XRHOLI) THEN 02696 PSNOWDZ(JJ,JST)=PSNOWDZ(JJ,JST)*ZSNOWRHO(JJ,JST)/XRHOLI 02697 ZSNOWRHO(JJ,JST)=XRHOLI 02698 ENDIF 02699 02700 ! 6. Update thickness and density and any freezing: 02701 ! ---------------------------------------------- 02702 ! 02703 PSNOWRHO(JJ,JST) = ZSNOWRHO(JJ,JST) 02704 PSNOWLIQ(JJ,JST) = ZSNOWLIQ(JJ,JST) 02705 ! 02706 ENDDO ! loop JST active snow layers 02707 ! 02708 ! Any remaining throughflow after freezing is available to 02709 ! the soil for infiltration or surface runoff (m). 02710 ! I.E. This is the amount of water leaving the snowpack: 02711 ! Rate water leaves the snowpack [kg/(m2 s)]: 02712 ! 02713 PTHRUFAL(JJ) = PTHRUFAL(JJ) + ZFLOWLIQ(JJ,INLVLS_USE(JJ))*XRHOLW/PTSTEP 02714 ENDDO ! loop JJ grid points 02715 IF (LHOOK) CALL DR_HOOK('SNOWCROREFRZ',1,ZHOOK_HANDLE) 02716 ! 02717 ! 02718 !------------------------------------------------------------------------------- 02719 ! 02720 END SUBROUTINE SNOWCROREFRZ 02721 !#################################################################### 02722 !#################################################################### 02723 !#################################################################### 02724 SUBROUTINE SNOWCROFLUX(PSNOWTEMP,PSNOWDZ,PEXNS,PEXNA, & 02725 PUSTAR2_IC, & 02726 PTSTEP,PALBT,PSW_RAD,PEMIST,PLWUPSNOW, & 02727 PLW_RAD,PTA,PSFCFRZ,PQA,PHPSNOW, & 02728 PSNOWTEMPO1,PSNOWFLUX,PCT,PRADSINK, & 02729 PQSAT,PDQSAT,PRSRA, & 02730 PRN,PH,PGFLUX,PLES3L,PLEL3L,PEVAP, & 02731 PUSTAR ) 02732 ! 02733 ! 02734 !! PURPOSE 02735 !! ------- 02736 ! Calculate the surface fluxes (atmospheric/surface). 02737 ! (Noilhan and Planton 1989; Noilhan and Mahfouf 1996) 02738 ! 02739 ! 02740 USE MODD_CSTS,ONLY : XSTEFAN, XCPD, XLSTT, XLVTT, XTT 02741 ! 02742 USE MODE_THERMOS 02743 ! 02744 IMPLICIT NONE 02745 ! 02746 !* 0.1 declarations of arguments 02747 ! 02748 REAL, INTENT(IN) :: PTSTEP 02749 ! 02750 REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ, PSNOWTEMPO1, PSNOWFLUX, PCT, 02751 PRADSINK, PEXNS, PEXNA 02752 ! 02753 REAL, DIMENSION(:), INTENT(IN) :: PALBT, PSW_RAD, PEMIST, PLW_RAD, 02754 PTA, PSFCFRZ, PQA, 02755 PHPSNOW, PQSAT, PDQSAT, PRSRA, 02756 PUSTAR2_IC 02757 ! 02758 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWTEMP 02759 ! 02760 REAL, DIMENSION(:), INTENT(OUT) :: PRN, PH, PGFLUX, PLES3L, PLEL3L, 02761 PEVAP, PLWUPSNOW, PUSTAR 02762 ! 02763 ! 02764 !* 0.2 declarations of local variables 02765 ! 02766 REAL, DIMENSION(SIZE(PSNOWDZ)) :: ZEVAPC, ZLE, ZSNOWTEMP, ZSMSNOW, ZGFLUX, 02767 ZDELTAT, ZSNOWTO3 02768 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02769 ! 02770 !------------------------------------------------------------------------------- 02771 ! 02772 ! 0. Initialize: 02773 ! -------------- 02774 ! 02775 IF (LHOOK) CALL DR_HOOK('SNOWCROFLUX',0,ZHOOK_HANDLE) 02776 ZSNOWTEMP(:) = PSNOWTEMP(:) 02777 ZLE(:) = 0.0 02778 ZSMSNOW(:) = 0.0 02779 ZGFLUX(:) = 0.0 02780 ! 02781 ZSNOWTO3(:) = PSNOWTEMPO1(:) ** 3 ! to save some CPU time, store this 02782 ! 02783 ! 1. Flux calculations when melt not occuring at surface (W/m2): 02784 ! -------------------------------------------------------------- 02785 ! 02786 ! 02787 ZDELTAT(:) = PSNOWTEMP(:) - PSNOWTEMPO1(:) ! surface T time change: 02788 ! 02789 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * ZSNOWTO3(:)*( PSNOWTEMPO1(:) + 4.* ZDELTAT(:) ) 02790 ! 02791 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02792 ! 02793 PH(:) = PRSRA(:) * XCPD * (PSNOWTEMP(:)/PEXNS(:) - PTA(:)/PEXNA(:)) 02794 ! 02795 ZEVAPC(:) = PRSRA(:) * ( (PQSAT(:) - PQA(:)) + PDQSAT(:)*ZDELTAT(:) ) 02796 ! 02797 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02798 ! 02799 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02800 ! 02801 ZLE(:) = PLES3L(:) + PLEL3L(:) 02802 ! 02803 PGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02804 ! 02805 ! 02806 ! 2. Initial melt adjustment 02807 ! -------------------------- 02808 ! If energy avalabile to melt snow, then recalculate fluxes 02809 ! at the freezing point and add residual heat to layer 02810 ! average heat. 02811 ! 02812 ! A) If temperature change is > 0 and passes freezing point this timestep, 02813 ! then recalculate fluxes at freezing point and excess energy 02814 ! will be used outside of this routine to change snow heat content: 02815 ! 02816 !write (*,*) 'attention test LFLUX traitement XTT supprime!' 02817 WHERE (PSNOWTEMP > XTT .AND. PSNOWTEMPO1 < XTT) 02818 ! 02819 ZDELTAT(:) = XTT - PSNOWTEMPO1(:) 02820 ! 02821 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * ZSNOWTO3(:)*( PSNOWTEMPO1(:) + 4.* ZDELTAT(:) ) 02822 ! 02823 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02824 ! 02825 PH(:) = PRSRA(:) * XCPD * (XTT/PEXNS(:) - PTA(:)/PEXNA(:)) 02826 ! 02827 ZEVAPC(:) = PRSRA(:) * ( (PQSAT(:) - PQA(:)) + PDQSAT(:)*ZDELTAT(:) ) 02828 ! 02829 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02830 ! 02831 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02832 ! 02833 ZLE(:) = PLES3L(:) + PLEL3L(:) 02834 ! 02835 ZGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02836 ! 02837 ZSMSNOW(:) = PGFLUX(:) - ZGFLUX(:) 02838 ! 02839 PGFLUX(:) = ZGFLUX(:) 02840 ! 02841 ! This will be used to change heat content of snow: 02842 ! 02843 ZSNOWTEMP(:) = PSNOWTEMP(:) - ZSMSNOW(:)*PTSTEP*PCT(:) 02844 ! 02845 END WHERE 02846 ! 02847 ! 3. Ongoing melt adjustment: explicit solution 02848 ! --------------------------------------------- 02849 ! If temperature change is 0 and at freezing point this timestep, 02850 ! then recalculate fluxes and surface temperature *explicitly* 02851 ! as this is *exact* for snow at freezing point (Brun, Martin) 02852 ! 02853 WHERE(PSNOWTEMP(:) > XTT .AND. PSNOWTEMPO1(:) >= XTT ) 02854 ! 02855 PLWUPSNOW(:) = PEMIST(:) * XSTEFAN * (XTT ** 4) 02856 ! 02857 PRN(:) = (1. - PALBT(:)) * PSW_RAD(:) + PEMIST(:) * PLW_RAD(:) - PLWUPSNOW(:) 02858 ! 02859 PH(:) = PRSRA(:) * XCPD * (XTT/PEXNS(:) - PTA(:)/PEXNA(:)) 02860 ! 02861 ZEVAPC(:) = PRSRA(:) * (PQSAT(:) - PQA(:)) 02862 ! 02863 PLES3L(:) = PSFCFRZ(:) * XLSTT * ZEVAPC(:) 02864 ! 02865 PLEL3L(:) = (1.-PSFCFRZ(:))* XLVTT * ZEVAPC(:) 02866 ! 02867 ZLE(:) = PLES3L(:) + PLEL3L(:) 02868 ! 02869 PGFLUX(:) = PRN(:) - PH(:) - ZLE(:) + PHPSNOW(:) 02870 ! 02871 ZSNOWTEMP(:) = XTT + PTSTEP*PCT(:)*(PGFLUX(:) + PRADSINK(:) - PSNOWFLUX(:)) 02872 ! 02873 END WHERE 02874 ! 02875 ! 4. Update surface temperature: 02876 ! ------------------------------ 02877 ! 02878 PSNOWTEMP(:) = ZSNOWTEMP(:) 02879 ! 02880 ! 5. Final evaporative flux (kg/m2/s) 02881 ! 02882 PEVAP(:) = ZEVAPC(:) 02883 !write(*,*) 'Flux surface:',PGFLUX(1),PRN(1),PH(1), ZLE(1), PHPSNOW(1) 02884 ! 02885 ! 5. Friction velocity 02886 ! -------------------- 02887 ! 02888 PUSTAR(:) = SQRT(PUSTAR2_IC(:)) 02889 ! 02890 IF (LHOOK) CALL DR_HOOK('SNOWCROFLUX',1,ZHOOK_HANDLE) 02891 ! 02892 !------------------------------------------------------------------------------- 02893 ! 02894 END SUBROUTINE SNOWCROFLUX 02895 !#################################################################### 02896 !#################################################################### 02897 !#################################################################### 02898 SUBROUTINE SNOWCROEVAPN(PLES3L,PTSTEP,PSNOWTEMP, & 02899 PSNOWRHO,PSNOWDZ,PEVAPCOR, & 02900 PSNOWHMASS ) 02901 ! 02902 ! 02903 !! PURPOSE 02904 !! ------- 02905 ! Remove mass from uppermost snow layer in response to 02906 ! evaporation (liquid) and sublimation. 02907 ! 02908 !! MODIFICATIONS 02909 !! ------------- 02910 !! Original A. Boone 02911 !! 05/2011: E. Brun Takes only into account sublimation and solid 02912 !! condensation. Evaporation and liquid condensation 02913 !! are taken into account in SNOWCROREFRZ 02914 ! 02915 USE MODD_CSTS, ONLY : XLVTT, XRHOLW, XLSTT, XLMTT, XCI, XTT 02916 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES 02917 ! 02918 IMPLICIT NONE 02919 ! 02920 !* 0.1 declarations of arguments 02921 ! 02922 REAL, INTENT(IN) :: PTSTEP 02923 ! 02924 REAL, DIMENSION(:), INTENT(IN) :: PSNOWTEMP 02925 ! 02926 REAL, DIMENSION(:), INTENT(IN) :: PLES3L ! (W/m2) 02927 ! 02928 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ, PSNOWHMASS, 02929 PEVAPCOR 02930 ! 02931 !* 0.2 declarations of local variables 02932 ! 02933 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWEVAPS, ZSNOWEVAP, ZSNOWEVAPX, 02934 ZSNOWDZ, ZEVAPCOR 02935 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02936 ! ZEVAPCOR = for vanishingy thin snow cover, 02937 ! allow any excess evaporation 02938 ! to be extracted from the soil 02939 ! to maintain an accurate water 02940 ! balance [kg/(m2 s)] 02941 ! 02942 ! 02943 !------------------------------------------------------------------------------- 02944 ! 02945 ! 0. Initialize: 02946 ! -------------- 02947 ! 02948 IF (LHOOK) CALL DR_HOOK('SNOWCROEVAPN',0,ZHOOK_HANDLE) 02949 ZEVAPCOR(:) = 0.0 02950 ZSNOWEVAPS(:) = 0.0 02951 ZSNOWEVAP(:) = 0.0 02952 ZSNOWEVAPX(:) = 0.0 02953 ZSNOWDZ(:) = 0.0 02954 ! 02955 ! 02956 ! 02957 WHERE(PSNOWDZ > 0.0) 02958 ! 02959 !! 1. Sublimation/condensation of snow ice 02960 ! ---------------------------------------- 02961 ! Reduce layer thickness and total snow depth 02962 ! if sublimation: add to correction term if potential 02963 ! sublimation exceeds available snow cover. 02964 ! 02965 ZSNOWEVAPS(:) = PLES3L(:)*PTSTEP/(XLSTT*PSNOWRHO(:)) 02966 ZSNOWDZ(:) = PSNOWDZ(:) - ZSNOWEVAPS(:) 02967 PSNOWDZ(:) = MAX(0.0, ZSNOWDZ(:)) 02968 ZEVAPCOR(:) = ZEVAPCOR(:) + MAX(0.0,-ZSNOWDZ(:))*PSNOWRHO(:)/PTSTEP 02969 ! 02970 ! Total heat content change due to snowfall and sublimation (added here): 02971 ! (for budget calculations): 02972 ! 02973 PSNOWHMASS(:) = PSNOWHMASS(:) - PLES3L(:)*(PTSTEP/XLSTT)* & 02974 ( XCI*(PSNOWTEMP(:)-XTT)- XLMTT) 02975 END WHERE 02976 ! 02977 ! 3. Update evaporation correction term: 02978 ! -------------------------------------- 02979 ! 02980 PEVAPCOR(:) = PEVAPCOR(:) + ZEVAPCOR(:) 02981 IF (LHOOK) CALL DR_HOOK('SNOWCROEVAPN',1,ZHOOK_HANDLE) 02982 ! 02983 !------------------------------------------------------------------------------- 02984 ! 02985 END SUBROUTINE SNOWCROEVAPN 02986 !#################################################################### 02987 !#################################################################### 02988 !#################################################################### 02989 SUBROUTINE SNOWCROGONE(PTSTEP,PLEL3L,PLES3L,PSNOWRHO, & 02990 PSNOWHEAT,PRADSINK_2D,PEVAPCOR,PTHRUFAL,PGRNDFLUX, & 02991 PGFLUXSNOW,PSNOWDZ,PSNOWLIQ,PSNOWTEMP,PRADXS, & 02992 PRR,INLVLS_USE) 02993 ! 02994 ! 02995 !! PURPOSE 02996 !! ------- 02997 ! Account for the case when the last trace of snow melts 02998 ! during a time step: ensure mass and heat balance of 02999 ! snow AND underlying surface. 03000 ! Original A. Boone 03001 ! 05/2011: E. Brun Takes into account sublimation and PGRNDFLUX 03002 ! Adds rain and evaporation/liquid condensation 03003 ! in PTHRUFAL 03004 ! 03005 USE MODD_CSTS,ONLY : XTT, XLSTT, XLVTT 03006 ! 03007 IMPLICIT NONE 03008 ! 03009 !* 0.1 declarations of arguments 03010 ! 03011 REAL, INTENT(IN) :: PTSTEP 03012 ! 03013 REAL, DIMENSION(:), INTENT(IN) :: PLEL3L, PLES3L, PGFLUXSNOW,PRR 03014 ! 03015 REAL, DIMENSION(:,:), INTENT(IN) :: PRADSINK_2D 03016 ! 03017 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWHEAT 03018 ! 03019 REAL, DIMENSION(:), INTENT(INOUT) :: PGRNDFLUX, PRADXS 03020 ! 03021 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWLIQ, PSNOWTEMP 03022 ! 03023 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL ! melt water [kg/(m2 s)] 03024 ! 03025 REAL, DIMENSION(:), INTENT(OUT) :: PEVAPCOR ! [kg/(m2 s)] 03026 ! PEVAPCOR = for vanishingy thin snow cover, 03027 ! allow any excess evaporation 03028 ! to be extracted from the soil 03029 ! to maintain an accurate water 03030 ! balance. 03031 ! 03032 INTEGER, DIMENSION(:), INTENT(INOUT) :: INLVLS_USE 03033 ! 03034 !* 0.2 declarations of local variables 03035 ! 03036 REAL, DIMENSION(SIZE(PLES3L)) :: PRADSINK 03037 ! 03038 INTEGER :: JJ 03039 ! 03040 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWHEATC 03041 INTEGER, DIMENSION(SIZE(PLES3L)) :: ISNOWGONE_DELTA 03042 REAL(KIND=JPRB) :: ZHOOK_HANDLE 03043 ! 03044 !------------------------------------------------------------------------------- 03045 ! 03046 ! 0. Initialize: 03047 ! -------------- 03048 ! 03049 IF (LHOOK) CALL DR_HOOK('SNOWCROGONE',0,ZHOOK_HANDLE) 03050 PEVAPCOR(:) = 0.0 03051 PTHRUFAL(:) = 0.0 03052 ! 03053 DO JJ=1, SIZE(PRADSINK) 03054 PRADSINK(JJ) = PRADSINK_2D(JJ,INLVLS_USE(JJ)) 03055 ZSNOWHEATC(JJ) = SUM(PSNOWHEAT(JJ,1:INLVLS_USE(JJ))) !total heat content (J m-2) 03056 END DO 03057 ! 03058 ISNOWGONE_DELTA(:) = 1 03059 ! 03060 ! 1. Simple test to see if snow vanishes: 03061 ! --------------------------------------- 03062 ! If so, set thicknesses (and therefore mass and heat) and liquid content 03063 ! to zero, and adjust fluxes of water, evaporation and heat into underlying 03064 ! surface. 03065 ! 03066 ! takes into account the heat content corresponding to the occasional 03067 ! sublimation and then PGRNDFLUX 03068 ! 03069 ZSNOWHEATC(:) = ZSNOWHEATC(:)+MAX(0.,PLES3L(:)*PTSTEP/XLSTT)* & 03070 XLMTT 03071 ! 03072 WHERE(PGFLUXSNOW(:) + PRADSINK(:) - PGRNDFLUX(:) >= (-ZSNOWHEATC(:)/PTSTEP) ) 03073 PGRNDFLUX(:) = PGFLUXSNOW(:) + (ZSNOWHEATC(:)/PTSTEP) 03074 PEVAPCOR(:) = PLES3L(:)/XLSTT 03075 PRADXS(:) = 0.0 03076 ISNOWGONE_DELTA(:) = 0 ! FLAG...if=0 then snow vanishes, else=1 03077 END WHERE 03078 ! 03079 ! 2. Final update of snow state and computation of corresponding flow 03080 ! Only if snow vanishes 03081 ! ----------------------------- 03082 PTHRUFAL(:) = 0. 03083 DO JJ=1, SIZE(PRADSINK) 03084 IF(ISNOWGONE_DELTA(JJ) == 0 ) THEN 03085 PTHRUFAL(JJ) = PTHRUFAL(JJ) + SUM(PSNOWRHO(JJ,1:INLVLS_USE(JJ))* & 03086 PSNOWDZ(JJ,1:INLVLS_USE(JJ)))/PTSTEP 03087 ! takes into account rain and condensation/evaporation 03088 PTHRUFAL(JJ)= PTHRUFAL(JJ) + PRR(JJ)-PLEL3L(JJ)/XLVTT 03089 PSNOWTEMP(JJ,:) = XTT 03090 PSNOWDZ(JJ,:) = 0. 03091 PSNOWLIQ(JJ,:) = 0. 03092 INLVLS_USE(JJ)=0 03093 ENDIF 03094 ENDDO 03095 IF (LHOOK) CALL DR_HOOK('SNOWCROGONE',1,ZHOOK_HANDLE) 03096 ! 03097 ! 03098 END SUBROUTINE SNOWCROGONE 03099 !#################################################################### 03100 !#################################################################### 03101 !#################################################################### 03102 SUBROUTINE SNOWCROEVAPGONE(PSNOWHEAT,PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ,& 03103 PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE, INLVLS_USE) 03104 ! 03105 !! PURPOSE 03106 !! ------- 03107 ! 03108 ! If all snow in uppermost layer evaporates/sublimates, re-distribute 03109 ! grid (below assumes very thin snowpacks so layer-thicknesses are 03110 ! constant). 03111 ! Original A. Boone 03112 ! 05/2011: E. Brun Takes into account previous changes in the energy 03113 ! content 03114 ! 03115 ! 03116 USE MODD_CSTS, ONLY : XTT, XRHOLW, XLMTT 03117 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN, XRHOSMAX_ES 03118 USE MODE_SNOW3L 03119 USE MODD_SNOW_METAMO 03120 USE MODD_TYPE_DATE_SURF 03121 ! 03122 IMPLICIT NONE 03123 ! 03124 !* 0.1 declarations of arguments 03125 ! 03126 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO ! snow density profile (kg/m3) 03127 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ ! snow layer thickness profile (m) 03128 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT ! snow heat content/enthalpy (J/m2) 03129 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWGRAN1 ! snow grain parameter 1 (-) 03130 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWGRAN2 ! snow grain parameter 2 (-) 03131 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHIST ! snow grain historical variable (-) 03132 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWAGE ! Snow grain age 03133 ! 03134 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP ! snow temperature profile (K) 03135 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWLIQ ! snow liquid water profile (m) 03136 ! 03137 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 03138 ! 03139 !* 0.2 declarations of local variables 03140 ! 03141 INTEGER :: JJ,JST ! loop control 03142 INTEGER :: JNLVLS 03143 ! 03144 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOWHEAT_1D ! total heat content (J/m2) 03145 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOWRHO_1D ! total snowpack average density (kg/m3) 03146 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOW ! total snow depth (m) 03147 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSCAP ! Snow layer heat capacity (J/K/m3) 03148 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZNDENT ! Number of dendritic layers (-) 03149 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZNVIEU ! Number of non dendritic layers (-) 03150 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOWAGE_1D ! total snowpack average 03151 !age (days) 03152 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) ::ZSNOWGRAN1N, 03153 ZSNOWGRAN2N,ZSNOWHISTN 03154 REAL(KIND=JPRB) :: ZHOOK_HANDLE 03155 03156 ! 03157 !------------------------------------------------------------------------------- 03158 ! 03159 ! Initialize: 03160 ! 03161 IF (LHOOK) CALL DR_HOOK('SNOWCROEVAPGONE',0,ZHOOK_HANDLE) 03162 JNLVLS = SIZE(PSNOWDZ,2) 03163 ZSNOWHEAT_1D(:) = 0.0; ZSNOWRHO_1D(:) = 0.0; ZSNOW(:) = 0.0; ZSCAP(:) = 0.0 03164 ! 03165 ! First, determine where uppermost snow layer has completely 03166 ! evaporated/sublimated (as it becomes thin): 03167 ! 03168 ZSNOWHEAT_1D(:) = 0. 03169 ZSNOW(:) = 0. 03170 ZSNOWRHO_1D(:) = 0. 03171 ZNDENT(:) = 0. 03172 ZNVIEU(:) = 0. 03173 ZSNOWAGE_1D(:) = 0. 03174 ! 03175 DO JJ=1,size(PSNOWRHO,1) 03176 IF(PSNOWDZ(JJ,1) == 0.0)THEN 03177 DO JST=2,INLVLS_USE(JJ) 03178 ZSNOWHEAT_1D(JJ) = ZSNOWHEAT_1D(JJ) + PSNOWDZ(JJ,JST)* & 03179 (SNOW3LSCAP(PSNOWRHO(JJ,JST))*(ZSNOWTEMP(JJ,JST)-XTT) & 03180 - XLMTT*PSNOWRHO(JJ,JST) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST) 03181 ZSNOW (JJ) = ZSNOW (JJ) + PSNOWDZ (JJ,JST) 03182 ZSNOWRHO_1D (JJ) = ZSNOWRHO_1D (JJ) + & 03183 PSNOWDZ (JJ,JST)*PSNOWRHO(JJ,JST) 03184 ZSNOWAGE_1D (JJ) = ZSNOWAGE_1D (JJ) + & 03185 PSNOWDZ (JJ,JST)*PSNOWRHO(JJ,JST)*PSNOWAGE(JJ,JST) 03186 ! snow grains 03187 IF(PSNOWGRAN1(JJ,JST)<-XEPSI)THEN ! Dendritic snow 03188 ZNDENT(JJ) = ZNDENT(JJ)+ 1.0 03189 ELSE ! Non dendritic snow 03190 ZNVIEU(JJ) = ZNVIEU(JJ)+1.0 03191 ENDIF 03192 ENDDO 03193 ENDIF 03194 END DO 03195 ZSNOWRHO_1D (:) = ZSNOWRHO_1D (:) / MAX(XSNOWDMIN,ZSNOW(:)) 03196 ZSNOWAGE_1D (:)= ZSNOWAGE_1D (:)/ MAX(XSNOWDMIN,ZSNOW(:)*ZSNOWRHO_1D (:)) 03197 ZSNOWRHO_1D (:) = MAX(XRHOSMIN_ES,MIN(XRHOSMAX_ES,ZSNOWRHO_1D(:))) 03198 ! 03199 ! Where uppermost snow layer has vanished, redistribute vertical 03200 ! snow mass and heat profiles (and associated quantities): 03201 ! 03202 CALL SNOW3LAVGRAIN(PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST, & 03203 ZSNOWGRAN1N, ZSNOWGRAN2N, ZSNOWHISTN,ZNDENT, ZNVIEU) 03204 ! 03205 DO JJ=1,SIZE(PSNOWRHO,1) 03206 IF(ZSNOW(JJ) /= 0.0) THEN 03207 PSNOWDZ(JJ,1:INLVLS_USE(JJ)) = ZSNOW(JJ)/INLVLS_USE(JJ) 03208 PSNOWHEAT(JJ,1:INLVLS_USE(JJ)) = ZSNOWHEAT_1D(JJ)/INLVLS_USE(JJ) 03209 PSNOWRHO(JJ,1:INLVLS_USE(JJ)) = ZSNOWRHO_1D(JJ) 03210 ! 03211 ZSCAP(JJ) = SNOW3LSCAP(ZSNOWRHO_1D(JJ)) 03212 PSNOWTEMP(JJ,1:INLVLS_USE(JJ)) = XTT + ( ((PSNOWHEAT(JJ,1:INLVLS_USE(JJ))/PSNOWDZ(JJ,1:INLVLS_USE(JJ))) & 03213 + XLMTT*PSNOWRHO(JJ,1:INLVLS_USE(JJ)))/ZSCAP(JJ) ) 03214 PSNOWLIQ(JJ,1:INLVLS_USE(JJ)) = MAX(0.0,PSNOWTEMP(JJ,1:INLVLS_USE(JJ))-XTT)*ZSCAP(JJ)* & 03215 PSNOWDZ(JJ,1:INLVLS_USE(JJ))/(XLMTT*XRHOLW) 03216 PSNOWTEMP(JJ,1:INLVLS_USE(JJ)) = MIN(XTT,PSNOWTEMP(JJ,1:INLVLS_USE(JJ))) 03217 ENDIF 03218 ENDDO 03219 IF (LHOOK) CALL DR_HOOK('SNOWCROEVAPGONE',1,ZHOOK_HANDLE) 03220 03221 ! 03222 END SUBROUTINE SNOWCROEVAPGONE 03223 ! 03224 ! 03225 ! 03226 03227 !#################################################################### 03228 !#################################################################### 03229 !#################################################################### 03230 SUBROUTINE SNOWNLFALL_UPGRID(TPTIME, OGLACIER, & 03231 PTSTEP,PSR,PTA,PVMOD, & 03232 PSNOW,PSNOWRHO,PSNOWDZ,PSNOWHEAT,PSNOWHMASS,PSNOWALB,PPERMSNOWFRAC, & 03233 PSNOWGRAN1,PSNOWGRAN2, GSNOWFALL,PSNOWDZN, & 03234 PSNOWRHOF,PSNOWDZF,PSNOWGRAN1F,PSNOWGRAN2F,PSNOWHISTF,PSNOWAGEF, & 03235 OMODIF_GRID,INLVLS_USE,OSNOWDRIFT) 03236 ! 03237 !! PURPOSE 03238 !! ------- 03239 ! Adds new snowfall and updates the vertical grid in order to keep an 03240 ! optimal discertisation 03241 ! 03242 !! AUTHOR 03243 !! ------ 03244 !! E. Brun * Meteo-France * 03245 !! 03246 ! 03247 USE MODD_CSTS, ONLY : XLMTT, XTT, XCI 03248 USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN, XANSMAX,XAGLAMAX,XSNOWCRITD 03249 USE MODD_SNOW_METAMO 03250 ! 03251 USE MODE_SNOW3L 03252 USE MODD_TYPE_DATE_SURF, ONLY: DATE_TIME 03253 ! 03254 IMPLICIT NONE 03255 ! 03256 !* 0.1 declarations of arguments 03257 ! 03258 TYPE(DATE_TIME), INTENT(IN) :: TPTIME ! current date and time 03259 LOGICAL, INTENT(IN) :: OGLACIER ! True = Over permanent snow and ice, 03260 ! initialise WGI=WSAT, 03261 ! Hsnow>=10m and allow 0.8<SNOALB<0.85 03262 ! False = No specific treatment 03263 ! 03264 REAL, INTENT(IN) :: PTSTEP 03265 ! 03266 REAL, DIMENSION(:), INTENT(IN) :: PSR, PTA, PVMOD, PPERMSNOWFRAC 03267 ! 03268 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOW, PSNOWALB 03269 ! 03270 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWDZ, PSNOWHEAT 03271 ! 03272 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWHMASS 03273 ! 03274 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWGRAN1, PSNOWGRAN2 03275 ! 03276 LOGICAL, DIMENSION(:), INTENT(INOUT):: GSNOWFALL 03277 ! 03278 ! Fresh snow characteristics 03279 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWRHOF, PSNOWDZF 03280 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWGRAN1F, PSNOWGRAN2F, PSNOWHISTF 03281 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWAGEF 03282 ! New vertical grid 03283 REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWDZN 03284 ! 03285 LOGICAL, DIMENSION(:), INTENT(OUT) ::OMODIF_GRID 03286 ! 03287 INTEGER, DIMENSION(:), INTENT(INOUT) :: INLVLS_USE 03288 03289 LOGICAL,INTENT(IN) :: OSNOWDRIFT ! if snowdrift then grain types are not modified by wind 03290 03291 !* 0.2 declarations of local variables 03292 ! 03293 INTEGER :: INLVLS, INLVLSMIN, INLVLSMAX, JJ, JST 03294 ! 03295 LOGICAL, DIMENSION(SIZE(PTA)) :: OAGREG_SURF 03296 ! 03297 REAL, DIMENSION(SIZE(PTA)) :: ZSNOWFALL,ZSNOWTEMP,ZSCAP,ZANSMAX 03298 REAL :: ZAGE_NOW 03299 ! 03300 ! ISBA-ES CROCUS (Pahaut 1976): snowfall density coefficients: 03301 ! 03302 REAL, PARAMETER :: ZSNOWFALL_A_SN = 109.0 ! kg/m3 03303 REAL, PARAMETER :: ZSNOWFALL_B_SN = 6.0 ! kg/(m3 K) 03304 REAL, PARAMETER :: ZSNOWFALL_C_SN = 26.0 ! kg/(m7/2 s1/2) 03305 ! 03306 ! Coefficients for the optimal vertical grid calculation 03307 REAL, PARAMETER :: DZ1=0.01 03308 REAL, PARAMETER :: DZ2=0.0125 03309 REAL, PARAMETER :: DZ3=0.015 03310 REAL, PARAMETER :: DZ3_bis=0.03 03311 REAL, PARAMETER :: DZ4=0.04 03312 REAL, PARAMETER :: DZ5=0.05 03313 REAL, PARAMETER :: DZ_base=0.02 03314 REAL, PARAMETER :: DZ_internal=0.07 03315 REAL, PARAMETER :: scale_cm=100. 03316 REAL,DIMENSION(5), PARAMETER :: DZMAX_internal=(/0.5,1.,2.,4.,10./) 03317 REAL, PARAMETER :: DZMIN_TOP_EXTREM=0.0001 03318 ! 03319 03320 ! Below this threshold of snowfall, new snowfall are aggregated with surface layer to avoid numerical problems 03321 ! (0.03 mm/h) 03322 REAL,PARAMETER :: XSNOWFALL_THRESHOLD=0.0333/3600. 03323 03324 ! The ratio between a new surface layer thickness and the second layer surface thickness is limited to 1/10 03325 REAL,PARAMETER :: XRATIO_NEWLAYER=0.1 03326 03327 ! Coefficients for cases with very thick snowpacks 03328 REAL, PARAMETER :: DEPTH_THRESHOLD1=3. 03329 REAL, PARAMETER :: DEPTH_THRESHOLD2=20. 03330 REAL, PARAMETER :: DEPTH_SURFACE=3. 03331 REAL :: PSNOW_UPPER ! snow depth treated normally (<= DEPTH_SURFACE) 03332 INTEGER :: NB_DEEP_LAYER,NB_UPPER_LAYER !separation between deep and upper layers 03333 ! if snow depth below DEPTH_SURFACE then NB_DEEP_LAYER=0 03334 REAL :: COEF_DEPTH !coefficient for repartition of deep snow above 3 meters 03335 INTEGER :: NB_MIN_LAYERS !why this test ? 03336 INTEGER :: NB_INTERMEDIATE ! number of intermediate layers (constant optimal gridding) 03337 INTEGER :: END_INTERMEDIATE ! layer indice for bottom of intermediate layers 03338 REAL :: THICKNESS_INTERMEDIATE 03339 INTEGER :: JSTDEEP,JSTEND 03340 ! 03341 ! Coefficients for computing the difference in 2 snow layer characteristics 03342 REAL, PARAMETER :: DIFF_1=20. 03343 REAL, PARAMETER :: DIFF_MAX=200. 03344 REAL, PARAMETER :: scale_DIFF=25. 03345 ! 03346 ! Coeefficients for snow layer splitting 03347 REAL, PARAMETER :: DZMIN_TOP=0.01 03348 REAL, PARAMETER :: DZMIN_TOP_BIS=0.005 03349 REAL, PARAMETER :: DZMIN_BOT=0.02 03350 REAL, PARAMETER :: SPLIT_COEF=8. 03351 ! 03352 ! Coeefficients for snow layer agregation 03353 REAL, PARAMETER :: AGREG_COEF_1=5. 03354 REAL, PARAMETER :: AGREG_COEF_2=4.5 03355 ! 03356 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) ::ZDZOPT 03357 03358 REAL :: ZPENALTY, ZDIFTYPE_INF,ZDIFTYPE_SUP,ZCRITSIZE, 03359 ZCRITSIZE_INF, ZCRITSIZE_SUP 03360 ! 03361 INTEGER JST_1, JJ_A_AGREG_SUP, JJ_A_AGREG_INF,JJ_A_DEDOUB 03362 03363 REAL :: ZSNOW2L 03364 03365 REAL(KIND=JPRB) :: ZHOOK_HANDLE 03366 ! 03367 !* 1.0 Initialization and snowage calculation for the present date 03368 ! 03369 IF (LHOOK) CALL DR_HOOK('SNOWNLFALL_UPGRID',0,ZHOOK_HANDLE) 03370 INLVLS = SIZE (PSNOWRHO(:,:),2) 03371 INLVLSMAX = SIZE (PSNOWRHO(:,:),2) 03372 INLVLSMIN = 3 03373 ZSNOWTEMP(:) = XTT 03374 GSNOWFALL(:) =.FALSE. 03375 ZSNOWFALL(:) = 0.0 !Matthieu Lafaysse 21/09/2012 03376 PSNOWHMASS(:) = 0.0 03377 PSNOWRHOF(:) = 0.0 03378 PSNOWDZF(:) = 0.0 03379 PSNOWGRAN1F(:) = 0.0 03380 PSNOWGRAN2F(:) = 0.0 03381 PSNOWHISTF(:) = 0.0 03382 PSNOWDZN(:,:) = PSNOWDZ(:,:) 03383 OAGREG_SURF(:)=.FALSE. 03384 OMODIF_GRID(:)=.FALSE. 03385 ! 03386 ! 03387 ! 03388 !* 1.1 Calculation of the optimal vertical grid size 03389 ! as a function of maximum number of layers and of current 03390 ! snow depth (modified 05/06/2012 by Matthieu Lafaysse) 03391 ! 03392 NB_MIN_LAYERS = 2 + INLVLSMAX/3 03393 ! 03394 DO JJ=1, SIZE(PSNOW(:)) 03395 03396 IF ((PSNOW(JJ) > DEPTH_THRESHOLD2).AND.( INLVLS_USE(JJ) > NB_MIN_LAYERS)) THEN 03397 03398 ! for very thick snowpack with enough snow layers 03399 ! special treatment 03400 ! we put the highest thickness in the lowest layers 03401 ! about 1/3 of layers for all snow except DEPTH_SURFACE=3 first meters 03402 03403 !number of "deep layers" 03404 NB_DEEP_LAYER= (INLVLS_USE(JJ)+INLVLSMAX)/6 03405 03406 !number of "upper layers" 03407 NB_UPPER_LAYER= INLVLS_USE(JJ)-NB_DEEP_LAYER 03408 !thickness of "upper layers" 03409 PSNOW_UPPER=DEPTH_SURFACE 03410 03411 !Arithmetic serie : 1+2+3+...+NB_DEEP_LAYER=NB_DEEP_LAYER*(NB_DEEP_LAYER+1)/2 03412 COEF_DEPTH = (PSNOW(JJ) - DEPTH_SURFACE)*2./((NB_DEEP_LAYER+1)*NB_DEEP_LAYER) 03413 03414 ! deep layers optimal thickness : 03415 ! increasing thickness with depth 03416 DO JSTDEEP=1,NB_DEEP_LAYER 03417 JST = NB_UPPER_LAYER + JSTDEEP 03418 ZDZOPT(JJ,JST) = COEF_DEPTH*JSTDEEP 03419 !This sum is equal to PSNOW(JJ)-DEPTH_SURFACE 03420 ENDDO 03421 03422 ELSE 03423 NB_UPPER_LAYER=INLVLS_USE(JJ) 03424 PSNOW_UPPER=PSNOW(JJ) 03425 END IF 03426 03427 ZDZOPT(JJ,1) = MIN(DZ1,PSNOW_UPPER/MAX(INLVLSMIN,NB_UPPER_LAYER)) 03428 ZDZOPT(JJ,2) = MIN(DZ2,PSNOW_UPPER/MAX(INLVLSMIN,NB_UPPER_LAYER)) 03429 ZDZOPT(JJ,3) = MIN(DZ3,PSNOW_UPPER/MAX(INLVLSMIN,NB_UPPER_LAYER)) 03430 ZDZOPT(JJ,4:5) = 0.0 ! initialization to avoid a problem in thickness_intermediate 03431 03432 IF (NB_UPPER_LAYER>0) THEN 03433 03434 !ML : replace > by >= on 12-12-20 because the last layer was not initialised in case of thick snowpacks 03435 IF (NB_UPPER_LAYER>=3) ZDZOPT(JJ,3) = MIN(DZ3_bis,PSNOW_UPPER/NB_UPPER_LAYER) 03436 IF (NB_UPPER_LAYER>=4) ZDZOPT(JJ,4) = MIN(DZ4,PSNOW_UPPER/NB_UPPER_LAYER) 03437 IF (NB_UPPER_LAYER>=5) ZDZOPT(JJ,5) = MIN(DZ5,PSNOW_UPPER/NB_UPPER_LAYER) 03438 03439 IF (NB_UPPER_LAYER==INLVLS_USE(JJ)) THEN 03440 ! last layer of' upper layers' : normal case : thin layer 03441 ZDZOPT(JJ,NB_UPPER_LAYER)= & 03442 MIN(DZ_base,PSNOW_UPPER/MAX(INLVLSMIN,NB_UPPER_LAYER)) 03443 03444 !remaining snow for remaining layers 03445 THICKNESS_INTERMEDIATE=PSNOW_UPPER - SUM(ZDZOPT(JJ,1:5))-ZDZOPT(JJ,NB_UPPER_LAYER) 03446 03447 IF (PSNOW_UPPER<DEPTH_THRESHOLD1) THEN 03448 NB_INTERMEDIATE=NB_UPPER_LAYER-6 03449 END_INTERMEDIATE=NB_UPPER_LAYER-1 03450 ELSE 03451 NB_INTERMEDIATE=NB_UPPER_LAYER-8 03452 END_INTERMEDIATE=NB_UPPER_LAYER-3 03453 IF (NB_INTERMEDIATE>0) THEN 03454 THICKNESS_INTERMEDIATE=THICKNESS_INTERMEDIATE*NB_INTERMEDIATE/FLOAT(NB_INTERMEDIATE+1) 03455 END IF 03456 END IF 03457 ELSE 03458 ! case with very thick snowpacks : 03459 ! the last layer of upper layers is not an exception 03460 THICKNESS_INTERMEDIATE=PSNOW_UPPER - SUM(ZDZOPT(JJ,1:5)) 03461 NB_INTERMEDIATE=NB_UPPER_LAYER-5 03462 END_INTERMEDIATE=NB_UPPER_LAYER 03463 END IF 03464 03465 ! For thick snowpack : add maximum value of optimal thickness to avoid too 03466 ! large differencies between layers 03467 IF (NB_INTERMEDIATE>0) THEN 03468 JSTEND=MIN(END_INTERMEDIATE,10) 03469 DO JST=6,JSTEND 03470 ZDZOPT(JJ,JST) = MIN(DZMAX_internal(JST-5),& 03471 MAX(DZ_internal,THICKNESS_INTERMEDIATE/NB_INTERMEDIATE)) 03472 END DO 03473 03474 IF (END_INTERMEDIATE>10) THEN 03475 DO JST=11,END_INTERMEDIATE 03476 ZDZOPT(JJ,JST) = MAX(DZ_internal,THICKNESS_INTERMEDIATE/NB_INTERMEDIATE) 03477 END DO 03478 END IF 03479 END IF 03480 03481 IF ((PSNOW_UPPER>DEPTH_THRESHOLD1).AND.(NB_UPPER_LAYER>=8)) THEN 03482 03483 !Linear interpolation of optimal thickness between layers N-3 and N : 03484 ZDZOPT(JJ,NB_UPPER_LAYER-2) = 0.34*ZDZOPT(JJ,NB_UPPER_LAYER)+& 03485 0.66*ZDZOPT(JJ,NB_UPPER_LAYER-3) 03486 ZDZOPT(JJ,NB_UPPER_LAYER-1) = 0.66*ZDZOPT(JJ,NB_UPPER_LAYER)+& 03487 0.34*ZDZOPT(JJ,NB_UPPER_LAYER-3) 03488 ENDIF 03489 END IF 03490 END DO 03491 03492 03493 !************************************************************************************ 03494 !This was the initial code for optimal layers until may 2012 03495 03496 03497 ! ! ! ! ! 03498 ! ! ! ! ! !* 1.1 Calculation of the optimal vertical grid size 03499 ! ! ! ! ! ! as a function of maximum number of layers and of current 03500 ! ! ! ! ! ! snow depth 03501 ! ! ! ! ! ! 03502 ! ! ! ! ! DO JJ=1, SIZE(PSNOW(:)) 03503 ! ! ! ! ! ZDZOPT(JJ,1) = MIN(DZ1,PSNOW(JJ)/MAX(INLVLSMIN,INLVLS_USE(JJ))) 03504 ! ! ! ! ! ZDZOPT(JJ,2) = MIN(DZ2,PSNOW(JJ)/MAX(INLVLSMIN,INLVLS_USE(JJ))) 03505 ! ! ! ! ! ZDZOPT(JJ,3) = MIN(DZ3,PSNOW(JJ)/MAX(INLVLSMIN,INLVLS_USE(JJ))) 03506 ! ! ! ! ! IF (INLVLS_USE(JJ)>3) ZDZOPT(JJ,3) = MIN(DZ3_bis,PSNOW(JJ)/INLVLS_USE(JJ)) 03507 ! ! ! ! ! IF (INLVLS_USE(JJ)>4) ZDZOPT(JJ,4) = MIN(DZ4,PSNOW(JJ)/INLVLS_USE(JJ)) 03508 ! ! ! ! ! IF (INLVLS_USE(JJ)>5) ZDZOPT(JJ,5) = MIN(DZ5,PSNOW(JJ)/INLVLS_USE(JJ)) 03509 ! ! ! ! ! IF (INLVLS_USE(JJ)>0) ZDZOPT(JJ,INLVLS_USE(JJ))= & 03510 ! ! ! ! ! MIN(DZ_base,PSNOW(JJ)/MAX(INLVLSMIN,INLVLS_USE(JJ))) 03511 ! ! ! ! ! DO JST=6,INLVLS_USE(JJ)-1,1 03512 ! ! ! ! ! ZDZOPT(JJ,JST) = MAX(DZ_internal,(PSNOW(JJ) - SUM(ZDZOPT(JJ,1:5))- & 03513 ! ! ! ! ! ZDZOPT(JJ,INLVLS_USE(JJ))) /(INLVLS_USE(JJ)-6)) 03514 ! ! ! ! ! END DO 03515 ! ! ! ! ! END DO 03516 ! ! ! ! ! ! 03517 ! ! ! ! ! 03518 03519 !************************************************************************************ 03520 03521 !* 2.0 Fresh snow characteristics 03522 ! 03523 ! 03524 ! 03525 ! Heat content of newly fallen snow (J/m2): 03526 ! NOTE for now we assume the snowfall has 03527 ! the temperature of the snow surface upon reaching the snow. 03528 ! This is done as opposed to using the air temperature since 03529 ! this flux is quite small and has little to no impact 03530 ! on the time scales of interest. If we use the above assumption 03531 ! then, then the snowfall advective heat flux is zero. 03532 !! 03533 DO JJ=1, SIZE(PSNOW(:)) 03534 IF (PSR(JJ) > UEPSI ) THEN 03535 ! newly fallen snow characteristics: 03536 IF (INLVLS_USE(JJ)>0) THEN !Case of new snowfall on a previously snow-free surface 03537 ZSCAP(JJ) = XCI*PSNOWRHO(JJ,1) 03538 ZSNOWTEMP(JJ) = XTT + (PSNOWHEAT(JJ,1) + & 03539 XLMTT*PSNOWRHO(JJ,1)*PSNOWDZ(JJ,1))/ & 03540 (ZSCAP(JJ)*MAX(XSNOWDMIN/INLVLS,PSNOWDZ(JJ,1))) 03541 ELSE ! case with bare ground 03542 ZSNOWTEMP(JJ)=PTA(JJ) 03543 ENDIF 03544 ZSNOWTEMP(JJ) = MIN(XTT, ZSNOWTEMP(JJ)) 03545 ! 03546 PSNOWHMASS(JJ) = PSR(JJ)*(XCI*(ZSNOWTEMP(JJ)-XTT)-XLMTT)*PTSTEP 03547 PSNOWRHOF(JJ) = MAX(XRHOSMIN_ES, ZSNOWFALL_A_SN + & 03548 ZSNOWFALL_B_SN*(PTA(JJ)-XTT)+ & 03549 ZSNOWFALL_C_SN*MIN(PVMOD(JJ),SQRT(PVMOD(JJ)))) 03550 ZSNOWFALL(JJ) = PSR(JJ)*PTSTEP/PSNOWRHOF(JJ) ! snowfall thickness (m) 03551 PSNOW(JJ) = PSNOW(JJ) + ZSNOWFALL(JJ) 03552 PSNOWDZF(JJ) = ZSNOWFALL(JJ) 03553 IF (OSNOWDRIFT) THEN 03554 PSNOWGRAN1F(JJ)=-XGRAN 03555 PSNOWGRAN2F(JJ)=XNSPH3 03556 ELSE 03557 PSNOWGRAN1F(JJ)= MAX(MIN(XNDEN1*PVMOD(JJ)-XNDEN2,XNDEN3),-XGRAN) 03558 PSNOWGRAN2F(JJ)= MIN(MAX(XNSPH1*PVMOD(JJ)+XNSPH2,XNSPH3),XNSPH4) 03559 END IF 03560 PSNOWHISTF(JJ) = 0.0 03561 PSNOWAGEF(JJ) = 0. 03562 GSNOWFALL(JJ) =.TRUE. 03563 OMODIF_GRID(JJ)=.TRUE. 03564 ! 03565 ENDIF 03566 ENDDO 03567 ! 03568 ! intialize the albedo: 03569 ! penser a changer 0.000001 par UEPSI 03570 IF(OGLACIER)THEN 03571 ZANSMAX(:) = XAGLAMAX * PPERMSNOWFRAC(:) + XANSMAX * (1.0-PPERMSNOWFRAC(:)) 03572 ELSE 03573 ZANSMAX(:) = XANSMAX 03574 ENDIF 03575 WHERE(GSNOWFALL(:) .AND. ABS(PSNOW(:)-ZSNOWFALL(:))< 0.000001) 03576 PSNOWALB(:) = ZANSMAX(:) 03577 END WHERE 03578 ! 03579 ! Computation of the new grid size 03580 ! It starts with successive exclusive cases 03581 ! Each case is described inside the corresponding condition 03582 ! 03583 ! cases with fresh snow 03584 ! 03585 DO JJ=1,SIZE(PSNOW(:)) ! grid point loop 03586 IF((.NOT.GSNOWFALL(JJ)).AND.PSNOW(JJ)>=XSNOWCRITD.AND. & 03587 INLVLS_USE(JJ)>=INLVLSMIN) THEN 03588 ! no fresh snow + deep enough snowpack + enough snow layers ==> no change 03589 ! 03590 ELSEIF(PSNOW(JJ)<XSNOWCRITD.OR.PSNOW(JJ)==ZSNOWFALL(JJ).OR. & 03591 INLVLS_USE(JJ)<INLVLSMIN ) THEN 03592 ! too shallow snowpack or only fresh snow or too few layers 03593 ! ==> uniform grid and identical snow layers / number depends on snow depth 03594 OMODIF_GRID(JJ)=.TRUE. 03595 INLVLS_USE(JJ)= MAX(INLVLSMIN,MIN(INLVLSMAX,INT(PSNOW(JJ)*scale_cm))) 03596 PSNOWDZN(JJ,1:INLVLS_USE(JJ))=PSNOW(JJ)/INLVLS_USE(JJ) 03597 ELSE 03598 ! fresh snow over snow covered ground + enough snow layers 03599 OMODIF_GRID(JJ)=.TRUE. 03600 ZDIFTYPE_SUP= SNOW3LDIFTYP(PSNOWGRAN1(JJ,1),PSNOWGRAN1F(JJ),& 03601 PSNOWGRAN2(JJ,1),PSNOWGRAN2F(JJ)) 03602 IF ((ZDIFTYPE_SUP<DIFF_1.AND.PSNOWDZ(JJ,1)< ZDZOPT(JJ,1)).OR. & 03603 (PSR(JJ)<XSNOWFALL_THRESHOLD .AND.PSNOWDZ(JJ,1) < 2.* ZDZOPT(JJ,1)) & 03604 .OR. PSNOWDZ(JJ,1)< DZMIN_TOP_EXTREM ) THEN 03605 ! Fresh snow is similar to a shallow surface layer (< ZDZOPT) 03606 ! or snowfall is very low and the surface layer not too deep (< 2*ZDZOPT) [NEW CONDITION 11/2012] 03607 ! or the surface layer is extremely thin (< DZMIN_TOP_EXTREM) [NEW CONDITION 11/2012] 03608 ! The two new conditions are necessary for forcings with very low precipitation 03609 ! (e.g. ERA interim reanalyses, or climate models) 03610 03611 ! ==> fresh snow is agregated to the surface layer 03612 PSNOWDZN(JJ,1)=PSNOWDZ(JJ,1)+PSNOWDZF(JJ) 03613 DO JST= INLVLS_USE(JJ), 2,-1 03614 PSNOWDZN(JJ, JST)=PSNOWDZ(JJ,JST) 03615 ENDDO 03616 ELSEIF(INLVLS_USE(JJ)<INLVLSMAX) THEN 03617 ! fresh snow is too different from the surface or the surface is too deep 03618 ! and there is room for extra layers ==> we create a new layer 03619 INLVLS_USE(JJ)=INLVLS_USE(JJ)+1 03620 03621 IF(PSNOWDZF(JJ) > XRATIO_NEWLAYER * PSNOWDZ(JJ,2)) THEN 03622 ! Snowfall is sufficient to create a new layer not lower than 1/10 of the second layer 03623 PSNOWDZN(JJ,1)=PSNOWDZF(JJ) 03624 DO JST= INLVLS_USE(JJ), 2,-1 03625 PSNOWDZN(JJ, JST)=PSNOWDZ(JJ,JST-1) 03626 ENDDO 03627 ELSE 03628 ! The ratio would be lower than 1/10 : [NEW : 11/2012] 03629 ! aggregate a part of the old layer with fresh snow to limit the ratio to 1/10. 03630 ZSNOW2L=PSNOWDZF(JJ)+PSNOWDZ(JJ,1) 03631 PSNOWDZN(JJ,1)=XRATIO_NEWLAYER*ZSNOW2L 03632 PSNOWDZN(JJ,2)=(1.-XRATIO_NEWLAYER)*ZSNOW2L 03633 DO JST= INLVLS_USE(JJ), 3,-1 03634 PSNOWDZN(JJ, JST)=PSNOWDZ(JJ,JST-1) 03635 ENDDO 03636 ENDIF 03637 03638 03639 ELSE 03640 ! fresh snow is too different from the surface or the surface is too deep 03641 ! and there is no room for extra layers 03642 ! ==> we agregate internal most similar snowlayers and create a new surface layer 03643 JJ_A_AGREG_SUP=1 03644 JJ_A_AGREG_INF=2 03645 ZDIFTYPE_INF= SNOW3LDIFTYP(PSNOWGRAN1(JJ,1),PSNOWGRAN1F(JJ),& 03646 PSNOWGRAN2(JJ,1),PSNOWGRAN2F(JJ)) 03647 ZCRITSIZE_INF= scale_DIFF*(PSNOWDZ(JJ,1)/ZDZOPT(JJ,1)+ & 03648 PSNOWDZ(JJ,2)/ZDZOPT(JJ,2)) 03649 ZPENALTY=ZDIFTYPE_INF+ZCRITSIZE_INF 03650 DO JST=2, INLVLS_USE(JJ)-1 03651 ZDIFTYPE_SUP= SNOW3LDIFTYP(PSNOWGRAN1(JJ,JST-1),PSNOWGRAN1(JJ, JST),& 03652 PSNOWGRAN2(JJ,JST-1),PSNOWGRAN2(JJ, JST)) 03653 ZDIFTYPE_INF= SNOW3LDIFTYP(PSNOWGRAN1(JJ,JST+1),PSNOWGRAN1(JJ, JST),& 03654 PSNOWGRAN2(JJ,JST+1),PSNOWGRAN2(JJ, JST)) 03655 ZCRITSIZE_SUP= scale_DIFF*(PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST)+ & 03656 PSNOWDZ(JJ,JST-1)/ZDZOPT(JJ,JST-1)) 03657 ZCRITSIZE_INF= scale_DIFF*(PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST)+ & 03658 PSNOWDZ(JJ,JST+1)/ZDZOPT(JJ,JST+1)) 03659 IF(ZDIFTYPE_SUP+ZCRITSIZE_SUP < ZPENALTY) THEN 03660 ZPENALTY=ZDIFTYPE_SUP+ZCRITSIZE_SUP 03661 JJ_A_AGREG_SUP=JST-1 03662 JJ_A_AGREG_INF=JST 03663 ENDIF 03664 IF(ZDIFTYPE_INF+ZCRITSIZE_INF < ZPENALTY) THEN 03665 ZPENALTY=ZDIFTYPE_INF+ZCRITSIZE_INF 03666 JJ_A_AGREG_SUP=JST 03667 JJ_A_AGREG_INF=JST+1 03668 ENDIF 03669 ENDDO 03670 ZDIFTYPE_SUP= SNOW3LDIFTYP(PSNOWGRAN1(JJ,INLVLS_USE(JJ)-1), & 03671 PSNOWGRAN1(JJ,INLVLS_USE(JJ)), PSNOWGRAN2(JJ,INLVLS_USE(JJ)-1), & 03672 PSNOWGRAN2(JJ, INLVLS_USE(JJ))) 03673 ZCRITSIZE_SUP=scale_DIFF*(PSNOWDZ(JJ,INLVLS_USE(JJ))/ZDZOPT(JJ,INLVLS_USE(JJ))+ & 03674 PSNOWDZ(JJ,INLVLS_USE(JJ)-1)/ZDZOPT(JJ,INLVLS_USE(JJ)-1)) 03675 IF(ZDIFTYPE_SUP+ZCRITSIZE_SUP < ZPENALTY) THEN 03676 JJ_A_AGREG_SUP=INLVLS_USE(JJ)-1 03677 JJ_A_AGREG_INF=INLVLS_USE(JJ) 03678 ENDIF 03679 ! agregation of the similar layers and shift of upper layers 03680 PSNOWDZN(JJ,JJ_A_AGREG_INF)=PSNOWDZ(JJ,JJ_A_AGREG_INF) + & 03681 PSNOWDZ(JJ,JJ_A_AGREG_SUP) 03682 DO JST= JJ_A_AGREG_SUP, 2,-1 03683 PSNOWDZN(JJ, JST)=PSNOWDZ(JJ,JST-1) 03684 ENDDO 03685 PSNOWDZN(JJ,1)=PSNOWDZF(JJ) 03686 03687 ! Limit the ratio between the new layer and the one beneath (ratio 1/10) 03688 ! [NEW : 11/2012] 03689 IF(PSNOWDZN(JJ,1) < XRATIO_NEWLAYER * PSNOWDZN(JJ,2)) THEN 03690 ZSNOW2L=PSNOWDZN(JJ,1) + PSNOWDZN(JJ,2) 03691 PSNOWDZN(JJ,1) = XRATIO_NEWLAYER*ZSNOW2L 03692 PSNOWDZN(JJ,2)= (1.-XRATIO_NEWLAYER)* ZSNOW2L 03693 ENDIF 03694 03695 ENDIF 03696 03697 ENDIF ! end of the case with fresh snow 03698 ENDDO ! end loop grid points 03699 ! 03700 ! cases with no fresh snow and no previous grid resize 03701 ! 03702 IF(INLVLSMIN==INLVLSMAX) THEN ! specific case with INLSVSMIN = INLVLSMAX (INLVLS) 03703 ! check if surface layer depth is too small 03704 ! in such a case looks for an other layer to be split 03705 DO JJ=1,SIZE(PSNOW(:)) ! loop grid points 03706 IF(.NOT.OMODIF_GRID(JJ)) THEN 03707 IF(PSNOWDZ(JJ,1)< DZMIN_TOP) THEN 03708 OMODIF_GRID(JJ)=.TRUE. 03709 OAGREG_SURF(JJ)=.TRUE. 03710 ZPENALTY=PSNOWDZ(JJ,2)/ZDZOPT(JJ,2) 03711 IF(PSNOWDZ(JJ,2)<DZMIN_TOP) ZPENALTY=0. 03712 JJ_A_DEDOUB=2 03713 DO JST=3, INLVLS 03714 ZCRITSIZE=PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST) 03715 IF(JST==INLVLS.AND.PSNOWDZ(JJ,JST)<DZMIN_BOT) ZCRITSIZE=0. 03716 IF (ZCRITSIZE > ZPENALTY) THEN 03717 ZPENALTY=ZCRITSIZE 03718 JJ_A_DEDOUB=JST 03719 ENDIF 03720 ENDDO 03721 IF(JJ_A_DEDOUB == 2) THEN 03722 PSNOWDZN(JJ,1)=0.5*(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 03723 PSNOWDZN(JJ,2)=0.5*(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 03724 ELSE 03725 PSNOWDZN(JJ,1)=(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 03726 DO JST=2, JJ_A_DEDOUB-2 03727 PSNOWDZN(JJ,JST)=PSNOWDZ(JJ,JST+1) 03728 ENDDO 03729 PSNOWDZN(JJ,JJ_A_DEDOUB-1)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03730 PSNOWDZN(JJ,JJ_A_DEDOUB)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03731 ENDIF 03732 ENDIF 03733 ENDIF 03734 ENDDO ! end grid points loop 03735 ! 03736 ! check if bottom layer depth is too small 03737 ! in such a case agregation with upper layer and 03738 ! looks for an other layer to be splitted 03739 DO JJ=1,SIZE(PSNOW(:)) ! loop grid points 03740 IF(.NOT.OMODIF_GRID(JJ)) THEN 03741 IF(PSNOWDZ(JJ,INLVLS)< DZMIN_TOP) THEN 03742 OMODIF_GRID(JJ)=.TRUE. 03743 ZPENALTY=PSNOWDZ(JJ,INLVLS-2)/ZDZOPT(JJ,INLVLS-2) 03744 JJ_A_DEDOUB=INLVLS-2 03745 DO JST= MAX(1,INLVLS-3),1,-1 03746 ZCRITSIZE=PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST) 03747 IF(JST==1.AND.PSNOWDZ(JJ,JST)<DZMIN_BOT) ZCRITSIZE=0. 03748 IF (ZCRITSIZE > ZPENALTY) THEN 03749 ZPENALTY=ZCRITSIZE 03750 JJ_A_DEDOUB=JST 03751 ENDIF 03752 ENDDO 03753 IF(JJ_A_DEDOUB == INLVLS-1) THEN 03754 PSNOWDZN(JJ,INLVLS)=0.5*(PSNOWDZ(JJ,INLVLS)+PSNOWDZ(JJ,INLVLS-1)) 03755 PSNOWDZN(JJ,INLVLS-1)=0.5*(PSNOWDZ(JJ,INLVLS)+PSNOWDZ(JJ,INLVLS-1)) 03756 ! write( *,*) '2 couches basales agregees pour creer 2 couches identiques' 03757 ELSE 03758 ! write( *,*) 'couches basales agregees et on dedouble ',JJ_A_DEDOUB 03759 PSNOWDZN(JJ,INLVLS)=(PSNOWDZ(JJ,INLVLS)+PSNOWDZ(JJ,INLVLS-1)) 03760 DO JST=INLVLS-1, JJ_A_DEDOUB+1,-1 03761 PSNOWDZN(JJ,JST)=PSNOWDZ(JJ,JST-1) 03762 ENDDO 03763 PSNOWDZN(JJ,JJ_A_DEDOUB)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03764 PSNOWDZN(JJ,JJ_A_DEDOUB+1)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03765 ENDIF 03766 ENDIF 03767 ENDIF 03768 ENDDO ! end grid points loop 03769 ENDIF ! end specific case INLSVSMIN = INLVLSMAX 03770 ! 03771 ! case without new snowfall and INVLS > INLVLSMIN 03772 ! 03773 ! check if surface layer depth is too small 03774 ! in such a case agregation with layer beneath 03775 ! in case of reaching INLVLSMIN, looks for an other layer to be splitted 03776 DO JJ=1,SIZE(PSNOW(:)) 03777 IF((.NOT.GSNOWFALL(JJ)).AND.PSNOW(JJ)> XSNOWCRITD.AND. & 03778 (.NOT.OMODIF_GRID(JJ))) THEN 03779 IF(PSNOWDZ(JJ,1)< DZMIN_TOP_BIS) THEN ! case shallow surface layer 03780 OMODIF_GRID(JJ)=.TRUE. 03781 OAGREG_SURF(JJ)=.TRUE. 03782 IF(INLVLS_USE(JJ)>INLVLSMIN) THEN ! case minimum not reached 03783 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 03784 PSNOWDZN(JJ,1)=PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2) 03785 DO JST=2,INLVLS_USE(JJ) 03786 PSNOWDZN(JJ,JST)=PSNOWDZ(JJ,JST+1) 03787 ENDDO 03788 ELSE ! case minimum reached 03789 ZPENALTY=PSNOWDZ(JJ,2)/ZDZOPT(JJ,2) 03790 IF(PSNOWDZ(JJ,2)<DZMIN_TOP) ZPENALTY=0. 03791 JJ_A_DEDOUB=2 03792 DO JST=3, INLVLS_USE(JJ) 03793 ZCRITSIZE=PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST) 03794 IF(JST==INLVLS_USE(JJ).AND.PSNOWDZ(JJ,JST)<DZMIN_BOT) ZCRITSIZE=0. 03795 IF (ZCRITSIZE > ZPENALTY) THEN 03796 ZPENALTY=ZCRITSIZE 03797 JJ_A_DEDOUB=JST 03798 ENDIF 03799 ENDDO 03800 IF(JJ_A_DEDOUB == 2) THEN ! case splitted layer == 2 03801 PSNOWDZN(JJ,1)=0.5*(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 03802 PSNOWDZN(JJ,2)=0.5*(PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2)) 03803 ELSE ! case splitted layer =/ 2 03804 PSNOWDZN(JJ,1)=PSNOWDZ(JJ,1)+PSNOWDZ(JJ,2) 03805 DO JST=2, JJ_A_DEDOUB-2 03806 PSNOWDZN(JJ,JST)=PSNOWDZ(JJ,JST+1) 03807 ENDDO 03808 PSNOWDZN(JJ,JJ_A_DEDOUB-1)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03809 PSNOWDZN(JJ,JJ_A_DEDOUB)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03810 ENDIF ! end case splitted layer =/ 2 03811 ENDIF ! end case minimum reached 03812 ENDIF ! end case shallow surface layer 03813 ENDIF 03814 ENDDO ! end grid points loop 03815 ! 03816 ! check if bottom layer depth is too small 03817 ! in such a case agregation with above layer 03818 ! in case of reaching INLVLSMIN, looks for an other layer to be splitted 03819 DO JJ=1,SIZE(PSNOW(:)) 03820 IF((.NOT.GSNOWFALL(JJ)).AND.(.NOT.OAGREG_SURF(JJ)).AND. & 03821 (.NOT.OMODIF_GRID(JJ)).AND. PSNOW(JJ)> XSNOWCRITD) THEN 03822 IF(PSNOWDZ(JJ,INLVLS_USE(JJ))< DZMIN_TOP) THEN ! case shallow bottom layer 03823 OMODIF_GRID(JJ)=.TRUE. 03824 IF(INLVLS_USE(JJ)>INLVLSMIN) THEN ! case minimum not reached 03825 OMODIF_GRID(JJ)=.TRUE. 03826 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 03827 PSNOWDZN(JJ,INLVLS_USE(JJ))=PSNOWDZ(JJ,INLVLS_USE(JJ))+ & 03828 PSNOWDZ(JJ,INLVLS_USE(JJ)+1) 03829 ELSE ! case minimum reached 03830 ZPENALTY=PSNOWDZ(JJ,INLVLS_USE(JJ)-2)/ZDZOPT(JJ,INLVLS_USE(JJ)-2) 03831 JJ_A_DEDOUB=INLVLS_USE(JJ)-2 03832 DO JST= MAX(1,INLVLS_USE(JJ)-3),1,-1 03833 ZCRITSIZE=PSNOWDZ(JJ,JST)/ZDZOPT(JJ,JST) 03834 IF(JST==1.AND.PSNOWDZ(JJ,JST)<DZMIN_BOT) ZCRITSIZE=0. 03835 IF (ZCRITSIZE > ZPENALTY) THEN 03836 ZPENALTY=ZCRITSIZE 03837 JJ_A_DEDOUB=JST 03838 ENDIF 03839 ENDDO 03840 IF(JJ_A_DEDOUB == INLVLS_USE(JJ)-1) THEN 03841 PSNOWDZN(JJ,INLVLS_USE(JJ))=0.5*(PSNOWDZ(JJ,INLVLS_USE(JJ))+ & 03842 PSNOWDZ(JJ,INLVLS_USE(JJ)-1)) 03843 PSNOWDZN(JJ,INLVLS_USE(JJ)-1)=0.5*(PSNOWDZ(JJ,INLVLS_USE(JJ))+ & 03844 PSNOWDZ(JJ,INLVLS_USE(JJ)-1)) 03845 ELSE 03846 PSNOWDZN(JJ,INLVLS_USE(JJ))=(PSNOWDZ(JJ,INLVLS_USE(JJ))+ & 03847 PSNOWDZ(JJ,INLVLS_USE(JJ)-1)) 03848 DO JST=INLVLS_USE(JJ)-1, JJ_A_DEDOUB+1,-1 03849 PSNOWDZN(JJ,JST)=PSNOWDZ(JJ,JST-1) 03850 ENDDO 03851 PSNOWDZN(JJ,JJ_A_DEDOUB)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03852 PSNOWDZN(JJ,JJ_A_DEDOUB+1)=0.5*PSNOWDZ(JJ,JJ_A_DEDOUB) 03853 ENDIF ! end case splitted layer =/ 2 03854 ENDIF ! end case minimum reached 03855 ENDIF ! end case shallow surface layer 03856 ENDIF 03857 ENDDO ! end grid points loop 03858 ! 03859 ! case whithout new snow fall and without a previous grid resize 03860 ! looks for a shallow layer to be splitted according to its depth and to 03861 ! the optimal grid size 03862 DO JJ=1,SIZE(PSNOW(:)) 03863 IF (.NOT.OMODIF_GRID(JJ))THEN 03864 DO JST=1,INLVLS-3 03865 IF (JST<=INLVLS_USE(JJ).AND.INLVLS_USE(JJ)<INLVLS-3.AND. & 03866 (.NOT.OMODIF_GRID(JJ))) THEN 03867 IF(PSNOWDZ(JJ,JST)>(SPLIT_COEF-FLOAT(INLVLS-INLVLS_USE(JJ))/& 03868 MAX(1.,FLOAT(INLVLS-INLVLSMIN)))*ZDZOPT(JJ,JST)) THEN 03869 DO JST_1=INLVLS_USE(JJ)+1, JST+2,-1 03870 PSNOWDZN(JJ,JST_1)=PSNOWDZ(JJ,JST_1-1) 03871 ZDZOPT(JJ,JST_1)=ZDZOPT(JJ,JST_1-1) 03872 ENDDO 03873 ! generale case : old layer divided in two equal layers 03874 IF(JST/=1.OR.PSNOWDZ(JJ,JST) < 3.*ZDZOPT(JJ,1))THEN 03875 PSNOWDZN(JJ,JST+1)=0.5*PSNOWDZ(JJ,JST) 03876 PSNOWDZN(JJ,JST)=0.5*PSNOWDZ(JJ,JST) 03877 ELSE 03878 ! if thick surface layer : force the surface layer to this value to avoid successive resizing 03879 ! [NEW : 11/2012] 03880 PSNOWDZN(JJ,1)= 1.5*ZDZOPT(JJ,1) 03881 PSNOWDZN(JJ,2)= PSNOWDZ(JJ,JST)-PSNOWDZN(JJ,1) 03882 ENDIF 03883 INLVLS_USE(JJ)=INLVLS_USE(JJ)+1 03884 OMODIF_GRID(JJ)=.TRUE. 03885 ENDIF 03886 03887 ENDIF 03888 ENDDO 03889 ENDIF 03890 ENDDO 03891 ! 03892 ! case whithout new snow fall and without a previous grid resize 03893 ! looks for a deep layer to be agregated to the layer beneath if similar 03894 ! according to its depth and to the optimal grid size 03895 03896 !NB : allow these changes for 5 layers and more [NEW] (before : 6 layers) 03897 03898 DO JJ=1,SIZE(PSNOW(:)) 03899 IF (.NOT.OMODIF_GRID(JJ))THEN 03900 DO JST=2,INLVLS 03901 IF (JST<=INLVLS_USE(JJ)-1.AND.INLVLS_USE(JJ)>INLVLSMIN+1.AND. & 03902 (.NOT.OMODIF_GRID(JJ))) THEN 03903 ZDIFTYPE_INF= SNOW3LDIFTYP(PSNOWGRAN1(JJ,JST+1),PSNOWGRAN1(JJ, JST),& 03904 PSNOWGRAN2(JJ,JST+1),PSNOWGRAN2(JJ, JST)) 03905 ZDIFTYPE_INF=MAX(DIFF_1,MIN(DIFF_MAX,ZDIFTYPE_INF)) 03906 IF(PSNOWDZ(JJ,JST)<ZDZOPT(JJ,JST)*AGREG_COEF_1/ZDIFTYPE_INF & 03907 .AND.PSNOWDZ(JJ,JST)+PSNOWDZ(JJ,JST+1)<AGREG_COEF_2*& 03908 MAX(ZDZOPT(JJ,JST),ZDZOPT(JJ,JST+1))) THEN 03909 PSNOWDZN(JJ,JST)= PSNOWDZ(JJ,JST)+PSNOWDZ(JJ,JST+1) 03910 ZDZOPT(JJ,JST)=ZDZOPT(JJ,JST+1) 03911 DO JST_1=JST+1,INLVLS_USE(JJ)-1 03912 PSNOWDZN(JJ,JST_1)=PSNOWDZ(JJ,JST_1+1) 03913 ZDZOPT(JJ,JST_1)=ZDZOPT(JJ,JST_1+1) 03914 ENDDO 03915 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 03916 OMODIF_GRID(JJ)=.TRUE. 03917 ENDIF 03918 ENDIF 03919 ENDDO 03920 ENDIF 03921 ENDDO 03922 03923 ! [NEW : 11/2012] 03924 ! In case of very low snow fall checks if a new internal snow layer is too shallow 03925 ! even if a the grid has already been resized in this time step 03926 ! starts from bottom to INLVS_USE-3 until old and new grid differ 03927 03928 DO JJ=1,SIZE(PSNOW(:)) 03929 IF (.NOT. GSNOWFALL(JJ).OR.INLVLS_USE(JJ)< INLVLSMIN+3) CYCLE ! go to next point 03930 IF(ABS(PSNOWDZN(JJ,INLVLS_USE(JJ))-PSNOWDZ(JJ,INLVLS_USE(JJ))) > UEPSI) & 03931 CYCLE ! go to next point 03932 ! bottom layer 03933 IF(PSNOWDZN(JJ,INLVLS_USE(JJ))< DZMIN_TOP) THEN ! case shallow bottom layer 03934 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 03935 PSNOWDZN(JJ,INLVLS_USE(JJ))=PSNOWDZN(JJ,INLVLS_USE(JJ))+ & 03936 PSNOWDZN(JJ,INLVLS_USE(JJ)+1) 03937 ELSE 03938 ! internal layer 03939 DO JST=INLVLS_USE(JJ)-1,4 03940 IF(ABS(PSNOWDZN(JJ,JST)-PSNOWDZ(JJ,JST)) & 03941 > UEPSI) EXIT ! old/new grid differ ==> go to next grid point 03942 ZDIFTYPE_SUP= SNOW3LDIFTYP(PSNOWGRAN1(JJ,JST-1),PSNOWGRAN1(JJ, JST),& 03943 PSNOWGRAN2(JJ,JST-1),PSNOWGRAN2(JJ, JST)) 03944 ZDIFTYPE_SUP=MAX(DIFF_1,MIN(DIFF_MAX,ZDIFTYPE_SUP)) 03945 IF(PSNOWDZN(JJ,JST)>ZDZOPT(JJ,JST)*AGREG_COEF_1/ZDIFTYPE_SUP & 03946 .OR.PSNOWDZN(JJ,JST)+PSNOWDZN(JJ,JST-1)>AGREG_COEF_2*& 03947 MAX(ZDZOPT(JJ,JST),ZDZOPT(JJ,JST-1))) CYCLE! cheks upper layer 03948 ! shallow internal layer to be merged with the upper layer 03949 PSNOWDZN(JJ,JST)= PSNOWDZN(JJ,JST)+PSNOWDZN(JJ,JST-1) 03950 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 03951 ! shifts the upper layers 03952 DO JST_1=JST-1,1 03953 PSNOWDZN(JJ,JST_1)=PSNOWDZ(JJ,JST_1-1) 03954 ZDZOPT(JJ,JST_1)=ZDZOPT(JJ,JST_1-1) 03955 ENDDO 03956 EXIT ! goto to next grid point 03957 ENDDO ! end loop internal layers 03958 ENDIF 03959 ENDDO ! end grid loops for checking shallow layers 03960 03961 03962 03963 03964 ! 03965 !final check of the consistensy of the new grid size 03966 ! 03967 DO JJ=1,SIZE(PSNOW(:)) 03968 IF(ABS(SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ)))-PSNOW(JJ)) > UEPSI) THEN 03969 write(*,*) 'error in grid resizing',JJ, SUM(PSNOWDZN(JJ,1:INLVLS_USE(JJ))), & 03970 PSNOW(JJ),INLVLS_USE(JJ) 03971 write( *,*) 'PSNOWDZ:', PSNOWDZ 03972 write( *,*) 'PSNOWDZN:', PSNOWDZN 03973 CALL ABOR1_SFX("SNOWCRO: error in grid resizing") 03974 ENDIF 03975 ENDDO 03976 IF (LHOOK) CALL DR_HOOK('SNOWNLFALL_UPGRID',1,ZHOOK_HANDLE) 03977 ! 03978 END SUBROUTINE SNOWNLFALL_UPGRID 03979 03980 !############################################################################### 03981 !############################################################################### 03982 !################################################################################ 03983 !################################################################################ 03984 ! 03985 SUBROUTINE SNOWNLGRIDFRESH_1D (PSNOW,PSNOWDZ,PSNOWDZN, & 03986 PSNOWRHO,PSNOWHEAT,PSNOWGRAN1,PSNOWGRAN2, & 03987 PSNOWHIST,PSNOWAGE,GSNOWFALL, & 03988 PSNOWRHOF, PSNOWDZF,PSNOWHEATF,PSNOWGRAN1F,& 03989 PSNOWGRAN2F, PSNOWHISTF,PSNOWAGEF, & 03990 INLVLS_USE) 03991 ! 03992 !! PURPOSE 03993 !! ------- 03994 ! Snow mass,heat and characteristics redistibution in case of 03995 ! grid resizing. Total mass and heat content of the overall snowpack 03996 ! unchanged/conserved within this routine. 03997 ! Grain size and type of mixed layers is deduced from the conservation 03998 ! of the average optical size 03999 ! 04000 !! AUTHOR 04001 !! ------ 04002 !! E. Brun * Meteo-France * 04003 !! 04004 ! 04005 USE MODD_SNOW_PAR, ONLY : XSNOWCRITD,XD1,XD2,XD3,XX,XVALB5,XVALB6 04006 USE MODD_SNOW_METAMO 04007 USE MODE_SNOW3L 04008 ! 04009 ! 04010 IMPLICIT NONE 04011 ! 04012 ! 04013 !* 0.1 declarations of arguments 04014 ! 04015 REAL, INTENT(IN) :: PSNOW 04016 ! 04017 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWHEAT, PSNOWRHO, PSNOWDZ, 04018 PSNOWDZN, PSNOWGRAN1, PSNOWGRAN2, 04019 PSNOWHIST 04020 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWAGE 04021 REAL, INTENT(IN) :: PSNOWRHOF, PSNOWDZF,PSNOWHEATF, 04022 PSNOWGRAN1F,PSNOWGRAN2F, PSNOWHISTF 04023 REAL, INTENT(IN) :: PSNOWAGEF 04024 ! 04025 INTEGER, INTENT(IN) :: INLVLS_USE 04026 ! 04027 LOGICAL, INTENT(IN) :: GSNOWFALL 04028 ! 04029 !* 0.2 declarations of local variables 04030 ! 04031 INTEGER JST 04032 ! 04033 INTEGER :: INLVLS_OLD, INLVLS_NEW 04034 ! 04035 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZSNOWRHON,ZSNOWGRAN1N,ZSNOWGRAN2N, 04036 ZSNOWHEATN,ZSNOWHISTN 04037 , ZSNOWZTOP_NEW, ZSNOWZBOT_NEW 04038 REAL,DIMENSION(SIZE(PSNOWRHO,1)) ::ZSNOWAGEN 04039 REAL, DIMENSION(SIZE(PSNOWRHO,1)+1) :: ZSNOWRHOO,ZSNOWGRAN1O,ZSNOWGRAN2O, 04040 ZSNOWHEATO,ZSNOWHISTO,ZSNOWDZO 04041 , ZSNOWZTOP_OLD, ZSNOWZBOT_OLD 04042 REAL,DIMENSION(SIZE(PSNOWRHO,1)+1) ::ZSNOWAGEO 04043 ! 04044 04045 04046 INTEGER :: JST_OLD,JST_NEW 04047 REAL :: ZDENTMOYN ,ZSPHERMOYN, ZALBMOYN, ZMASTOTN,ZMASTOTO, ZSNOWHEAN, 04048 ZSNOWHEAO,ZPROPOR,ZMASDZ_OLD, ZDIAM,ZHISTMOYN 04049 REAL :: ZAGEMOYN 04050 REAL :: ZPSNOW_OLD, ZPSNOW_NEW 04051 REAL(KIND=JPRB) :: ZHOOK_HANDLE 04052 04053 !------------------------------------------------------------------------------- 04054 ! 04055 ! 0. Initialization: 04056 ! ------------------ 04057 ! 04058 ! starts by checking the consistency between both vertical grid sizes 04059 ! 04060 IF (LHOOK) CALL DR_HOOK('SNOWNLGRIDFRESH_1D',0,ZHOOK_HANDLE) 04061 INLVLS_NEW = INLVLS_USE 04062 INLVLS_OLD = -1 04063 ZPSNOW_NEW=0. 04064 ZPSNOW_OLD=0. 04065 ! 04066 DO JST_NEW=1,INLVLS_NEW 04067 ZPSNOW_NEW=ZPSNOW_NEW+PSNOWDZN(JST_NEW) 04068 ENDDO 04069 IF (ABS(ZPSNOW_NEW-PSNOWDZF)<UEPSI) THEN 04070 INLVLS_OLD=0 04071 ELSE 04072 DO JST_OLD=1,SIZE(PSNOWRHO) 04073 IF(PSNOWDZ(JST_OLD)>=UEPSI) THEN 04074 ZPSNOW_OLD=ZPSNOW_OLD+PSNOWDZ(JST_OLD) 04075 IF (ABS(ZPSNOW_NEW-PSNOWDZF-ZPSNOW_OLD)<UEPSI) THEN 04076 INLVLS_OLD=JST_OLD 04077 ENDIF 04078 ENDIF 04079 ENDDO 04080 IF (INLVLS_OLD==-1) THEN 04081 write(*,*)'pb INLVLS_OLD INLVLS_NEW=',INLVLS_NEW 04082 write(*,*)'pb INLVLS_OLD',PSNOWDZF 04083 write(*,*)'pb INLVLS_OLD',PSNOWDZN 04084 write(*,*)'pb INLVLS_OLD',PSNOWDZ 04085 CALL ABOR1_SFX('SNOWCRO: Error INLVLS_OLD') 04086 ENDIF 04087 ENDIF 04088 ! 04089 IF (GSNOWFALL) INLVLS_OLD=INLVLS_OLD+1 04090 ZPSNOW_OLD=PSNOW 04091 ZPSNOW_NEW=ZPSNOW_OLD 04092 04093 ! initialization of variables describing the initial snowpack + new snowfall 04094 04095 IF(GSNOWFALL) THEN 04096 DO JST_OLD=2, INLVLS_OLD 04097 ZSNOWDZO(JST_OLD) =PSNOWDZ(JST_OLD-1) 04098 ZSNOWRHOO(JST_OLD) =PSNOWRHO(JST_OLD-1) 04099 ZSNOWHEATO(JST_OLD) = PSNOWHEAT(JST_OLD-1) 04100 ZSNOWGRAN1O(JST_OLD) = PSNOWGRAN1(JST_OLD-1) 04101 ZSNOWGRAN2O(JST_OLD) = PSNOWGRAN2(JST_OLD-1) 04102 ZSNOWHISTO(JST_OLD) = PSNOWHIST(JST_OLD-1) 04103 ZSNOWAGEO(JST_OLD) = PSNOWAGE(JST_OLD-1) 04104 ENDDO 04105 ZSNOWDZO(1) =PSNOWDZF 04106 ZSNOWRHOO(1) =PSNOWRHOF 04107 ZSNOWHEATO(1) = PSNOWHEATF 04108 ZSNOWGRAN1O(1) = PSNOWGRAN1F 04109 ZSNOWGRAN2O(1) = PSNOWGRAN2F 04110 ZSNOWHISTO(1) = PSNOWHISTF 04111 ZSNOWAGEO(1) = PSNOWAGEF 04112 ELSE 04113 DO JST_OLD=1, INLVLS_OLD 04114 ZSNOWDZO(JST_OLD) =PSNOWDZ(JST_OLD) 04115 ZSNOWRHOO(JST_OLD) =PSNOWRHO(JST_OLD) 04116 ZSNOWHEATO(JST_OLD) = PSNOWHEAT(JST_OLD) 04117 ZSNOWGRAN1O(JST_OLD) = PSNOWGRAN1(JST_OLD) 04118 ZSNOWGRAN2O(JST_OLD) = PSNOWGRAN2(JST_OLD) 04119 ZSNOWHISTO(JST_OLD) = PSNOWHIST(JST_OLD) 04120 ZSNOWAGEO(JST_OLD) = PSNOWAGE(JST_OLD) 04121 ENDDO 04122 ENDIF 04123 04124 04125 04126 ! 04127 ! 1. Calculate vertical grid limits (m): 04128 ! -------------------------------------- 04129 ! 04130 ZSNOWZTOP_OLD(1) = ZPSNOW_OLD 04131 ZSNOWZTOP_NEW(1) = ZPSNOW_NEW 04132 ZSNOWZBOT_OLD(1) = ZSNOWZTOP_OLD(1)-ZSNOWDZO(1) 04133 ZSNOWZBOT_NEW(1) = ZSNOWZTOP_NEW(1)-PSNOWDZN(1) 04134 ! 04135 DO JST_OLD=2,INLVLS_OLD 04136 ZSNOWZTOP_OLD(JST_OLD) = ZSNOWZBOT_OLD(JST_OLD-1) 04137 ZSNOWZBOT_OLD(JST_OLD) = ZSNOWZTOP_OLD(JST_OLD)-ZSNOWDZO(JST_OLD) 04138 ENDDO 04139 ! Check consistency 04140 IF (ABS(ZSNOWZBOT_OLD(INLVLS_OLD)) > 0.00001) write (*,*) 'Error bottom OLD' 04141 ! 04142 ZSNOWZBOT_OLD(INLVLS_OLD)=0. 04143 DO JST_NEW=2,INLVLS_NEW 04144 ZSNOWZTOP_NEW(JST_NEW) = ZSNOWZBOT_NEW(JST_NEW-1) 04145 ZSNOWZBOT_NEW(JST_NEW) = ZSNOWZTOP_NEW(JST_NEW)-PSNOWDZN(JST_NEW) 04146 ENDDO 04147 ! Check consistency 04148 if (ABS(ZSNOWZBOT_NEW(INLVLS_NEW)) > 0.00001) write (*,*) 'Error bottom NEW' 04149 ZSNOWZBOT_NEW(INLVLS_NEW)=0. 04150 ! 04151 ! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing: 04152 ! -------------------------------------------------------------------- 04153 ! 04154 04155 ! loop over the new snow layers 04156 ! Summ or avergage of the constituting quantities of the old snow layers 04157 ! which are totally or partially inserted in the new snow layer 04158 04159 DO JST_NEW=1,INLVLS_NEW 04160 ZSNOWHEAN=0. 04161 ZMASTOTN=0. 04162 ZDENTMOYN=0. 04163 ZSPHERMOYN=0. 04164 ZALBMOYN=0. 04165 ZDIAM=0. 04166 ZHISTMOYN=0. 04167 ZAGEMOYN=0. 04168 04169 ! lopp over the ols snow layers 04170 DO JST_OLD=1, INLVLS_OLD 04171 IF( ZSNOWZTOP_OLD(JST_OLD) <= ZSNOWZBOT_NEW(JST_NEW)) THEN 04172 ! JST_OLD lower than JJ_NEW ==> no contribution 04173 ELSEIF ( ZSNOWZBOT_OLD(JST_OLD) >= ZSNOWZTOP_NEW(JST_NEW)) THEN 04174 ! JST_OLD higher than JJ_NEW ==> no contribution 04175 ELSE 04176 ! old layer contributes to the new one 04177 ! computation of its contributing ratio and mass/heat 04178 ZPROPOR= (MIN(ZSNOWZTOP_OLD(JST_OLD), ZSNOWZTOP_NEW(JST_NEW))& 04179 - MAX(ZSNOWZBOT_OLD(JST_OLD), ZSNOWZBOT_NEW(JST_NEW))) & 04180 / ZSNOWDZO(JST_OLD) 04181 ZMASDZ_OLD= ZPROPOR*ZSNOWRHOO(JST_OLD)*ZSNOWDZO(JST_OLD) 04182 ZMASTOTN=ZMASTOTN + ZMASDZ_OLD 04183 ZSNOWHEAN=ZSNOWHEAN+ZPROPOR*ZSNOWHEATO(JST_OLD) 04184 ! contribution to the grain optical size and then to the albedo 04185 IF(ZSNOWGRAN1O(JST_OLD)<0.) THEN 04186 ZDIAM=-ZSNOWGRAN1O(JST_OLD)*XD1/XX+(1.+ZSNOWGRAN1O(JST_OLD)/XX)* & 04187 (ZSNOWGRAN2O(JST_OLD)*XD2/XX+(1.-ZSNOWGRAN2O(JST_OLD)/XX)*XD3) 04188 ZDIAM=ZDIAM/10000. 04189 ZDENTMOYN= ZDENTMOYN-ZMASDZ_OLD*ZSNOWGRAN1O(JST_OLD)/XX 04190 ZSPHERMOYN=ZSPHERMOYN+ZMASDZ_OLD*ZSNOWGRAN2O(JST_OLD)/XX 04191 ELSE 04192 ZDIAM=ZSNOWGRAN2O(JST_OLD) 04193 ZDENTMOYN= ZDENTMOYN+ZMASDZ_OLD*0. 04194 ZSPHERMOYN=ZSPHERMOYN+ZMASDZ_OLD*ZSNOWGRAN1O(JST_OLD)/XX 04195 ENDIF 04196 ZALBMOYN=ZALBMOYN+MAX(0.,ZMASDZ_OLD*(XVALB5-XVALB6*SQRT(ZDIAM))) 04197 ZHISTMOYN=ZHISTMOYN+ZMASDZ_OLD*ZSNOWHISTO(JST_OLD) 04198 ZAGEMOYN=ZAGEMOYN+ZMASDZ_OLD*ZSNOWAGEO(JST_OLD) 04199 ENDIF 04200 ENDDO 04201 ! the new layer inherits from the weihted average properties of the old ones 04202 ! heat and mass 04203 ZSNOWHEATN(JST_NEW)= ZSNOWHEAN 04204 ZSNOWRHON(JST_NEW)= ZMASTOTN/PSNOWDZN(JST_NEW) 04205 ! grain type and size decuced from the average albedo 04206 ZALBMOYN=ZALBMOYN/ZMASTOTN 04207 ZSPHERMOYN=MAX(0.,ZSPHERMOYN/ZMASTOTN) 04208 ZDENTMOYN=MAX(0.,ZDENTMOYN/ZMASTOTN) 04209 ZDIAM=((XVALB5-ZALBMOYN)/XVALB6)**2 04210 IF (ZDIAM <XD2/10000. - 0.0000001) THEN 04211 ! dendricity is preserved if possible and sphericity is adjusted 04212 ZSNOWGRAN1N(JST_NEW)=-XX*ZDENTMOYN 04213 IF(ABS(ZSNOWGRAN1N(JST_NEW)+XX)< 0.01) THEN 04214 ZSNOWGRAN2N(JST_NEW)=XX*ZSPHERMOYN 04215 ELSEIF (ABS(ZSNOWGRAN1N(JST_NEW))< 0.0001) THEN 04216 ! dendritic snow 04217 ZSNOWGRAN1N(JST_NEW)=XX*ZSPHERMOYN 04218 ZSNOWGRAN2N(JST_NEW)=ZDIAM 04219 ELSE ! non dendritic 04220 ZSNOWGRAN2N(JST_NEW)=XX*((ZDIAM*10000.+ZSNOWGRAN1N(JST_NEW)*XD1/XX) & 04221 / (1.+ZSNOWGRAN1N(JST_NEW)/XX)-XD3)/(XD2-XD3) 04222 IF (ZSNOWGRAN2N(JST_NEW)<0.) THEN 04223 ZSNOWGRAN2N(JST_NEW)=0. 04224 ENDIF 04225 IF (ZSNOWGRAN2N(JST_NEW)> XX + 0.0000001) THEN 04226 ZSNOWGRAN2N(JST_NEW)=XX 04227 ENDIF 04228 ENDIF 04229 ELSEIF (ZDIAM >XD3/10000.) THEN 04230 ZSNOWGRAN1N(JST_NEW)=XX*ZSPHERMOYN 04231 ZSNOWGRAN2N(JST_NEW)=ZDIAM 04232 ELSEIF(ZDENTMOYN<= 0.+0.0000001) THEN 04233 ! size between D2 and D3 and dendricity == 0 04234 ZSNOWGRAN1N(JST_NEW)=XX*ZSPHERMOYN 04235 ZSNOWGRAN2N(JST_NEW)=ZDIAM 04236 ELSE 04237 ! size between D2 and D3 and dendricity < 0 04238 ! sphericity is firts preserved, if possible. If not, 04239 ! denditricity =0 04240 ZSNOWGRAN1N(JST_NEW)=-XX*ZDENTMOYN 04241 ZSNOWGRAN2N(JST_NEW)=XX*((ZDIAM*10000.+ZSNOWGRAN1N(JST_NEW)*XD1/XX) & 04242 / (1.+ZSNOWGRAN1N(JST_NEW)/XX)-XD3)/(XD2-XD3) 04243 IF ( ZSNOWGRAN2N(JST_NEW)<0..OR. ZSNOWGRAN2N(JST_NEW)> XX) THEN 04244 ! inconsistency with ZDIAM ==> dendricity = 0 04245 ZSNOWGRAN1N(JST_NEW)=XX*ZSPHERMOYN 04246 ZSNOWGRAN2N(JST_NEW)=ZDIAM 04247 ENDIF 04248 ENDIF 04249 ZSNOWHISTN(JST_NEW)=NINT(ZHISTMOYN/ZMASTOTN) 04250 ZSNOWAGEN(JST_NEW)=ZAGEMOYN/ZMASTOTN 04251 ENDDO 04252 ! 04253 ! check of consistency between new and old snowpacks 04254 ZSNOWHEAN=0. 04255 ZMASTOTN=0. 04256 ZSNOWHEAO=0. 04257 ZMASTOTO=0. 04258 ZPSNOW_NEW=0. 04259 ZPSNOW_OLD=0. 04260 04261 DO JST=1,INLVLS_NEW 04262 ZSNOWHEAN=ZSNOWHEAN+ZSNOWHEATN(JST) 04263 ZMASTOTN=ZMASTOTN+ZSNOWRHON(JST)*PSNOWDZN(JST) 04264 ZPSNOW_NEW=ZPSNOW_NEW +PSNOWDZN(JST) 04265 ENDDO 04266 DO JST=1,INLVLS_OLD 04267 ZSNOWHEAO=ZSNOWHEAO+ZSNOWHEATO(JST) 04268 ZMASTOTO=ZMASTOTO+ZSNOWRHOO(JST)*ZSNOWDZO(JST) 04269 ZPSNOW_OLD=ZPSNOW_OLD +ZSNOWDZO(JST) 04270 ENDDO 04271 if (abs(ZSNOWHEAN-ZSNOWHEAO)>0.0001.OR.ABS(ZMASTOTN-ZMASTOTO)>0.0001 & 04272 .OR. ABS(ZPSNOW_NEW-ZPSNOW_OLD)> 0.0001) THEN 04273 write(*,*) 'Warning diff', ZSNOWHEAN-ZSNOWHEAO,ZMASTOTN-ZMASTOTO,ZPSNOW_NEW-ZPSNOW_OLD 04274 ENDIF 04275 04276 04277 ! 04278 ! 5. Update mass (density and thickness) and heat: 04279 ! ------------------------------------------------ 04280 ! 04281 PSNOWDZ(:) = PSNOWDZN(:) 04282 PSNOWRHO(:) = ZSNOWRHON(:) 04283 PSNOWHEAT(:) = ZSNOWHEATN(:) 04284 PSNOWGRAN1(:) = ZSNOWGRAN1N(:) 04285 PSNOWGRAN2(:) = ZSNOWGRAN2N(:) 04286 PSNOWHIST(:) = ZSNOWHISTN(:) 04287 PSNOWAGE(:) = ZSNOWAGEN(:) 04288 IF (LHOOK) CALL DR_HOOK('SNOWNLGRIDFRESH_1D',1,ZHOOK_HANDLE) 04289 ! 04290 ! 04291 !-------------------------------------------------------------------------------- 04292 ! 04293 END SUBROUTINE SNOWNLGRIDFRESH_1D 04294 !#################################################################### 04295 !#################################################################### 04296 !################################################################### 04297 SUBROUTINE SNOWDRIFT(PTSTEP, PVMOD, PSNOWRHO,PSNOWDZ,PSNOW, & 04298 PSNOWGRAN1, PSNOWGRAN2,PSNOWHIST,INLVLS_USE & 04299 ,PTA,PQA,PPS,PRHOA) 04300 ! 04301 !! PURPOSE 04302 !! ------- 04303 ! Snow compaction and metamorphism due to drift 04304 ! Mass is unchanged: layer thickness is reduced 04305 ! in proportion to density increases. Method inspired from 04306 ! Brun et al. (1997) and Guyomarch 04307 ! 04308 ! - computes a mobility index of each snow layer from its grains, density 04309 ! and history 04310 ! - computes a drift index of each layer from its mobility and wind speed 04311 ! - computes a transport index with an exponential decay taking into 04312 ! account its depth and the mobility of upper layers 04313 ! - increases density and changes grains in case of transport 04314 ! 04315 ! HISTORY: 04316 ! Basic parameterization from Crocus/ARPEGE Coupling (1997) 04317 ! Implementation in V5 04318 ! Insertion in V6 of grains type evolution in case of dendritic snow (V. 04319 ! Vionnet) 04320 ! 07/2012 (for V7.3): E. Brun, M. Lafaysse : optional sublimation of drifted snow 04321 ! 2012-09-20 : bug correction : ZFF was not computed if LSNOWDRIFT_SUBLIM=FALSE. 04322 04323 USE MODD_CSTS,ONLY : XTT 04324 USE MODE_THERMOS 04325 04326 USE MODD_SNOW_PAR, ONLY : LSNOWDRIFT_SUBLIM 04327 ! 04328 IMPLICIT NONE 04329 ! 04330 !* 0.1 declarations of arguments 04331 ! 04332 REAL, INTENT(IN) :: PTSTEP 04333 ! 04334 REAL, DIMENSION(:), INTENT(IN) :: PVMOD 04335 ! 04336 04337 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ,PSNOWGRAN1, 04338 PSNOWGRAN2,PSNOWHIST 04339 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW 04340 ! 04341 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE 04342 ! 04343 !* 0.2 declarations of local variables 04344 ! 04345 INTEGER :: JJ,JST ! looping indexes 04346 ! 04347 ! 04348 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHO2 04349 ! 04350 REAL :: ZPROFEQU, ZRMOB, ZRDRIFT, ZRT, ZDRO, 04351 ZDGR1, ZDGR2 04352 ! modif_EB pour blowing snow sublimation 04353 REAL, DIMENSION(:), INTENT(IN) :: PTA, PQA,PPS, PRHOA 04354 REAL :: ZVT ! wind speed threshold for surface 04355 !transport 04356 REAL :: ZQS_EFFECT ! effect of QS on snow 04357 REAL :: ZWIND_EFFECT ! effect of wind on snow 04358 REAL :: ZDRIFT_EFFECT ! effect of QS and wind on snow 04359 ! transformation 04360 REAL :: ZQS !Blowing snow sublimation (kg/m2/s) 04361 REAL :: ZRHI 04362 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZQSATI,ZFF ! QS wrt ice, gust speed 04363 ! 04364 ! Calibration coefficients 04365 REAL, PARAMETER ::VTIME= 48*3600. ! characteristic time for 04366 !compaction and metamorphism by wind drift 04367 ! 04368 REAL, PARAMETER ::VROMAX= 350. ! maximum density for 04369 ! drift compaction UNIT : kg m-3 04370 REAL, PARAMETER ::VROMIN= 50. ! minimum density for 04371 ! mobility computation UNIT : kg m-3 04372 REAL, PARAMETER ::VMOB1= 0.295 ! coefficient for computing 04373 ! the mobility index 04374 REAL, PARAMETER ::VMOB2= 0.833 ! coefficient for computing 04375 ! the mobility index 04376 REAL, PARAMETER ::VMOB3= 0.583 ! coefficient for computing 04377 ! the mobility index 04378 REAL, PARAMETER ::VMOB4= -0.0583 ! coefficient for computing 04379 ! the mobility index 04380 REAL, PARAMETER ::VDRIFT1= 2.868 ! coefficient for computing 04381 ! the drift index 04382 REAL, PARAMETER ::VDRIFT2= 0.085 ! coefficient for computing 04383 ! the drift index 04384 REAL, PARAMETER ::VDRIFT3= 3.25 ! coefficient for computing 04385 ! the drift index 04386 REAL, PARAMETER ::VSIZEMIN=3.E-4 ! minimum size decrease 04387 ! by drift UNIT = m 04388 ! 04389 ! modif_EB pour sublim 04390 ! a pour but de tenir compte du fait que le vent moyen est > rafales 04391 ! on en tient compte egalement pour diminuer la duree de l'effet 04392 REAL, PARAMETER :: COEF_FF=1.25 ! coefficient for gust diagnosis from average wind 04393 REAL, PARAMETER :: COEF_EFFECT=1.0 ! coefficient for impact on density du drift 04394 REAL, PARAMETER :: ZQS_REF=2.E-5 ! valeur de reference de ZQS pour effet neige 04395 REAL(KIND=JPRB) :: ZHOOK_HANDLE 04396 04397 !------------------------------------------------------------------------------- 04398 ! 04399 ! 0. Initialization: 04400 ! ------------------ 04401 IF (LHOOK) CALL DR_HOOK('SNOWDRIFT',0,ZHOOK_HANDLE) 04402 ! 04403 ! PRINT*,LSNOWDRIFT_SUBLIM 04404 04405 DO JJ=1, SIZE(PSNOW) 04406 DO JST=1, INLVLS_USE(JJ) 04407 ZSNOWRHO2(JJ,JST) = PSNOWRHO(JJ,JST) 04408 ENDDO 04409 ENDDO 04410 ! 04411 04412 IF (LSNOWDRIFT_SUBLIM) THEN 04413 ZQSATI(:)= QSATI(PTA(:),PPS(:)) 04414 END IF 04415 ! 1. Computation of drift and induced settling and metamorphism 04416 ! ------------------ 04417 ! 04418 DO JJ=1, SIZE(PSNOW) 04419 ! 04420 ! gust speed 04421 ZFF(JJ)=PVMOD(JJ)*COEF_FF 04422 04423 ! initialization decay coeff 04424 ZPROFEQU=0. 04425 DO JST=1, INLVLS_USE(JJ) 04426 ! mobility index computation of a layer as a function of its properties 04427 IF(PSNOWGRAN1(JJ,JST) < 0.) THEN 04428 ! dendritic case 04429 ZRMOB=0.34*(0.5+0.75*(-PSNOWGRAN1(JJ,JST)/99.)-0.5* & 04430 PSNOWGRAN2(JJ,JST)/99.)+ & 04431 0.66*(1.25-1.25*(MAX(PSNOWRHO(JJ,JST),VROMIN)/1000.-VROMIN/1000.) / VMOB1) 04432 ELSE 04433 ! non dendritic case 04434 ZRMOB=0.34*(VMOB2-VMOB2*PSNOWGRAN1(JJ,JST)/99.-VMOB3* & 04435 PSNOWGRAN2(JJ,JST)*10000./100.)+ & 04436 0.66*(1.25-1.25*(MAX(PSNOWRHO(JJ,JST),VROMIN)/1000.-VROMIN/1000.) / VMOB1) 04437 ENDIF 04438 ! correction in case of former wet snow 04439 IF (PSNOWHIST(JJ,JST) >= 2) ZRMOB = MIN(ZRMOB, VMOB4) 04440 ! 04441 ! computation of drift index supposing no overburden snow 04442 ZRDRIFT = ZRMOB - (VDRIFT1 *exp(-VDRIFT2*ZFF(JJ))-1.) 04443 ! update the decay coeff by half the current layer 04444 ZPROFEQU = ZPROFEQU +0.5 * PSNOWDZ(JJ,JST) *0.1 * (VDRIFT3-ZRDRIFT) 04445 ! computation of the drift index inclunding the decay by overburden snow 04446 ZRT = max(0.,ZRDRIFT * exp(-ZPROFEQU*100)) 04447 ! modif_EB exit loop if there is no drift 04448 IF (ZRDRIFT <= 0.) EXIT 04449 04450 IF (LSNOWDRIFT_SUBLIM) THEN 04451 ! 04452 IF (JST==1) THEN 04453 ! modif_EB cas specifique surface pour blowing snow sublimation 04454 04455 ! computation of wind speed threshold QSATI and RH withe respect to ice 04456 ZVT=-LOG((ZRMOB+1.)/VDRIFT1)/VDRIFT2 04457 ZRHI=PQA(JJ)/ZQSATI(JJ) 04458 ! computation of sublimation rate according to Gordon's PhD 04459 ZQS=0.0018*(XTT/PTA(JJ))**4*ZVT*PRHOA(JJ)*ZQSATI(JJ)* & 04460 (1.-ZRHI)*(ZFF(JJ)/ZVT)**3.6 04461 ! write(*,*) 'surface Vt vent*coef ZRDRIFT ZRMOB :',ZVT,& 04462 ! ZFF(JJ),ZRDRIFT,ZRMOB 04463 ! write(*,*) 'V>Vt ZQS :',ZQS 04464 !surface depth decrease in case of blowing snow sublimation 04465 ! write(*,*) 'V>Vt DSWE DZ Z:',- MAX(0.,ZQS)*PTSTEP/COEF_FF,& 04466 ! - MAX(0.,ZQS)*PTSTEP/COEF_FF/PSNOWRHO(JJ,JST),PSNOWDZ(JJ,JST) 04467 ! 2 lignes ci-dessous a valider pour avoir sublim drift 04468 PSNOWDZ(JJ,JST)=MAX(0.5*PSNOWDZ(JJ,JST),PSNOWDZ(JJ,JST) - & 04469 MAX(0.,ZQS)*PTSTEP/COEF_FF/PSNOWRHO(JJ,JST)) 04470 04471 ENDIF 04472 ! modif_EB enhancement of snow transformation in case of sublimation 04473 ZQS_EFFECT= (MIN(3.,MAX(0.,ZQS)/ZQS_REF)) * ZRT 04474 ZWIND_EFFECT= COEF_EFFECT * ZRT 04475 ZDRIFT_EFFECT= (ZQS_EFFECT+ZWIND_EFFECT)* PTSTEP/COEF_FF/ VTIME 04476 ! write(*,*) 'ZQS_EFFECT,ZWIND_EFFECT,ZDRIFT_EFFECT:',ZQS_EFFECT,ZWIND_EFFECT,ZDRIFT_EFFECT 04477 ENDIF 04478 04479 ! settling by wind transport only in case of not too dense snow 04480 IF(PSNOWRHO(JJ,JST) < VROMAX) THEN 04481 ZDRO = ZDRIFT_EFFECT*(VROMAX - PSNOWRHO(JJ,JST)) 04482 PSNOWRHO(JJ,JST) = MIN(VROMAX , PSNOWRHO(JJ,JST)+ZDRO) 04483 PSNOWDZ(JJ,JST) = PSNOWDZ(JJ,JST)*ZSNOWRHO2(JJ,JST)/ & 04484 PSNOWRHO(JJ,JST) 04485 ENDIF 04486 ! 04487 ! metamorphism induced by snow drift 04488 IF(PSNOWGRAN1(JJ,JST) < 0.) THEN 04489 ! dendritic case 04490 ZDGR1 = ZDRIFT_EFFECT*(-PSNOWGRAN1(JJ,JST))*0.5 04491 ZDGR1 = MIN(ZDGR1,-0.99*PSNOWGRAN1(JJ,JST)) 04492 PSNOWGRAN1(JJ,JST) = PSNOWGRAN1(JJ,JST)+ZDGR1 04493 ! modif_VV_140910 04494 ZDGR2 = ZDRIFT_EFFECT*(99.-PSNOWGRAN2(JJ,JST)) 04495 PSNOWGRAN2(JJ,JST) = MIN(99.,PSNOWGRAN2(JJ,JST)+ZDGR2) 04496 ! fin modif_VV_140910 04497 ELSE 04498 ! non dendritic case 04499 ZDGR1 = ZDRIFT_EFFECT*(99.-PSNOWGRAN1(JJ,JST)) 04500 ZDGR2 = ZDRIFT_EFFECT*5./10000. 04501 PSNOWGRAN1(JJ,JST) = MIN(99.,PSNOWGRAN1(JJ,JST)+ZDGR1) 04502 PSNOWGRAN2(JJ,JST) = MAX(VSIZEMIN,PSNOWGRAN2(JJ,JST)-ZDGR2) 04503 ENDIF 04504 04505 ! update the decay coeff by half the current layer 04506 ZPROFEQU = ZPROFEQU +0.5 * PSNOWDZ(JJ,JST) *0.1 * (VDRIFT3-ZRDRIFT) 04507 ! 04508 ENDDO ! snow layers loop 04509 ENDDO ! grid points loop 04510 ! 04511 ! 2. Update total snow depth: 04512 ! ----------------------------------------------- 04513 ! 04514 ! Compaction of total snowpack depth 04515 ! 04516 DO JJ=1, SIZE(PSNOWDZ,1) 04517 PSNOW(JJ)=SUM(PSNOWDZ(JJ,1:INLVLS_USE(JJ))) 04518 ENDDO 04519 IF (LHOOK) CALL DR_HOOK('SNOWDRIFT',1,ZHOOK_HANDLE) 04520 ! 04521 END SUBROUTINE SNOWDRIFT 04522 !#################################################################### 04523 !################################################################### 04524 SUBROUTINE TRIDIAG_GROUND_SNOWCRO(PA,PB,PC,PY,PX,INLVLS_USE,IDIFLOOP) 04525 ! ######################################### 04526 ! 04527 ! 04528 !!**** *TRIDIAG_GROUND* - routine to solve a time implicit scheme 04529 !! 04530 !! 04531 !! PURPOSE 04532 !! ------- 04533 ! The purpose of this routine is to resolve the linear system: 04534 ! 04535 ! A.X = Y 04536 ! 04537 ! where A is a tridiagonal matrix, and X and Y two vertical vectors. 04538 ! However, the computations are performed at the same time for all 04539 ! the verticals where an inversion of the system is necessary. 04540 ! This explain the dimansion of the input variables. 04541 ! 04542 !!** METHOD 04543 !! ------ 04544 !! 04545 !! Then, the classical tridiagonal algorithm is used to invert the 04546 !! implicit operator. Its matrix is given by: 04547 !! 04548 !! ( b(1) c(1) 0 0 0 0 0 0 ) 04549 !! ( a(2) b(2) c(2) 0 ... 0 0 0 0 ) 04550 !! ( 0 a(3) b(3) c(3) 0 0 0 0 ) 04551 !! ....................................................................... 04552 !! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 ) 04553 !! ....................................................................... 04554 !! ( 0 0 0 0 0 ... a(n-1) b(n-1) c(n-1)) 04555 !! ( 0 0 0 0 0 ... 0 a(n) b(n) ) 04556 !! 04557 !! 04558 !! All these computations are purely vertical and vectorizations are 04559 !! easely achieved by processing all the verticals in parallel. 04560 !! 04561 !! EXTERNAL 04562 !! -------- 04563 !! 04564 !! NONE 04565 !! 04566 !! IMPLICIT ARGUMENTS 04567 !! ------------------ 04568 !! 04569 !! REFERENCE 04570 !! --------- 04571 !! 04572 !! AUTHOR 04573 !! ------ 04574 !! V. Masson 04575 !! 04576 !! MODIFICATIONS 04577 !! ------------- 04578 !! Original May 13, 1998 04579 !! 05/2011: Brun Special treatment to tackle the variable number 04580 !! of snow layers 04581 !! In case of second call, a shift of 1 snow layer 04582 !! is applied in the control loops. 04583 !! --------------------------------------------------------------------- 04584 ! 04585 !* 0. DECLARATIONS 04586 ! 04587 IMPLICIT NONE 04588 ! 04589 ! 04590 !* 0.1 declarations of arguments 04591 ! 04592 REAL, DIMENSION(:,:), INTENT(IN) :: PA ! lower diag. elements of A matrix 04593 REAL, DIMENSION(:,:), INTENT(IN) :: PB ! main diag. elements of A matrix 04594 REAL, DIMENSION(:,:), INTENT(IN) :: PC ! upper diag. elements of A matrix 04595 REAL, DIMENSION(:,:), INTENT(IN) :: PY ! r.h.s. term 04596 ! 04597 REAL, DIMENSION(:,:), INTENT(OUT) :: PX ! solution of A.X = Y 04598 ! 04599 INTEGER, DIMENSION(:), INTENT(IN) :: INLVLS_USE ! number of effective layers 04600 ! 04601 INTEGER :: IDIFLOOP ! shift in control loops: 0 or 1 04602 !* 0.2 declarations of local variables 04603 ! 04604 INTEGER :: JK ! vertical loop control 04605 INTEGER :: IN ! number of vertical levels 04606 ! 04607 REAL, DIMENSION(SIZE(PA,1) ) :: ZDET ! work array 04608 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: ZW ! work array 04609 REAL(KIND=JPRB) :: ZHOOK_HANDLE 04610 ! --------------------------------------------------------------------------- 04611 ! 04612 IF (LHOOK) CALL DR_HOOK('TRIDIAG_GROUND_SNOWCRO',0,ZHOOK_HANDLE) 04613 IN=SIZE(PX,2) 04614 ! 04615 ! 04616 !* 1. levels going up 04617 ! --------------- 04618 ! 04619 !* 1.1 first level 04620 ! ----------- 04621 ! 04622 ZDET(:) = PB(:,1) 04623 ! 04624 PX (:,1) = PY(:,1) / ZDET(:) 04625 ! 04626 !* 1.2 other levels 04627 ! ------------ 04628 ! 04629 DO JK=2,IN 04630 WHERE (JK<=INLVLS_USE(:)-IDIFLOOP) 04631 ZW(:,JK) = PC(:,JK-1)/ZDET(:) 04632 ZDET(:) = PB(:,JK ) - PA(:,JK)*ZW(:,JK) 04633 04634 PX (:,JK) = ( PY(:,JK) - PA(:,JK)*PX(:,JK-1) ) / ZDET(:) 04635 ENDWHERE 04636 END DO 04637 ! 04638 !------------------------------------------------------------------------------- 04639 ! 04640 !* 2. levels going down 04641 ! ----------------- 04642 ! 04643 DO JK=IN-1,1,-1 04644 WHERE (JK<=INLVLS_USE(:)-1-IDIFLOOP) 04645 PX (:,JK) = PX(:,JK) - ZW(:,JK+1)*PX(:,JK+1) 04646 ENDWHERE 04647 END DO 04648 IF (LHOOK) CALL DR_HOOK('TRIDIAG_GROUND_SNOWCRO',1,ZHOOK_HANDLE) 04649 ! 04650 !------------------------------------------------------------------------------- 04651 ! 04652 END SUBROUTINE TRIDIAG_GROUND_SNOWCRO 04653 04654 04655 !#################################################################### 04656 !#################################################################### 04657 !#################################################################### 04658 SUBROUTINE SNOWCROLAYER_GONE(PTSTEP,PSCAP,PSNOWTEMP,PSNOWDZ, & 04659 PSNOWRHO,PSNOWLIQ,PSNOWGRAN1,PSNOWGRAN2, & 04660 PSNOWHIST,PSNOWAGE,PLES3L, INLVLS_USE) 04661 ! 04662 ! 04663 !! PURPOSE 04664 ! Account for the case when one or several snow layers melt 04665 ! during a time step: 04666 ! in that case, merge these layers with the underlying layer 04667 ! except for the bottom layer which is merged to the abovelying layer 04668 ! energy and mass are conserved 04669 ! a new merged layer keeps the grain, histo and age properties of the 04670 ! non-melted layer 04671 ! 04672 USE MODD_CSTS,ONLY : XTT, XLMTT, XRHOLW, XRHOLI, XLVTT 04673 ! 04674 USE MODE_SNOW3L 04675 ! 04676 IMPLICIT NONE 04677 ! 04678 !* 0.1 declarations of arguments 04679 ! 04680 REAL, INTENT(IN) :: PTSTEP 04681 ! 04682 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSCAP 04683 ! 04684 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWRHO, 04685 PSNOWLIQ 04686 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE 04687 ! 04688 INTEGER, DIMENSION(:), INTENT(INOUT) :: INLVLS_USE ! 04689 ! 04690 REAL, DIMENSION(:), INTENT(IN) :: PLES3L 04691 ! 04692 !* 0.2 declarations of local variables 04693 ! 04694 REAL :: ZHEAT, ZMASS, ZDZ, ZLIQ, ZSNOWLWE 04695 ! 04696 INTEGER :: JJ,JST,JST_1, JST_2,JST_MAX,IDIFF_LAYER ! loop counter 04697 ! 04698 REAL, PARAMETER :: ZSNOWDZMIN = 0.0001 ! (m) 04699 REAL(KIND=JPRB) :: ZHOOK_HANDLE 04700 ! ZSNOWDZMIN = minimum snow layer thickness 04701 !------------------------------------------------------------------------------- 04702 ! 04703 IF (LHOOK) CALL DR_HOOK('SNOWCROLAYER_GONE',0,ZHOOK_HANDLE) 04704 DO JJ=1,SIZE(PSNOWRHO,1) ! loop on gridpoints 04705 04706 JST_MAX= INLVLS_USE(JJ) 04707 IDIFF_LAYER=0 ! used as shift counter of previously melted layers 04708 ! 04709 DO JST_1=JST_MAX,1-JST_MAX,-1 ! loop on 2 x layers in case of multi melt 04710 ! 04711 JST= JST_1+IDIFF_LAYER 04712 ! Merge is possible only in case of 2 active layers or more 04713 IF(JST>=1.AND.INLVLS_USE(JJ)>1) THEN 04714 ! Total Liquid equivalent water content of snow (m): 04715 ZSNOWLWE = PSNOWRHO(JJ,JST)*PSNOWDZ(JJ,JST)/XRHOLW 04716 ! Consideration of sublimation if any 04717 IF(JST==1) ZSNOWLWE = ZSNOWLWE- MAX(0.,PLES3L(JJ)*PTSTEP/(XLSTT*XRHOLW)) 04718 ! 04719 ! Test if avalaible energy exceeds total latent heat 04720 IF (PSCAP(JJ,JST)*MAX(0.0, PSNOWTEMP(JJ,JST) - XTT)* PSNOWDZ(JJ,JST)>= & 04721 (ZSNOWLWE-PSNOWLIQ(JJ,JST))*XLMTT*XRHOLW-UEPSI) THEN 04722 IF(JST==INLVLS_USE(JJ)) THEN 04723 ! Case of a total melt of the bottom layer: merge with above layer 04724 ! which keeps its grain, histo and age properties 04725 ZHEAT= PSNOWDZ(JJ,JST)*(PSCAP(JJ,JST)*(PSNOWTEMP(JJ,JST)-XTT) & 04726 - XLMTT*PSNOWRHO(JJ,JST) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST) & 04727 + PSNOWDZ(JJ,JST-1)*(PSCAP(JJ,JST-1)*(PSNOWTEMP(JJ,JST-1)-XTT)& 04728 - XLMTT*PSNOWRHO(JJ,JST-1) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST-1) 04729 ZMASS= PSNOWDZ(JJ,JST)*PSNOWRHO(JJ,JST) & 04730 + PSNOWDZ(JJ,JST-1)*PSNOWRHO(JJ,JST-1) 04731 ZDZ = PSNOWDZ(JJ,JST) + PSNOWDZ(JJ,JST-1) 04732 ZLIQ = PSNOWLIQ(JJ,JST)+PSNOWLIQ(JJ,JST-1) 04733 PSNOWDZ(JJ,JST-1)=ZDZ 04734 PSNOWRHO(JJ,JST-1)=ZMASS/ZDZ 04735 PSNOWLIQ(JJ,JST-1)=ZLIQ 04736 ! Temperature of the merged layer is deduced from the heat content 04737 PSCAP(JJ,JST-1)= SNOW3LSCAP(PSNOWRHO(JJ,JST-1)-PSNOWLIQ(JJ,JST-1) & 04738 * XRHOLW / MAX(PSNOWDZ(JJ,JST-1),ZSNOWDZMIN)) 04739 PSNOWTEMP(JJ,JST-1)=XTT +((((ZHEAT-XLMTT*XRHOLW*PSNOWLIQ(JJ,JST-1)) & 04740 /PSNOWDZ(JJ,JST-1)) + XLMTT*PSNOWRHO(JJ,JST-1))/PSCAP(JJ,JST-1) ) 04741 ! Decrease the number of active snow layers 04742 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 04743 ! 04744 ELSE 04745 ! Case of a total melt of the bottom layer: merge with beneath layer 04746 ! which keeps its grain, histo and age properties 04747 ! 04748 ZHEAT= PSNOWDZ(JJ,JST)*(PSCAP(JJ,JST)*(PSNOWTEMP(JJ,JST)-XTT) & 04749 - XLMTT*PSNOWRHO(JJ,JST) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST) & 04750 + PSNOWDZ(JJ,JST+1)*(PSCAP(JJ,JST+1)*(PSNOWTEMP(JJ,JST+1)-XTT)& 04751 - XLMTT*PSNOWRHO(JJ,JST+1) ) + XLMTT*XRHOLW*PSNOWLIQ(JJ,JST+1) 04752 ZMASS= PSNOWDZ(JJ,JST)*PSNOWRHO(JJ,JST) & 04753 + PSNOWDZ(JJ,JST+1)*PSNOWRHO(JJ,JST+1) 04754 ZDZ = PSNOWDZ(JJ,JST) + PSNOWDZ(JJ,JST+1) 04755 ZLIQ = PSNOWLIQ(JJ,JST)+PSNOWLIQ(JJ,JST+1) 04756 PSNOWDZ(JJ,JST)=ZDZ 04757 PSNOWRHO(JJ,JST)=ZMASS/ZDZ 04758 PSNOWLIQ(JJ,JST)=ZLIQ 04759 PSNOWGRAN1(JJ,JST)=PSNOWGRAN1(JJ,JST+1) 04760 PSNOWGRAN2(JJ,JST)=PSNOWGRAN2(JJ,JST+1) 04761 PSNOWHIST(JJ,JST)=PSNOWHIST(JJ,JST+1) 04762 PSNOWAGE(JJ,JST)=PSNOWAGE(JJ,JST+1) 04763 ! Temperature of the merged layer is deduced from the heat content 04764 PSCAP(JJ,JST)= SNOW3LSCAP(PSNOWRHO(JJ,JST)-PSNOWLIQ(JJ,JST) & 04765 * XRHOLW / MAX(PSNOWDZ(JJ,JST),ZSNOWDZMIN)) 04766 PSNOWTEMP(JJ,JST) = XTT + ( (((ZHEAT-XLMTT*XRHOLW*PSNOWLIQ(JJ,JST)) & 04767 /PSNOWDZ(JJ,JST)) + XLMTT*PSNOWRHO(JJ,JST))/PSCAP(JJ,JST) ) 04768 ! Shift the above layers 04769 DO JST_2=JST+1,INLVLS_USE(JJ)-1 04770 PSNOWTEMP(JJ,JST_2)=PSNOWTEMP(JJ,JST_2+1) 04771 PSCAP(JJ,JST_2)=PSCAP(JJ,JST_2+1) 04772 PSNOWDZ(JJ,JST_2)=PSNOWDZ(JJ,JST_2+1) 04773 PSNOWRHO(JJ,JST_2)=PSNOWRHO(JJ,JST_2+1) 04774 PSNOWLIQ(JJ,JST_2)=PSNOWLIQ(JJ,JST_2+1) 04775 PSNOWGRAN1(JJ,JST_2)=PSNOWGRAN1(JJ,JST_2+1) 04776 PSNOWGRAN2(JJ,JST_2)=PSNOWGRAN2(JJ,JST_2+1) 04777 PSNOWHIST(JJ,JST_2)=PSNOWHIST(JJ,JST_2+1) 04778 PSNOWAGE(JJ,JST_2)=PSNOWAGE(JJ,JST_2+1) 04779 ENDDO ! loop JST_2 04780 ! Decrease the number of active snow layers 04781 INLVLS_USE(JJ)=INLVLS_USE(JJ)-1 04782 ! Update the shift counter IDIFF_LAYER 04783 IDIFF_LAYER=IDIFF_LAYER+1 04784 ENDIF ! end test of bottom layer 04785 ENDIF ! end test on availibility of energy 04786 ENDIF ! end test on the number of remaining active layers 04787 ENDDO ! end loop on the snow layers 04788 ENDDO ! end loop gridpoints 04789 IF (LHOOK) CALL DR_HOOK('SNOWCROLAYER_GONE',1,ZHOOK_HANDLE) 04790 ! 04791 END SUBROUTINE SNOWCROLAYER_GONE 04792 !#################################################################### 04793 !################################################################### 04794 04795 !#################################################################### 04796 !################################################################### 04797 SUBROUTINE SNOWCROPRINTPROFILE(CINFO,NLAYERS,LPRINTGRAN,PSNOWDZ,PSNOWRHO,PSNOWTEMP,& 04798 PSNOWLIQ,PSNOWHEAT,PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE) 04799 04800 ! Matthieu Lafaysse 08/06/2012 04801 ! This routine prints the snow profile of a given point for debugging 04802 04803 IMPLICIT NONE 04804 04805 CHARACTER(*),INTENT(IN)::CINFO 04806 LOGICAL,INTENT(IN)::LPRINTGRAN 04807 INTEGER,INTENT(IN)::NLAYERS 04808 REAL,DIMENSION(:),INTENT(IN)::PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ, 04809 PSNOWHEAT,PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE 04810 04811 04812 INTEGER::JST 04813 04814 IF (LPRINTGRAN) THEN 04815 04816 WRITE(*,*) 04817 WRITE(*,*)TRIM(CINFO) 04818 WRITE(*,'(9(A12,"|"))')"------------","------------","------------",& 04819 "------------","------------","------------","------------",& 04820 "------------","------------" 04821 WRITE(*,'(9(A12,"|"))')"PSNOWDZ","PSNOWRHO","PSNOWTEMP","PSNOWLIQ","PSNOWHEAT",& 04822 "PSNOWGRAN1","PSNOWGRAN2","PSNOWHIST","PSNOWAGE" 04823 WRITE(*,'(9(A12,"|"))')"------------","------------","------------",& 04824 "------------","------------","------------","------------",& 04825 "------------","------------" 04826 DO JST=1,NLAYERS 04827 WRITE(*,'(9(ES12.3,"|")," L",I2.2)')PSNOWDZ(JST),PSNOWRHO(JST),PSNOWTEMP(JST),& 04828 PSNOWLIQ(JST),PSNOWHEAT(JST),PSNOWGRAN1(JST),PSNOWGRAN2(JST),& 04829 PSNOWHIST(JST),PSNOWAGE(JST),JST 04830 ENDDO 04831 WRITE(*,'(9(A12,"|"))')"------------","------------","------------",& 04832 "------------","------------","------------","------------",& 04833 "------------","------------" 04834 WRITE(*,*) 04835 04836 ELSE 04837 04838 WRITE(*,*) 04839 WRITE(*,*)TRIM(CINFO) 04840 WRITE(*,'(5(A12,"|"))')"------------","------------","------------",& 04841 "------------","------------" 04842 WRITE(*,'(5(A12,"|"))')"PSNOWDZ","PSNOWRHO","PSNOWTEMP","PSNOWLIQ","PSNOWHEAT" 04843 WRITE(*,'(5(A12,"|"))')"------------","------------","------------",& 04844 "------------","------------" 04845 DO JST=1,NLAYERS 04846 WRITE(*,'(5(ES12.3,"|")," L",I2.2)')PSNOWDZ(JST),PSNOWRHO(JST),PSNOWTEMP(JST),& 04847 PSNOWLIQ(JST),PSNOWHEAT(JST),JST 04848 ENDDO 04849 WRITE(*,'(5(A12,"|"))')"------------","------------","------------",& 04850 "------------","------------" 04851 WRITE(*,*) 04852 04853 END IF 04854 04855 04856 END SUBROUTINE SNOWCROPRINTPROFILE 04857 !#################################################################### 04858 !################################################################### 04859 SUBROUTINE SNOWCROPRINTATM(CINFO,PTA,PQA,PVMOD,PRR,PSR,PSW_RAD,PLW_RAD,& 04860 PTG, PSOILCOND, PD_G, PPSN3L) 04861 04862 ! Matthieu Lafaysse 08/06/2012 04863 ! This routine prints the atmospheric forcing of a given point for debugging 04864 ! and ground data 04865 04866 IMPLICIT NONE 04867 04868 CHARACTER(*),INTENT(IN)::CINFO 04869 REAL,INTENT(IN)::PTA,PQA,PVMOD,PRR,PSR,PSW_RAD,PLW_RAD 04870 REAL,INTENT(IN)::PTG, PSOILCOND, PD_G, PPSN3L 04871 04872 04873 INTEGER::JST 04874 04875 CALL SNOWCROPRINTDATE() 04876 WRITE(*,*) 04877 WRITE(*,*)TRIM(CINFO) 04878 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04879 "------------" 04880 WRITE(*,'(4(A12,"|"))')"PTA","PQA","PRR","PSR" 04881 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04882 "------------" 04883 WRITE(*,'(4(ES12.3,"|")," meteo1")')PTA,PQA,PRR,PSR 04884 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04885 "------------" 04886 WRITE(*,'(3(A12,"|"))')"------------","------------","------------" 04887 WRITE(*,'(3(A12,"|"))')"PSW_RAD","PLW_RAD","PVMOD" 04888 WRITE(*,'(3(A12,"|"))')"------------","------------","------------" 04889 WRITE(*,'(3(ES12.3,"|")," meteo2")')PSW_RAD,PLW_RAD,PVMOD 04890 WRITE(*,'(3(A12,"|"))')"------------","------------","------------" 04891 WRITE(*,*) 04892 WRITE(*,*)"Ground :" 04893 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04894 "------------" 04895 WRITE(*,'(4(A12,"|"))')"PTG","PSOILCOND","PD_G","PPSN3L" 04896 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04897 "------------" 04898 WRITE(*,'(4(ES12.3,"|")," soil")')PTG,PSOILCOND,PD_G,PPSN3L 04899 WRITE(*,'(4(A12,"|"))')"------------","------------","------------",& 04900 "------------" 04901 04902 04903 04904 04905 END SUBROUTINE SNOWCROPRINTATM 04906 04907 SUBROUTINE SNOWCROSTOPBALANCE(PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN,& 04908 PSR,PRR,PTHRUFAL,PEVAP,PEVAPCOR,PGRNDFLUX,PHSNOW,PRNSNOW,PLEL3L,PLES3L,& 04909 PHPSNOW,PSNOWHMASS,PSNOWDZ,PTSTEP) 04910 04911 ! stop if energy and mass balances are not closed 04912 04913 IMPLICIT NONE 04914 04915 REAL , DIMENSION(:), INTENT(IN) :: PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN 04916 REAL , DIMENSION(:), INTENT(IN) :: PSR,PRR,PTHRUFAL,PEVAP,PEVAPCOR 04917 REAL , DIMENSION(:), INTENT(IN) :: PGRNDFLUX,PHSNOW,PRNSNOW,PLEL3L,PLES3L,PHPSNOW,PSNOWHMASS 04918 REAL , DIMENSION(:), INTENT(IN) :: PSNOWDZ !first layer 04919 REAL , INTENT(IN) :: PTSTEP !time step 04920 04921 REAL,PARAMETER :: XWARNING_MASSBALANCE=1.E-4 04922 REAL,PARAMETER :: XWARNING_ENERGYBALANCE=1.E-4 04923 04924 04925 REAL,DIMENSION(SIZE(PSR)) :: ZMASSBALANCE,ZENERGYBALANCE 04926 04927 04928 ZMASSBALANCE= PSUMMASS_FIN-PSUMMASS_INI- & 04929 (PSR+PRR-PTHRUFAL-PEVAP)*PTSTEP - & 04930 PEVAPCOR*PTSTEP 04931 04932 ZENERGYBALANCE=(PSUMHEAT_FIN-PSUMHEAT_INI)/PTSTEP- & 04933 (-PGRNDFLUX-PHSNOW+PRNSNOW- & 04934 PLEL3L-PLES3L+PHPSNOW) & 04935 -PSNOWHMASS/PTSTEP & 04936 -PEVAPCOR*PSNOWDZ/MAX(UEPSI,PSNOWDZ)* & 04937 (ABS(PLEL3L)*XLVTT/MAX(UEPSI,ABS(PLEL3L))+ & 04938 ABS(PLES3L)*XLSTT/MAX(UEPSI,ABS(PLES3L))) 04939 04940 IF (ANY(ZMASSBALANCE>XWARNING_MASSBALANCE)) STOP "WARNING MASS BALANCE !" 04941 IF (ANY(ZENERGYBALANCE>XWARNING_ENERGYBALANCE)) STOP "WARNING ENERGY BALANCE !" 04942 04943 04944 04945 END SUBROUTINE SNOWCROSTOPBALANCE 04946 04947 !#################################################################### 04948 !################################################################### 04949 SUBROUTINE SNOWCROPRINTBALANCE(PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN,& 04950 PSR,PRR,PTHRUFAL,PEVAP,PEVAPCOR,PGRNDFLUX,PHSNOW,PRNSNOW,PLEL3L,PLES3L,& 04951 PHPSNOW,PSNOWHMASS,PSNOWDZ,PTSTEP) 04952 04953 ! Matthieu Lafaysse / Eric Brun 03/10/2012 04954 ! Print energy and mass balances. 04955 04956 IMPLICIT NONE 04957 04958 REAL , INTENT(IN) :: PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN 04959 REAL , INTENT(IN) :: PSR,PRR,PTHRUFAL,PEVAP,PEVAPCOR 04960 REAL , INTENT(IN) :: PGRNDFLUX,PHSNOW,PRNSNOW,PLEL3L,PLES3L,PHPSNOW,PSNOWHMASS 04961 REAL , INTENT(IN) :: PSNOWDZ !first layer 04962 REAL , INTENT(IN) :: PTSTEP !time step 04963 04964 REAL :: ZMASSBALANCE,ZENERGYBALANCE 04965 04966 write(*,*) ' ' 04967 write(*,FMT='(A1,67("+"),A1)')"+","+" 04968 CALL SNOWCROPRINTDATE() 04969 write(*,*) ' ' 04970 04971 ! print des residus de bilan et des differents termes pour le point 04972 write (*,FMT="(A25,1x,E17.10)") 'final mass (kg/m2) =',PSUMMASS_FIN 04973 write (*,FMT="(A25,1x,E17.10)") 'final energy (J/m2) =',ZSUMHEAT_FIN 04974 write(*,*) ' ' 04975 ZMASSBALANCE= PSUMMASS_FIN-PSUMMASS_INI- & 04976 (PSR+PRR-PTHRUFAL-PEVAP)*PTSTEP - & 04977 PEVAPCOR*PTSTEP 04978 write(*,FMT="(A25,1x,E17.10)") 'mass balance (kg/m2) =', ZMASSBALANCE 04979 04980 write(*,*) ' ' 04981 write(*,FMT="(A35)") 'mass balance contribution (kg/m2) ' 04982 write(*,FMT="(A51,1x,E17.10)") 'delta mass:',(PSUMMASS_FIN-PSUMMASS_INI) 04983 write(*,FMT="(A51,1x,E17.10)") 'hoar or condensation (>0 towards snow):',-PEVAP*PTSTEP 04984 write(*,FMT="(A51,1x,E17.10)") 'rain:',PRR*PTSTEP 04985 write(*,FMT="(A51,1x,E17.10)") 'snow:',PSR*PTSTEP 04986 write(*,FMT="(A51,1x,E17.10)") 'run-off:',PTHRUFAL*PTSTEP 04987 write(*,FMT="(A51,1x,E17.10)") 'evapcor:',PEVAPCOR*PTSTEP 04988 04989 write(*,FMT='(A1,55("-"),A1)')"+","+" 04990 write(*,*) ' ' 04991 ZENERGYBALANCE=(PSUMHEAT_FIN-PSUMHEAT_INI)/PTSTEP- & 04992 (-PGRNDFLUX-PHSNOW+PRNSNOW- & 04993 PLEL3L-PLES3L+PHPSNOW) & 04994 -PSNOWHMASS/PTSTEP & 04995 -PEVAPCOR*PSNOWDZ/MAX(UEPSI,PSNOWDZ)* & 04996 (ABS(PLEL3L)*XLVTT/MAX(UEPSI,ABS(PLEL3L))+ & 04997 ABS(PLES3L)*XLSTT/MAX(UEPSI,ABS(PLES3L))) 04998 write(*,FMT="(A25,4(1x,E17.10))") 'energy balance (W/m2)=',ZENERGYBALANCE 04999 05000 write(*,*) ' ' 05001 write(*,FMT="(A55)") 'energy balance contribution (W/m2) >0 towards snow :' 05002 write(*,FMT="(A51,1x,E17.10)") 'delta heat:',(ZSUMHEAT_FIN-ZSUMHEAT_INI)/PTSTEP 05003 write(*,FMT="(A51,1x,E17.10)") 'radiation (LW + SW):',PRNSNOW 05004 write(*,FMT="(A51,1x,E17.10)") 'sensible flux :',-PHSNOW 05005 write(*,FMT="(A51,1x,E17.10)") 'ground heat flux :',-PGRNDFLUX 05006 write(*,FMT="(A51,1x,E17.10)") 'liquid latent flux:',-PLEL3L 05007 write(*,FMT="(A51,1x,E17.10)") 'solid latent flux:',-PLES3L 05008 write(*,FMT="(A51,1x,E17.10)") 'rain sensible heat:',PHPSNOW 05009 write(*,FMT="(A51,1x,E17.10)") 'snowfall/hoar heat (sensible + melt heat):',PSNOWHMASS/PTSTEP 05010 write(*,FMT="(A51,1x,E17.10)") 'evapcor:',PEVAPCOR* & 05011 PSNOWDZ/MAX(UEPSI,PSNOWDZ)* & 05012 (ABS(PLEL3L)*XLVTT/MAX(UEPSI,ABS(PLEL3L))+ & 05013 ABS(PLES3L)*XLSTT/MAX(UEPSI,ABS(PLES3L))) 05014 write(*,FMT='(A1,67("+"),A1)')"+","+" 05015 05016 05017 05018 05019 END SUBROUTINE SNOWCROPRINTBALANCE 05020 05021 !#################################################################### 05022 !################################################################### 05023 SUBROUTINE SNOWCROPRINTDATE() 05024 05025 WRITE(*,FMT='(I4.4,2("-",I2.2)," Hour=",F5.2)') & 05026 TPTIME%TDATE%YEAR,TPTIME%TDATE%MONTH,TPTIME%TDATE%DAY,TPTIME%TIME/3600 05027 05028 END SUBROUTINE SNOWCROPRINTDATE 05029 !#################################################################### 05030 !################################################################### 05031 05032 ! 05033 END SUBROUTINE SNOWCRO
1.8.0
