! ================================================================================================================================ ! MODULE : explicitsnow ! ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr ! ! LICENCE : IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF Computes hydrologic Xsnow processes on continental points. !! !!\n DESCRIPTION: This module computes hydrologic snow processes on continental points. !! !! RECENT CHANGES: !! !! REFERENCES : None !! !! SVN : !! $HeadURL$ !! $Date$ !! $Revision$ !! \n !_ ================================================================================================================================ MODULE explicitsnow USE ioipsl_para USE constantes_soil USE constantes USE pft_parameters USE pft_parameters_var !! -SC- Ajout pour ajout throughfall dans zrainfall USE qsat_moisture USE sechiba_io_p USE xios_orchidee USE grid !! cdc Ajout pour explicitsnow_read_reftempicefile IMPLICIT NONE ! Public routines : PRIVATE PUBLIC :: explicitsnow_main, explicitsnow_initialize, explicitsnow_finalize CONTAINS !================================================================================================================================ !! SUBROUTINE : explicitsnow_initialize !! !>\BRIEF Read variables for explictsnow module from restart file !! !! DESCRIPTION : Read variables for explictsnow module from restart file !! !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_initialize( kjit, kjpindex, rest_id, snowrho, & snowtemp, snowdz, snowheat, snowgrain, & icetemp, icedz) !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier !! 0.2 Output variables REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho !! Snow density (Kg/m^3) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz !! Snow layer thickness REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat !! Snow heat content (J/m^2) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain !! Snow grainsize (m) REAL(r_std), DIMENSION (kjpindex,nice), INTENT(out) :: icetemp !! Ice temperature REAL(r_std), DIMENSION (kjpindex,nice), INTENT(out) :: icedz !! Ice layer thickness [m] !! local variables integer(i_std) :: ji !! 1. Read from restart file CALL restget_p (rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit,.TRUE.,snowrho, "gather", nbp_glo, index_g) CALL setvar_p (snowrho, val_exp, 'Snow Density profile', xrhosmin) CALL restget_p (rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit,.TRUE.,snowtemp, "gather", nbp_glo, index_g) CALL setvar_p (snowtemp, val_exp, 'Snow Temperature profile', tp_00) CALL restget_p (rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit,.TRUE.,snowdz, "gather", nbp_glo, index_g) CALL setvar_p (snowdz, val_exp, 'Snow depth profile', 0.0) CALL restget_p (rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit,.TRUE.,snowheat, "gather", nbp_glo, index_g) CALL setvar_p (snowheat, val_exp, 'Snow Heat profile', 0.0) CALL restget_p (rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit,.TRUE.,snowgrain, "gather", nbp_glo, index_g) CALL setvar_p (snowgrain, val_exp, 'Snow grain profile', 0.0) IF (ok_ice_sheet) THEN CALL restget_p (rest_id, 'icetemp', nbp_glo, nice, 1, kjit,.TRUE.,icetemp, "gather", nbp_glo, index_g) IF (ALL(icetemp(:,:)==val_exp)) THEN ! icetemp was not found in restart file IF (read_reftempice) THEN ! Read variable ptn from file CALL explicitsnow_read_reftempicefile(kjpindex,lalo,icetemp) ELSE ! Initialize icetemp with a constant value which can be set in run.def !Config Key = EXPLICITSNOW_TICEPRO !Config Desc = Initial ice temperature profile if not found in restart !Config Def = 273. !Config If = OK_SECHIBA !Config Help = The initial value of the temperature profile in the soil if !Config its value is not found in the restart file. Here !Config we only require one value as we will assume a constant !Config throughout the column. !Config Units = Kelvin [K] CALL setvar_p (icetemp, val_exp,'EXPLICITSNOW_TICEPRO',273._r_std) END IF END IF do ji=1,kjpindex icedz(ji,:) = ZSICOEF1(:) enddo ENDIF END SUBROUTINE explicitsnow_initialize !================================================================================================================================ !! SUBROUTINE : explicitsnow_main !! !>\BRIEF Main call for snow calculations !! !! DESCRIPTION : Main routine for calculation of the snow processes with explictsnow module. !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_main(kjpindex, precip_rain, precip_snow, temp_air, pb, & ! in u, v, temp_sol_new, soilcap, pgflux, & ! in frac_nobio, totfrac_nobio,frac_snow_nobio,gtemp, & ! in lambda_snow, cgrnd_snow, dgrnd_snow, & ! in lambda_ice, cgrnd_ice, dgrnd_ice, njsc, & ! in vevapsno, snow_age, snow_nobio_age,snow_nobio, snowrho, & ! inout snowgrain, snowdz, snowtemp, snowheat, snow, & ! inout temp_sol_add, icetemp, icedz, & ! inout snowliq, subsnownobio, grndflux, snowmelt, tot_melt, & ! output subsinksoil, zrainfall, frac_snow_veg, veget, veget_max) ! output !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: u,v !! Horizontal wind speed REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pgflux !! Net energy into snowpack REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ... REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ... REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegetation REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: gtemp !! First soil layer temperature REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_ice !! Coefficient of the linear extrapolation of ice surface temperature REAL(r_std), DIMENSION (kjpindex,nice), INTENT (in) :: cgrnd_ice !! Integration coefficient for ice numerical scheme REAL(r_std), DIMENSION (kjpindex,nice), INTENT (in) :: dgrnd_ice !! Integration coefficient for ice numerical scheme REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) !! 0.2 Output fields REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowliq !! Snow liquid content (m) REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...) REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: grndflux !! Net flux into soil [W/m2] REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowmelt !! Snow melt [mm/dt_sechiba] REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt from ice and snow REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil !cdc ajout zrainfall pour calcul smb REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: zrainfall !! Rain precipitation on snow !! @tex $(kg m^{-2})$ @endtex !! 0.3 Modified fields REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ... REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain !! Snow grainsize (m) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow layer thickness [m] REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat !! Snow heat content (J/m^2) REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Ice temperature REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Ice layer thickness [m] REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) !! 0.4 Local declaration INTEGER(i_std) :: ji, iv, jj,m,jv REAL(r_std),DIMENSION (kjpindex) :: snow_depth_tmp REAL(r_std),DIMENSION (kjpindex) :: snowmelt_from_maxmass REAL(r_std) :: snowdzm1 REAL(r_std), DIMENSION (kjpindex) :: thrufal !! Water leaving snow pack [kg/m2/s] REAL(r_std), DIMENSION (kjpindex) :: snowmelt_tmp,snowmelt_ice,temp_sol_new_old REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowdz_old REAL(r_std), DIMENSION (kjpindex) :: ZLIQHEATXS REAL(r_std) :: grndflux_tmp REAL(r_std), DIMENSION (nsnow) :: snowtemp_tmp REAL(r_std) :: s2flux_tmp,fromsoilflux REAL(r_std), DIMENSION (kjpindex,nsnow) :: pcapa_snow REAL(r_std), DIMENSION (kjpindex) :: psnowhmass REAL(r_std), PARAMETER :: XP00 = 1.E5 REAL(r_std), DIMENSION (kjpindex) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba] !! 1. Initialization temp_sol_new_old = temp_sol_new snowmelt_ice(:) = zero tot_melt(:) = zero snowmelt(:) = zero snowmelt_from_maxmass(:) = zero !! 2. on Vegetation ! 2.1 Snow fall CALL explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,& & snowheat,snowgrain,snowtemp,psnowhmass) ! 2.2 calculate the new snow discretization snow_depth_tmp(:) = SUM(snowdz(:,:),2) snowdz_old = snowdz ! update snowdz CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz) ! 2.3 Snow heat redistribution CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain) ! 2.4 Diagonize water portion of the snow from snow heat content: DO ji=1, kjpindex IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN snowtemp(ji,:) = snow3ltemp_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:)) snowliq(ji,:) = snow3lliq_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:)) ELSE snowliq(ji,:) = zero snowtemp(ji,:) = tp_00 ENDIF END DO ! 2.5 snow compaction !cdc snow compaction : original scheme CALL explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz) !cdc snow compaction resulting from changes in snow viscosity : compaction Decharme 2016 ! CALL explicitsnow_compactn_up(kjpindex,u,v,snowtemp,snowrho,snowdz,snowliq) !cdc snow compaction due to snowdrift (used with explicitsnow_compactn_up) ! CALL explicitsnow_drift(kjpindex,u,v,snowrho,snowdz) ! Update snow heat DO ji = 1, kjpindex snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:)) ENDDO !2.6 Calculate the snow temperature profile based on heat diffusion CALL explicitsnow_profile (kjpindex,cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add) !2.7 Test whether snow is existed on the ground or not grndflux(:)=0.0 CALL explicitsnow_gone(kjpindex,pgflux,& & snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt) !2.8 Calculate snow melt/refreezing processes CALL explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,totfrac_nobio,frac_snow_nobio, & & snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux, temp_air,frac_snow_veg,veget,veget_max, zrainfall) IF (ok_ice_sheet) THEN !cdc Calculate ice temperature CALL explicitsnow_iceprofile(kjpindex, cgrnd_ice,dgrnd_ice,lambda_ice,temp_sol_new,icedz, & & njsc, snowdz, cgrnd_snow, dgrnd_snow, snowtemp, temp_sol_add, icetemp) !cdc Calculate ice melt CALL explicitsnow_icemelt(kjpindex,njsc,snowdz,precip_snow, zrainfall, snowmelt, vevapsno, icetemp,icedz, & & grndflux,ice_sheet_melt) !cdc Update icetemp and icedz CALL explicitsnow_icelevels(kjpindex,njsc,ice_sheet_melt,icetemp,icedz) ENDIF ! 2.9 Snow sublimation changing snow thickness CALL explicitsnow_subli(kjpindex,snowrho,snowdz,vevapsno, & snowliq,snowtemp,snow,subsnownobio,subsinksoil) ! 2.9.1 Snowmelt_from_maxmass CALL explicitsnow_maxmass(kjpindex,snowrho,soilcap,snow,snowdz,snowmelt_from_maxmass) !2.10 Calculate snow grain size using the updated thermal gradient CALL explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain) !2.11 Update snow heat ! Update the heat content (variable stored each time step) ! using current snow temperature and liquid water content: ! ! First, make check to make sure heat content not too large ! (this can result due to signifigant heating of thin snowpacks): ! add any excess heat to ground flux: ! DO ji=1,kjpindex DO jj=1,nsnow ZLIQHEATXS(ji) = MAX(0.0, snowliq(ji,jj)*ph2o - 0.10*snowdz(ji,jj)*snowrho(ji,jj))*chalfu0/dt_sechiba snowliq(ji,jj) = snowliq(ji,jj) - ZLIQHEATXS(ji)*dt_sechiba/(ph2o*chalfu0) snowliq(ji,jj) = MAX(0.0, snowliq(ji,jj)) grndflux(ji) = grndflux(ji) + ZLIQHEATXS(ji) ENDDO ENDDO snow(:) = 0.0 DO ji=1,kjpindex !domain size snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:)) ENDDO DO ji = 1, kjpindex snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:)) ENDDO !! 4. On other surface types - not done yet IF ( nnobio .GT. 1 ) THEN WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW' WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES' CALL ipslerr_p(3,'explicitsnow_main','Cannot treat snow on these surface types.','','') ENDIF !! 5. Computes snow age on land and land ice (for albedo) call explicitsnow_age(kjpindex,snow,precip_snow,precip_rain,frac_snow_nobio, & temp_sol_new,snow_age,snow_nobio_age) !! 6. Check the snow on land DO ji=1,kjpindex IF (snow(ji) .EQ. 0) THEN snowrho(ji,:)=50.0 snowgrain(ji,:)=0.0 snowdz(ji,:)=0.0 snowliq(ji,:)=0.0 ENDIF ENDDO ! Snow melt only if there is more than a given mass : maxmass_snow ! Here I suggest to remove the snow based on a certain threshold of snow ! depth instead of snow mass because it is quite difficult for ! explictsnow module to calculate other snow properties following the ! removal of snow mass ! to define the threshold of snow depth based on old snow density (330 ! kg/m3) ! maxmass_snowdepth=maxmass_snow/sn_dens ! snowmelt_from_maxmass(:) = 0.0 !! 7. compute total melt DO ji=1,kjpindex !cdc tot_melt(ji) = icemelt(ji) + snowmelt(ji) + snowmelt_ice(ji) + snowmelt_from_maxmass(ji) tot_melt(ji) = snowmelt(ji) + snowmelt_from_maxmass(ji) ENDDO !cdc modif Fabienne pour TWBR snow_depth_tmp(:) = SUM(snowdz(:,:),2) snowdz_old = snowdz ! update snowdz CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz) ! snow heat redistribution CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain) IF (printlev>=3) WRITE(numout,*) 'explicitsnow_main done' ! Write output variables CALL xios_orchidee_send_field("snowliq",snowliq) CALL xios_orchidee_send_field("snowrho",snowrho) CALL xios_orchidee_send_field("snowheat",snowheat) !cdc CALL xios_orchidee_send_field("snowgrain",snowgrain) ! plantage parfois avec sortie snowgrain CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass/dt_sechiba) CALL xios_orchidee_send_field("soilcap",soilcap) END SUBROUTINE explicitsnow_main !================================================================================================================================ !! SUBROUTINE : explicitsnow_finalize !! !>\BRIEF Write variables for explictsnow module to restart file !! !! DESCRIPTION : Write variables for explictsnow module to restart file !! !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, & snowtemp, snowdz, snowheat, snowgrain, icetemp) !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content (J/m^2) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grainsize (m) REAL(r_std), DIMENSION (kjpindex,nice), INTENT(in) :: icetemp !! Ice temperature (K) !! 1. Write to restart file CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g) CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g) IF ( ok_ice_sheet) CALL restput_p(rest_id, 'icetemp', nbp_glo, nice, 1, kjit, icetemp, 'scatter', nbp_glo, index_g) END SUBROUTINE explicitsnow_finalize !================================================================================================================================ !! SUBROUTINE : explicitsnow_grain !! !>\BRIEF Compute evolution of snow grain size !! !! DESCRIPTION : !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : R. Jordan (1991) !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth (m) REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: gtemp !! First soil layer temperature REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pb !! Surface pressure (hpa) !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m) !! 0.4 Local variables REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowdz,zdz,ztheta REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: ztemp,zdiff,ztgrad,zthetaa,zfrac,& zexpo,zckt_liq,zckt_ice,zckt REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrhomin,zgrainmin INTEGER(i_std) :: ji,jj !! 0.5 Local parameters REAL(r_std), PARAMETER :: ztheta_crit = 0.02 !! m3 m-3 REAL(r_std), PARAMETER :: zc1_ice = 8.047E+9 !! kg m-3 K REAL(r_std), PARAMETER :: zc1_liq = 5.726E+8 !! kg m-3 K REAL(r_std), PARAMETER :: zdeos = 0.92E-4 !! effective diffusion !! coef for water vapor in snow !! at 0C and 1000 mb (m2 s-1) REAL(r_std), PARAMETER :: zg1 = 5.0E-7 !! m4 kg-1 REAL(r_std), PARAMETER :: zg2 = 4.0E-12 !! m2 s-1 REAL(r_std), PARAMETER :: ztheta_w = 0.05 !! m3 m-3 REAL(r_std), PARAMETER :: ztheta_crit_w = 0.14 !! m3 m-3 REAL(r_std), PARAMETER :: zdzmin = 0.01 !! m : minimum thickness !! for thermal gradient evaluation: !! to prevent excessive gradients !! for vanishingly thin snowpacks. REAL(r_std), PARAMETER :: xp00=1.E5 !! 1. initialize DO ji=1,kjpindex zsnowdz(ji,:) = MAX(xsnowdmin/nsnow, snowdz(ji,:)) DO jj=1,nsnow-1 zdz(ji,jj) = zsnowdz(ji,jj) + zsnowdz(ji,jj+1) ENDDO zdz(ji,nsnow) = zsnowdz(ji,nsnow) ! compute interface average volumetric water content (m3 m-3): ! first, layer avg VWC: ! ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:)) ! at interfaces: zthetaa(ji,0) = ztheta(ji,1) DO jj=1,nsnow-1 zthetaa(ji,jj) = (zsnowdz(ji,jj) *ztheta(ji,jj) + & zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj) ENDDO zthetaa(ji,nsnow) = ztheta(ji,nsnow) ! compute interface average temperatures (K): ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ztemp(ji,0) = snowtemp(ji,1) DO jj=1,nsnow-1 ztemp(ji,jj) = (zsnowdz(ji,jj) *snowtemp(ji,jj) + & zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj) ENDDO ztemp(ji,nsnow) = snowtemp(ji,nsnow) ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute variation of saturation vapor pressure with temperature ! for solid and liquid phases: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - zexpo(ji,:) = chalsu0/(xrv*ztemp(ji,:)) zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)**2)*(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:)) ! zexpo(ji,:) = chalev0/(xrv*ztemp(ji,:)) zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)**2)*(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:)) ! ! compute the weighted ice/liquid total variation (N m-2 K): ! zfrac(ji,:) = MIN(1.0, zthetaa(ji,:)/ztheta_crit) zckt(ji,:) = zfrac(ji,:)*zckt_liq(ji,:) + (1.0 - zfrac(ji,:))*zckt_ice(ji,:) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Compute effective diffusion coefficient (m2 s-1): ! -diffudivity relative to a reference diffusion at 1000 mb and freezing point ! multiplied by phase energy coefficient ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! DO jj=0,nsnow zdiff(ji,jj) = zdeos*(xp00/(pb(ji)*100.))*((ztemp(ji,jj)/tp_00)**6)*zckt(ji,jj) ENDDO ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Temperature gradient (K m-1): ztgrad(ji,0) = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal DO jj=1,nsnow-1 ztgrad(ji,jj) = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj)) ENDDO ! ! assume at base of snow, temperature is in equilibrium with soil ! (but obviously must be at or below freezing point!) ! ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow)) ! prognostic grain size (m) equation: !------------------------------------------------------------------- ! first compute the minimum grain size (m): ! zrhomin(ji,:) = xrhosmin zgrainmin(ji,:) = snow3lgrain_1d(zrhomin(ji,:)) ! dry snow: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! DO jj=1,nsnow IF(ztheta(ji,jj) == 0.0) THEN ! only allow growth due to vapor flux INTO layer: ! aab add sublimation(only condensation) as upper BC...? snowgrain(ji,jj) = snowgrain(ji,jj) + & (dt_sechiba*zg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj)))* & ( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) - & zdiff(ji,jj) *MIN(0.0,ztgrad(ji,jj)) ) ELSE ! wet snow ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! snowgrain(ji,jj) = snowgrain(ji,jj) + & (dt_sechiba*zg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj)))* & MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w) END IF ENDDO ENDDO END SUBROUTINE explicitsnow_grain !================================================================================================================================ !! SUBROUTINE : explicitsnow_compactn !! !>\BRIEF Compute Compaction/Settling !! !! DESCRIPTION : !! Snow compaction due to overburden and settling. !! Mass is unchanged: layer thickness is reduced !! in proportion to density increases. Method !! of Anderson (1976): see Loth and Graf, 1993, !! J. of Geophys. Res., 98, 10,451-10,464. !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): snowrho, snowdz !! !! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976) !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth !! 0.4 Local variables REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle REAL(r_std),DIMENSION(kjpindex) :: zsmassc !! cummulative snow mass (kg/m2) REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit INTEGER(i_std) :: ji,jj !! 1. initialize zsnowrho2 = snowrho zsettle(:,:) = ZSNOWCMPCT_ACM zviscocity(:,:) = ZSNOWCMPCT_V0 !! 2. Calculating Cumulative snow mass (kg/m2): DO ji=1, kjpindex IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:) zsmassc (ji)= 0.0 DO jj=1,nsnow zsmass(ji,jj) = zsmassc(ji) + zwsnowdz(ji,jj) zsmassc(ji) = zsmassc(ji) + zwsnowdz(ji,jj) ENDDO !! 3. Computing compaction/Settling ! ---------------------- ! Compaction/settling if density below upper limit ! (compaction is generally quite small above ~ 500 kg m-3): ! DO jj=1,nsnow IF (snowrho(ji,jj) .LT. xrhosmax) THEN ! ! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile) ! zsettle(ji,jj) = ZSNOWCMPCT_ACM*EXP( & -ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj))) & -ZSNOWCMPCT_CCM*MAX(0.0, & snowrho(ji,jj)-ZSNOWCMPCT_RHOD)) ! ! Snow viscocity: ! !cdc test MAR-ice60 viscosite plus faible si snowtemp proche de 0° !cdc Si T<264K on reprend l'ancienne fonction IF (snowtemp(ji,jj).GT.264.) THEN zviscocity(ji,jj) = ZSNOWCMPCT_W0*EXP( ZSNOWCMPCT_WT*(tp_00-MIN(tp_00,snowtemp(ji,jj))) + & ZSNOWCMPCT_WR*snowrho(ji,jj) ) ELSE zviscocity(ji,jj) = ZSNOWCMPCT_V0*EXP( ZSNOWCMPCT_VT*(tp_00-MIN(tp_00,snowtemp(ji,jj))) + & ZSNOWCMPCT_VR*snowrho(ji,jj) ) ENDIF ! Calculate snow density: compaction from weight/over-burden ! Anderson 1976 method: zsnowrho2(ji,jj) = snowrho(ji,jj) + snowrho(ji,jj)*dt_sechiba*( & (cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) & + zsettle(ji,jj) ) ! Conserve mass by decreasing grid thicknesses in response ! to density increases ! snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj)) ENDIF ENDDO ! Update density (kg m-3): snowrho(ji,:) = zsnowrho2(ji,:) ENDIF ENDDO END SUBROUTINE explicitsnow_compactn !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_compactn_up !! !>\BRIEF Compute Compaction/Settling !! !! DESCRIPTION : !! Snow compaction due to overburden and settling. !! Mass is unchanged: see Decharme et al. (2016) !! !! RECENT CHANGE(S) : 12/07/2021 !! !! MAIN OUTPUT VARIABLE(S): snowrho, snowdz !! !! REFERENCE(S) : Decharme et al. (2016)vi !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_compactn_up(kjpindex,u,v,snowtemp,snowrho,snowdz,snowliq) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth !! 0.4 Local variables REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsmass, zsnowrho2, zviscocity REAL(r_std),DIMENSION(kjpindex,nsnow) :: ztemp !! temperature dependance for snow viscosity REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwholdmax !! max water liquid content REAL(r_std),DIMENSION(kjpindex,nsnow) :: Fwliq !! describes the decrease of viscosity in presence of liquid water INTEGER(i_std) :: ji,jj !! 1. initialize zsnowrho2(:,:) = snowrho(:,:) zsmass(:,:) = 0.0 !! 2. Calculating Cumulative snow mass (kg/m2): DO ji=1, kjpindex IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN ! 1. Cumulative snow mass (kg/m2): DO jj=2,nsnow zsmass(ji,jj) = zsmass(ji,jj-1) + snowdz(ji,jj) * snowrho(ji,jj-1) ENDDO ! 1st layer : half the mass of the uppermost layer applied to itself zsmass(ji,1)= 0.5 * snowdz(ji,1) * snowrho(ji,1) ! 2. Compaction ! Liquid water effect zwholdmax(ji,:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) zwholdmax(ji,:) = max(1.e-10, zwholdmax(ji,:)) Fwliq(ji,:) = 1./ (1. + 10.*MIN(1.0,snowliq(ji,:)/zwholdmax(ji,:))) ! Snow viscocity, density and grid thicknesses DO jj=1,nsnow IF (snowrho(ji,jj) .LT. xrhosmax) THEN ! Temperature dependence limited to 5K: Schleef et al. (2014) ztemp(ji,jj) = ZSNOWCMPCT_XT * MIN(DTref,tp_00-MIN(tp_00,snowtemp(ji,jj))) ! Calculate snow viscocity: Brun et al. (1989), Vionnet et al. (2012) zviscocity(ji,jj) = ZSNOWCMPCT_X0 * Fwliq(ji,jj) * EXP(ztemp(ji,jj) & + ZSNOWCMPCT_XR * snowrho(ji,jj)) * snowrho(ji,jj)/xrhosref ! Calculate new snow density: zsnowrho2(ji,jj) = snowrho(ji,jj) + dt_sechiba & *(snowrho(ji,jj)*cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) ! Conserve mass by decreasing grid thicknesses in response to density increases snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj)) ENDIF ENDDO ! Update density (kg m-3): snowrho(ji,:) = zsnowrho2(ji,:) ENDIF ENDDO END SUBROUTINE explicitsnow_compactn_up !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_drift !! !>\BRIEF Compute Compaction due to snow drift : wind induced densification of near-surface snow layers !! !! DESCRIPTION : !! Snow compaction due to snowdrift !! Mass is unchanged: see Decharme et al. (2016) !! !! RECENT CHANGE(S) : 04/10/2021 !! !! MAIN OUTPUT VARIABLE(S): snowrho, snowdz !! !! REFERENCE(S) : Decharme et al. (2016) !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_drift(kjpindex,u,v,snowrho,snowdz) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth !! 0.4 Local variables REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrho2, zsettle REAL(r_std),DIMENSION(kjpindex) :: zsmobc !! cumulative mobility index REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsmob !! mobility index REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwmob !! for mobility index REAL(r_std),DIMENSION(kjpindex,nsnow) :: Zmob !! Mobility index REAL(r_std),DIMENSION(kjpindex,nsnow) :: Zwind !! Wind-driven compaction index REAL(r_std),DIMENSION(kjpindex,nsnow) :: Ztau !! Function used for the calculation of teh compaction rate REAL(r_std),DIMENSION(kjpindex,nsnow) :: tau_cmpct !! Compaction rate (s) REAL(r_std) :: speed !! Wind speed INTEGER(i_std) :: ji,jj LOGICAL :: Zwindup !! True if Zwind(jj-1) > 0 !! 1. initialize zsnowrho2(:,:) = snowrho(:,:) zsettle(:,:) = 0. zsmobc(:) = 0.0 !! 2. Calculating Cumulative snow mass (kg/m2): DO ji=1, kjpindex Zwindup = .true. IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN ! Computation of zsmob : Analogous to zsmass (used in Equation B3 Decharme 2016) DO jj=1,nsnow ! Mobility index : Equation B1 dans Appendix B Decharme 2016 Zmob(ji,jj) = Amob*(1-MAX(0.,(snowrho(ji,jj)-xrhosmin)/xrhosmob)) ! Amob= 1.25, xrhosmin= 50 kg.m-3, xrhosmob= 295 kg.m-3 zwmob(ji,jj) = snowdz(ji,jj)*(b_tau - Zmob(ji,jj)) zsmob(ji,jj) = zsmobc(ji) + zwmob(ji,jj) zsmobc(ji) = zsmobc(ji) + zwmob(ji,jj) ENDDO ! Snow viscocity, density and grid thicknesses DO jj=1,nsnow IF (snowrho(ji,jj) .LT. xrhosmax) THEN ! Calculate wind densification : Brun 1997 ! Mobility index : Equation B1 dans Appendix B Decharme 2016 Zmob(ji,jj) = Amob*(1-MAX(0.,(snowrho(ji,jj)-xrhosmin)/xrhosmob)) ! Amob= 1.25, xrhosmin= 50 kg.m-3, xrhosmob= 295 kg.m-3 ! Computation of zsmob : Analogous to zsmass (used in Equation B3 Decharme 2016) zwmob(ji,jj) = snowdz(ji,jj)*(b_tau - Zmob(ji,jj)) zsmob(ji,jj) = zsmobc(ji) + zwmob(ji,jj) zsmobc(ji) = zsmobc(ji) + zwmob(ji,jj) speed=SQRT(u(ji)*u(ji)+v(ji)*v(ji)) ! wind speed ! Wind driven compaction index : Equation B2 in Appendix B Decharme 2016 Zwind(ji,jj) = 1 - ACMPCT*EXP(-BCMPCT*KCMPCT*speed) + Zmob(ji,jj) ! ACMPCT=2.868, BCMPCT=0.085 s.m-1, KCMPCT=1.25 IF ((Zwind(ji,jj).GT.0.).AND.Zwindup) THEN ! snowdrift occurs only if Zwind>0 Ztau(ji,jj) = MAX(0., Zwind(ji,jj))*EXP(a_tau*zsmob(ji,jj)) ! Wind compaction rate (Equation B3 in Appendix B) tau_cmpct(ji,jj) = 2*KCMPCT*one_day/Ztau(ji,jj) ! zsettle : Second term in the right-hand side of Equation 13 in Decharme et al. 2016 zsettle(ji,jj) = MAX(0.,(xrhoswind-snowrho(ji,jj))/tau_cmpct(ji,jj)) !xrhoswind=350 ELSE Zwindup = .false. Ztau(ji,jj) = 0. tau_cmpct(ji,jj) = 0. zsettle(ji,jj) = 0. ENDIF ! Calculate new snow density: ! settling by wind transport only in case of not too dense snow IF (snowrho(ji,jj).LT.xrhoswind) THEN ! settling by wind cannot lead to densities above xrhoswind zsnowrho2(ji,jj) = min(xrhoswind, snowrho(ji,jj) + dt_sechiba * zsettle(ji,jj)) ! Conserve mass by decreasing grid thicknesses in response to density increases snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj)) ENDIF ENDIF ENDDO ! Update density (kg m-3): snowrho(ji,:) = zsnowrho2(ji,:) ENDIF ENDDO END SUBROUTINE explicitsnow_drift !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_transf !! !>\BRIEF Computing snow mass and heat redistribution due to grid thickness configuration resetting !! !! DESCRIPTION : Snow mass and heat redistibution due to grid thickness !! configuration resetting. Total mass and heat content !! of the overall snowpack unchanged/conserved within this routine. !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz_old !! Snow depth at the previous time step !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth (m) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2) !! 0.4 Local varibles REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzo REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzn REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowddz REAL(r_std),DIMENSION(kjpindex,nsnow) :: zdelta REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrhon,zsnowheatn,zsnowgrainn REAL(r_std),DIMENSION(kjpindex) :: zsnowmix_delta REAL(r_std),DIMENSION(kjpindex) :: zsumheat, zsumswe, zsumgrain INTEGER(i_std),DIMENSION(nsnow,2) :: locflag REAL(r_std) :: psnow INTEGER(i_std) :: ji,jj,jjj ! Initialization zsumheat(:) = 0.0 zsumswe(:) = 0.0 zsumgrain(:) = 0.0 zsnowmix_delta(:) = 0.0 locflag(:,:) = 0 DO ji=1, kjpindex psnow = SUM(snowdz(ji,:)) IF (psnow .GE. xsnowcritd .AND. snowdz_old(ji,1) .NE. 0 .AND. snowdz_old(ji,2) .NE. 0 .AND. snowdz_old(ji,3) .NE. 0) THEN zsnowzo(ji,0) = 0. zsnowzn(ji,0) = 0. zsnowzo(ji,1) = snowdz_old(ji,1) zsnowzn(ji,1) = snowdz(ji,1) DO jj=2,nsnow zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj) zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj) ENDDO DO jj=1,nsnow IF (jj .EQ. 1) THEN locflag(jj,1)=1 ELSE DO jjj=nsnow,1,-1 !upper bound of the snow layer IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN locflag(jj,1) = jjj ENDIF ENDDO ENDIF IF (jj .EQ. nsnow) THEN locflag(jj,2)=nsnow ELSE DO jjj=nsnow,1,-1 !lower bound of the snow layer IF (zsnowzn(ji,jj) .LE. zsnowzo(ji,jjj)) THEN locflag(jj,2) = jjj ENDIF ENDDO ENDIF !to interpolate ! when heavy snow occurred IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1)) zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) ELSE !snow density zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * & (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + & snowrho(ji,locflag(jj,2)) * & (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1)) DO jjj=locflag(jj,1),locflag(jj,2)-1 zsnowrhon(ji,jj) = zsnowrhon(ji,jj) + & (jjj-locflag(jj,1))*snowrho(ji,jjj)*snowdz_old(ji,jjj) ENDDO zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj) !snow heat IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * & (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + & snowheat(ji,locflag(jj,2)) * & (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2)) ELSE zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2))*(zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1))) ENDIF DO jjj=locflag(jj,1),locflag(jj,2)-1 zsnowheatn(ji,jj) = zsnowheatn(ji,jj) + & (jjj-locflag(jj,1))*snowheat(ji,jjj) ENDDO !snow grain zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * & (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + & snowgrain(ji,locflag(jj,2)) * & (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1)) DO jjj=locflag(jj,1),locflag(jj,2)-1 zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) + & (jjj-locflag(jj,1))*snowgrain(ji,jjj)*snowdz_old(ji,jjj) ENDDO zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj) ENDIF ENDDO snowrho(ji,:) = zsnowrhon(ji,:) snowheat(ji,:) = zsnowheatn(ji,:) snowgrain(ji,:) = zsnowgrainn(ji,:) ENDIF ! Vanishing or very thin snowpack check: ! ----------------------------------------- ! ! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous): ! this avoids problems related to heat and mass exchange for ! thin layers during heavy snowfall or signifigant melt: one ! new/old layer can exceed the thickness of several old/new layers. ! Therefore, mix (conservative): ! ! modified by Tao Wang IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) & & .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN zsumheat(ji) = SUM(snowheat(ji,:)) zsumswe(ji) = SUM(snowrho(ji,:)*snowdz_old(ji,:)) zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:)) zsnowmix_delta(ji) = 1.0 DO jj=1,nsnow zsnowheatn(ji,jj) = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow) snowdz(ji,jj) = zsnowmix_delta(ji)*(psnow/nsnow) zsnowrhon(ji,jj) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow) zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow) ENDDO ! Update mass (density and thickness), heat and grain size: ! ------------------------------------------------------------ snowrho(ji,:) = zsnowrhon(ji,:) snowheat(ji,:) = zsnowheatn(ji,:) snowgrain(ji,:) = zsnowgrainn(ji,:) ENDIF ENDDO END SUBROUTINE explicitsnow_transf !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_fall !! !>\BRIEF Computes snowfall !! !! DESCRIPTION : Computes snowfall !routine. !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,& & snowheat,snowgrain,snowtemp,psnowhmass) !! 0.1 Input variables INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba) REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! Horizontal wind speed REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio !! 0.2 Output variables REAL(r_std), DIMENSION(kjpindex),INTENT(out) :: psnowhmass !! Heat content of snowfall (J/m2) !! 0.3 Modified variables REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density profile (kg/m3) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow layer thickness profile (m) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2) REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m) !! 0.4 Local variables REAL(r_std), DIMENSION(kjpindex) :: rhosnew !! Snowfall density REAL(r_std), DIMENSION(kjpindex) :: dsnowfall !! Snowfall thickness (m) REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowdz_old REAL(r_std), DIMENSION(kjpindex) :: snow_depth,snow_depth_old,newgrain REAL(r_std) :: snowfall_delta,speed INTEGER(i_std) :: ji,jj !! 1. initialize the variables snowdz_old = snowdz DO ji=1,kjpindex snow_depth(ji) = SUM(snowdz(ji,:)) ENDDO snow_depth_old = snow_depth snowfall_delta = 0.0 !! 2. incorporate snowfall into snowpack DO ji = 1, kjpindex speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) ! new snow fall on snowpack ! NOTE: when the surface temperature is zero, it means that snowfall has no ! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug) psnowhmass(ji) = 0.0 IF ( (precip_snow(ji) .GT. 0.0) ) THEN ! calculate !cdc psnowhmass(ji) = precip_snow(ji)*(xci*(snowtemp(ji,1)-tp_00)-chalfu0) psnowhmass(ji) = precip_snow(ji)*(xci*(temp_air(ji)-tp_00)-chalfu0) ! Snowfall density: Following CROCUS (Pahaut 1976) rhosnew(ji) = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00) + & snowfall_c_sn*SQRT(speed)) ! Augment total pack depth: dsnowfall(ji) = precip_snow(ji)/rhosnew(ji) !snowfall thickness (m) snow_depth(ji) = snow_depth(ji) + dsnowfall(ji) ! Fresh snowfall changes the snowpack density and liquid content in uppermost layer IF (dsnowfall(ji) .GT. zero) THEN snowrho(ji,1) = (snowdz(ji,1)*snowrho(ji,1) + dsnowfall(ji)*rhosnew(ji))/ & (snowdz(ji,1)+dsnowfall(ji)) snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji) ! Add energy of snowfall to snowpack: ! Update heat content (J/m2) (therefore the snow temperature ! and liquid content): snowheat(ji,1) = snowheat(ji,1) + psnowhmass(ji) ! Incorporate snowfall grain size: newgrain(ji) = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji))) snowgrain(ji,1) = (snowdz_old(ji,1)*snowgrain(ji,1) + dsnowfall(ji)*newgrain(ji))/ & snowdz(ji,1) ENDIF ELSE dsnowfall(ji) = 0. ENDIF ! new snow fall on snow free surface. ! we use the linearization for the new snow fall on snow-free ground IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN snowfall_delta = 1.0 DO jj=1,nsnow snowdz(ji,jj) = snowfall_delta*(dsnowfall(ji)/nsnow) + & (1.0-snowfall_delta)*snowdz(ji,jj) snowheat(ji,jj) = snowfall_delta*(psnowhmass(ji)/nsnow) + & (1.0-snowfall_delta)*snowheat(ji,jj) snowrho(ji,jj) = snowfall_delta*rhosnew(ji) + & (1.0-snowfall_delta)*snowrho(ji,jj) snowgrain(ji,jj) = snowfall_delta*newgrain(ji) + & (1.0-snowfall_delta)*snowgrain(ji,jj) ENDDO ENDIF ENDDO END SUBROUTINE explicitsnow_fall !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_gone !! !>\BRIEF Check whether snow is gone !! !! DESCRIPTION : If so, set thickness (and therefore mass and heat) and liquid !! content to zero, and adjust fluxes of water, evaporation and !! heat into underlying surface. !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_gone(kjpindex,pgflux,& snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt) !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pgflux !! Net energy into snow pack(w/m2) REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowheat !! Snow heat content (J/m^2) !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content REAL(r_std),DIMENSION(kjpindex), INTENT(inout) :: grndflux !! Soil/snow interface heat flux (W/m2) REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: snowmelt !! Snow melt REAL(r_std),DIMENSION(kjpindex) :: thrufal !! Water leaving snowpack(kg/m2/s) !! 0.4 Local variables INTEGER(i_std) :: ji,jj REAL(r_std),DIMENSION(kjpindex) :: snowgone_delta REAL(r_std),DIMENSION (kjpindex) :: totsnowheat !!snow heat content at each layer REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit ! first caculate total snowpack snow heat content snowgone_delta(:) = un thrufal(:)=0.0 snowmelt(:)=0 totsnowheat(:) = SUM(snowheat(:,:),2) DO ji = 1, kjpindex IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN IF ( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba) thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:)) snowgone_delta(ji) = 0.0 snowmelt(ji) = snowmelt(ji)+thrufal(ji) ENDIF ! update of snow state (either still present or not) DO jj=1,nsnow snowdz(ji,jj) = snowdz(ji,jj) *snowgone_delta(ji) snowliq(ji,jj) = snowliq(ji,jj) *snowgone_delta(ji) snowtemp(ji,jj) = (1.0-snowgone_delta(ji))*tp_00 + snowtemp(ji,jj)*snowgone_delta(ji) ENDDO ELSE snowdz(ji,:)=0. snowliq(ji,:)=0. snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:)) !cdc MODIF ICE grndflux(ji) = pgflux(ji) ENDIF ENDDO END SUBROUTINE explicitsnow_gone !================================================================================================================================ !! SUBROUTINE : explicitsnow_melt_refrz !! !>\BRIEF Computes snow melt and refreezing processes within snowpack !! !! DESCRIPTION : !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): None !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,totfrac_nobio,frac_snow_nobio, & snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air, frac_snow_veg, veget, veget_max, zrainfall) !! 0.1 Input variables INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size REAL(r_std),DIMENSION (kjpindex,nsnow) :: pcapa_snow !! Heat capacity for snow REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pgflux !! Net energy into snowpack(w/m2) REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ... REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegetation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type !! 0.2 Output variables REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: zrainfall !! Rain precipitation on snow !! @tex $(kg m^{-2})$ @endtex !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow layer density (Kg/m^3) REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: snowmelt !! Snowmelt !! 0.4 Local variables REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil REAL(r_std) :: enerin,melttot,hrain REAL(r_std),DIMENSION (nsnow) :: zsnowlwe REAL(r_std),DIMENSION (nsnow) :: flowliq REAL(r_std),DIMENSION (kjpindex) :: snowmass REAL(r_std),DIMENSION (nsnow) :: zphase !! Phase change (from ice to water) (J/m2) REAL(r_std),DIMENSION (nsnow) :: zphase2 !! Phase change (from water to ice) REAL(r_std),DIMENSION (nsnow) :: zphase3 !! Phase change related with net energy input to snowpack REAL(r_std),DIMENSION (nsnow) :: zsnowdz !! Snow layer depth REAL(r_std),DIMENSION (nsnow) :: zsnowmelt !! Snow melt (liquid water) (m) REAL(r_std),DIMENSION (nsnow) :: zsnowtemp REAL(r_std),DIMENSION (nsnow) :: zmeltxs !! Excess melt REAL(r_std),DIMENSION (nsnow) :: zwholdmax !! Maximum liquid water holding (m) REAL(r_std),DIMENSION (nsnow) :: zcmprsfact !! Compression factor due to densification from melting REAL(r_std),DIMENSION (nsnow) :: zscap !! Snow heat capacity (J/m3 K) REAL(r_std),DIMENSION (nsnow) :: zsnowliq !! (m) REAL(r_std),DIMENSION (nsnow) :: snowtemp_old REAL(r_std),DIMENSION (0:nsnow) :: zflowliqt !!(m) REAL(r_std),DIMENSION (kjpindex) :: frac_rain_veg REAL(r_std) :: zpcpxs REAL(r_std) :: ztotwcap REAL(r_std),DIMENSION(kjpindex,nsnow) :: snowdz_old,snowliq_old INTEGER(i_std) :: jj,ji, iv, jv REAL(r_std),DIMENSION(nsnow) :: snowdz_old2 REAL(r_std),DIMENSION(nsnow) :: zsnowrho REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrefrz !! (m) !initialize snowdz_old = snowdz snowliq_old = snowliq zrainfall(:) = 0. frac_rain_veg(:) = 0. zrefrz(:,:) = 0. DO ji = 1, kjpindex snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:)) IF ((snowmass(ji) .GT. min_sechiba)) THEN !! 1 snow melting due to positive snowpack snow temperature !! 1.0 total liquid equivalent water content of each snow layer zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o !! 1.1 phase change (J/m2) pcapa_snow(ji,:) = snowrho(ji,:)*xci zphase(:) = MIN(pcapa_snow(ji,:)*MAX(0.0, snowtemp(ji,:)-tp_00)* & snowdz(ji,:), & MAX(0.0,zsnowlwe(:)-snowliq(ji,:))*chalfu0*ph2o) !! 1.2 update snow liq water and temperature if melting zsnowmelt(:) = zphase(:)/(chalfu0*ph2o) !! 1.3 cool off the snow layer temperature due to melt zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:)) snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:)) zmeltxs(:) = (zsnowtemp(:)-snowtemp(ji,:))*pcapa_snow(ji,:)*snowdz(ji,:) !! 1.4 loss of snowpack depth and liquid equivalent water zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension !!-SC- Modification de zwholdmax selon Vionnet_et_al_GMD_2012 !!-SC- zwholdmax(:)= 0.05*snowdz(ji,:)*(1-(snowrho(ji,:)/920)) zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ & (zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:))) snowdz(ji,:) = snowdz(ji,:)*zcmprsfact(:) snowrho(ji,:) = zsnowlwe(:)*ph2o/snowdz(ji,:) snowliq(ji,:) = snowliq(ji,:) + zsnowmelt(:) !! 2 snow refreezing process !! 2.1 freeze liquid water in any layer zscap(:) = snowrho(ji,:)*xci !J K-1 m-3 zphase2(:) = MIN(zscap(:)* & MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:), & snowliq(ji,:)*chalfu0*ph2o) ! warm layer and reduce liquid if freezing occurs zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:)) snowtemp_old(:) = snowtemp(ji,:) snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:)) ! Reduce liquid portion if freezing occurs: snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))* & zscap(:)*zsnowdz(:)/(chalfu0*ph2o) ) !!-SC- Refreezing : zrefrz(ji,:) = (snowtemp(ji,:)-snowtemp_old(:))*zscap(:)*zsnowdz(:)/(chalfu0*ph2o) !!-SC- snowliq(ji,:) = MAX(snowliq(ji,:), 0.0) !! 3. thickness change due to snowmelt in excess of holding capacity zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:))) snowliq(ji,:) = snowliq(ji,:) - flowliq(:) snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:) ! to prevent possible very small negative values (machine ! prescision as snow vanishes snowdz(ji,:) = MAX(0.0, snowdz(ji,:)) !! 4. liquid water flow ztotwcap = SUM(zwholdmax(:)) ! Rain entering snow (m): !cdc version svn zrainfall = precip_rain(ji)/ph2o*frac_snow_nobio(ji,iice)*totfrac_nobio(ji) ! rainfall (m) !!-SC-- zrainfall-06 (Le bon) !~ DO jv = 1,nvm !~ zrainfall = (frac_snow_nobio(ji,iice)*totfrac_nobio(ji) & !~ + frac_snow_veg(ji)*throughfall_by_pft(jv)*veget(ji,jv) & !~ + frac_snow_veg(ji)*(veget_max(ji,jv) - veget(ji,jv)) ) & !~ * precip_rain(ji)/ph2o !~ ENDDO !!-SC-- fin zrainfall-06 !cdc version Christophe zrainfall tableau et somme sur jv DO jv = 1,nvm frac_rain_veg(ji) = frac_snow_veg(ji)*throughfall_by_pft(jv)*veget(ji,jv) & + frac_snow_veg(ji)*(veget_max(ji,jv) - veget(ji,jv)) & + frac_rain_veg(ji) ENDDO zrainfall(ji) = (frac_snow_nobio(ji,iice)*totfrac_nobio(ji) & + frac_rain_veg(ji) )* precip_rain(ji)/ph2o zflowliqt(0) = MIN(zrainfall(ji),ztotwcap) ! Rain assumed to directly pass through the pack to runoff (m): zpcpxs = zrainfall(ji) - zflowliqt(0) ! DO jj=1,nsnow zflowliqt(jj) = flowliq(jj) ENDDO ! translated into a density increase: flowliq(:) = 0.0 ! clear this array for work zsnowliq(:) = snowliq(ji,:) ! reset liquid water content ! DO jj=1,nsnow snowliq(ji,jj) = snowliq(ji,jj) + zflowliqt(jj-1) flowliq(jj) = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj))) snowliq(ji,jj) = snowliq(ji,jj) - flowliq(jj) !Modified by Tao Wang based on previous ISBA-ES scheme snowrho(ji,jj) = snowrho(ji,jj) + (snowliq(ji,jj) - zsnowliq(jj))* & & ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj)) zflowliqt(jj) = zflowliqt(jj) + flowliq(jj) ENDDO snowmelt(ji) = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) * ph2o ! excess heat from melting, using it to warm underlying ground to conserve energy meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2) ! energy flux into the soil grndflux(ji) = grndflux(ji) + meltxs(ji) ELSE snowdz(ji,:)=0. snowliq(ji,:)=0. snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:)) !This addition is to get the precipitation that falls on the snow in case the snow has just totally disppeared in explicitsnow_gone. snowmelt(ji)=snowmelt(ji)+ precip_rain(ji)*frac_snow_nobio(ji,iice)*totfrac_nobio(ji) zrefrz(ji,:) = 0. ENDIF ENDDO !cdc sortie du zrainfall : CALL xios_orchidee_send_field("zrainfall",zrainfall) !!-SC-Refreezing: CALL xios_orchidee_send_field("zrefrz",zrefrz) END SUBROUTINE explicitsnow_melt_refrz !================================================================================================================================ !! SUBROUTINE : explicitsnow_icemelt !! !>\BRIEF Computes ice melt on ice sheet area (no refreezing process) !! !! DESCRIPTION : !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): ice_sheet_melt, icetemp, icedz, grndflux !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_icemelt(kjpindex,njsc,snowdz, precip_snow, zrainfall, snowmelt, vevapsno, icetemp,icedz,& grndflux, ice_sheet_melt) !! 0.1 Input variables INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba) REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: zrainfall !! Rain precipitation on snow !! @tex $(kg m^{-2})$ @endtex REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snowmelt !! Snowmelt REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex !! 0.2 Output variables REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba] !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Snow depth REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil !! 0.4 Local variables REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil REAL(r_std),DIMENSION (nice) :: zphase_ice !! Phase change (from ice to water) (J/m2) REAL(r_std),DIMENSION (nice) :: zicemelt !! Ice melt (liquid water) (m) REAL(r_std),DIMENSION (nice) :: zicetemp REAL(r_std),DIMENSION (nice) :: zmeltxs !! Excess melt INTEGER(i_std) :: ji REAL(r_std),DIMENSION (kjpindex,nice) :: pcapa_ice !! Heat capacity for ice REAL(r_std), DIMENSION (kjpindex) :: surf_massbal !! surface mass balance [mm/dt_sechiba] !initialize ice_sheet_melt(:) = 0. DO ji = 1, kjpindex IF (njsc(ji).EQ.13.) THEN ! ice grid point !! 1 Ice melting due to positive ice temperature !! 1.1 phase change (J/m2) ! pcapa_ice(ji,:) = 917. * xci pcapa_ice(ji,:) = rho_ice * (ZICETHRMHEAT1 + ZICETHRMHEAT2*(icetemp(ji,:)-tp_00)) ! formulation as in GRISLI prop_th_icetemp.f90 zphase_ice(:) = pcapa_ice(ji,:)*MAX(0.0, icetemp(ji,:)-tp_00) * icedz(ji,:) !! 1.2 ice melt J/m2/(J/kg)=kg/m2 = mm zicemelt(:) = zphase_ice(:)/chalfu0 !! 1.3 cool off the snow layer temperature due to melt zicetemp(:) = icetemp(ji,:) - zphase_ice(:)/(pcapa_ice(ji,:)* icedz(ji,:)) icetemp(ji,:) = MIN(tp_00, zicetemp(:)) zmeltxs(:) = (zicetemp(:)-icetemp(ji,:))*pcapa_ice(ji,:)*icedz(ji,:) icedz(ji,:) = max(0.,icedz(ji,:) - zicemelt(:)/ rho_ice) ice_sheet_melt(ji) = sum(zicemelt(:)) ! excess heat from melting, using it to warm underlying ground to conserve energy meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2) ! energy flux into the soil grndflux(ji) = grndflux(ji) + meltxs(ji) ELSE ice_sheet_melt(ji) = 0. ENDIF ENDDO surf_massbal(:) = precip_snow(:) + zrainfall(:) - snowmelt(:) - ice_sheet_melt(:) - vevapsno(:) CALL xios_orchidee_send_field("surf_massbal",surf_massbal/dt_sechiba) CALL xios_orchidee_send_field("ice_sheet_melt",ice_sheet_melt/dt_sechiba) END SUBROUTINE explicitsnow_icemelt !================================================================================================================================ !! SUBROUTINE : explicitsnow_icelevels !! !>\BRIEF Update icelevels (constant) and icetemp !! !! DESCRIPTION : !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): icetemp, icedz !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_icelevels(kjpindex,njsc,ice_sheet_melt,icetemp,icedz) !! 0.1 Input variables INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba] !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Snow depth !! 0.4 Local variables INTEGER(i_std) :: ji,jj,xx,locxx REAL(r_std),DIMENSION(nice) :: znt_tmp !! Ice grid layer boundary after melt REAL(r_std),DIMENSION (kjpindex,nice) :: icetemp_new !! Update ice temperature REAL(r_std) :: totdz !! melt cumulated on all ice layers DO ji = 1, kjpindex ! on ne recalcule icetemp que pour le type Ice (classification usda n°13) IF (njsc(ji).EQ.13) THEN IF (ice_sheet_melt(ji).GT.min_sechiba) THEN totdz = 0. DO jj=1,nice-1 icetemp_new(ji,jj) = icetemp(ji,jj)*icedz(ji,jj)/ZSICOEF1(jj) + icetemp(ji,jj+1)*(ZSICOEF1(jj)-icedz(ji,jj))/ZSICOEF1(jj) totdz = totdz + ZSICOEF1(jj)-icedz(ji,jj) ! melt cumulated on all ice layers ENDDO icetemp_new(ji,nice) = icetemp(ji,nice)*(ZSICOEF1(nice) - totdz)/ZSICOEF1(nice) + tp_00 * totdz/ZSICOEF1(nice) icetemp(ji,:)=icetemp_new(ji,:) icedz(ji,:)=ZSICOEF1(:) !cdc temperature base glace = 0°C : ! icetemp(ji,nice)= tp_00 !cdc test en forcant la glace à 0°C sur toute l'epaisseur !cdc icetemp(ji,:)=tp_00 ENDIF ENDIF ENDDO END SUBROUTINE explicitsnow_icelevels !================================================================================================================================ !! SUBROUTINE : explicitsnow_levels !! !>\BRIEF Computes snow discretization based on given total snow depth !! !! DESCRIPTION : !! RECENT CHANGE(S) : compatible with 3 and 12 layers !! !! MAIN OUTPUT VARIABLE(S): snowdz !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz) !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_thick !! Total snow depth !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth !! 0.4 Local variables INTEGER(i_std) :: ji, jj !! parameters for 3 layers snow model REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF1 = (/0.25, 0.50, 0.25/) !! Snow grid parameters REAL(r_std), PARAMETER, DIMENSION(2) :: ZSGCOEF2 = (/0.05, 0.34/) !! Snow grid parameters REAL(r_std), PARAMETER :: ZSNOWTRANS = 0.20 !! Minimum total snow depth at which surface layer thickness is constant (m) REAL(r_std), PARAMETER :: XSNOWCRITD = 0.03 !! (m) !! parameters for 12 layers snon model REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF11 = (/0.3, 0.4, 0.3/) !! Snow grid parameters : distribution of snow on layers 6, 7 and 8 REAL(r_std), PARAMETER, DIMENSION(12) :: Zmax = (/0.01, 0.05, 0.15, 0.5, 1., 0., 0., 0., 1., 0.5, 0.1, 0.02/) !! Snow grid parameters : max thickness of each layer REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF22 = (/0.3, 0.3, 0.3/) LOGICAL, PARAMETER, DIMENSION(12) :: mask_d_r = (/.true.,.true.,.true.,.true.,.true.,.false.,.false.,.false.,.true.,.true.,.true.,.true./) ! True for nsnow=1,5 and nsnow=9,12 REAL(r_std), DIMENSION(kjpindex) :: d_r ! sum of the snow depth of layers 6, 7 and 8. IF (nsnow .eq. 3) THEN DO ji=1,kjpindex IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow) snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow) snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3) ENDIF ENDDO WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. & snow_thick(:) .GT. (XSNOWCRITD+0.01) ) snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1) snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2) snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3) END WHERE DO ji = 1,kjpindex IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN snowdz(ji,1) = ZSGCOEF2(1) snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1) ! When using simple finite differences, limit the thickness ! factor between the top and 2nd layers to at most 10 snowdz(ji,2) = MIN(10*ZSGCOEF2(1), snowdz(ji,2) ) snowdz(ji,3) = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1) ENDIF ENDDO ELSE IF (nsnow .eq. 12) THEN ! snow layering as in Decharme 2016 3.1.1 DO ji=1,kjpindex !! test if snowdz must be updated IF ((snowdz(ji,1).LT.(0.5*min(Zmax(1),snow_thick(ji)/12)).OR.(snowdz(ji,1).GT.(1.5*min(Zmax(1),snow_thick(ji)/12)))) & .OR. (snowdz(ji,2).LT.(0.5*min(Zmax(2),snow_thick(ji)/12)).OR.(snowdz(ji,2).GT.(1.5*min(Zmax(2),snow_thick(ji)/12)))) & .OR. (snowdz(ji,12).LT.(0.5*min(Zmax(12),snow_thick(ji)/12)).OR.(snowdz(ji,12).GT.(1.5*min(Zmax(12),snow_thick(ji)/12))))) THEN DO jj=1,5 snowdz(ji,jj) = min(Zmax(jj),snow_thick(ji)/12) ENDDO DO jj=9,nsnow snowdz(ji,jj) = min(Zmax(jj),snow_thick(ji)/12) ENDDO ! d_r(ji) = snow_thick(ji) - sum(snowdz(ji,:),mask=mask_d_r) ! version code std 12 couches !cdc version compatible avec 3 couches a la compilation d_r(ji) = snow_thick(ji) do jj=1,nsnow if (mask_d_r(jj)) then d_r(ji) = d_r(ji) - snowdz(ji,jj) endif enddo snowdz(ji,6) = ZSGCOEF11(1)*d_r(ji) - min(0., ZSGCOEF22(1)*d_r(ji) - snowdz(ji,5)) snowdz(ji,7) = ZSGCOEF11(2)*d_r(ji) + min(0., ZSGCOEF22(2)*d_r(ji) - snowdz(ji,5)) + min(0.,ZSGCOEF22(2)*d_r(ji) - snowdz(ji,9)) snowdz(ji,8) = ZSGCOEF11(3)*d_r(ji) - min(0., ZSGCOEF22(3)*d_r(ji) - snowdz(ji,9)) ENDIF ENDDO ENDIF END SUBROUTINE explicitsnow_levels !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_profile !! !>\BRIEF !! !! DESCRIPTION : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile. !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): snowtemp, temp_sol_add !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add) !! 0. Variables and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! skin temperature REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout) :: snowtemp REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) !! 0.4 Local variables INTEGER(i_std) :: ji, jg !_ !================================================================================================================================ !! 1. Computes the snow temperatures ptn. DO ji = 1,kjpindex IF (SUM(snowdz(ji,:)) .GT. 0) THEN snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / & (lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un) temp_sol_add(ji) = zero DO jg = 1,nsnow-1 snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg) ENDDO ENDIF ENDDO IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done ' END SUBROUTINE explicitsnow_profile !! !================================================================================================================================ !! SUBROUTINE : explicitsnow_iceprofile !! !>\BRIEF !! !! DESCRIPTION : In this routine solves the numerical ice thermal scheme, ie calculates the new ice temperature profile. !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): icetemp, temp_sol_add !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ !================================================================================================================================ SUBROUTINE explicitsnow_iceprofile (kjpindex, cgrnd_ice,dgrnd_ice,lambda_ice,temp_sol_new,icedz,& njsc, snowdz, cgrnd_snow, dgrnd_snow, snowtemp, temp_sol_add, icetemp) !! 0. Variables and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! skin temperature REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: cgrnd_ice !! Integration coefficient for ice numerical scheme REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: dgrnd_ice !! Integration coefficient for ice numerical scheme REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_ice !! Coefficient of the linear extrapolation of surface temperature REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: icedz !! ice layer thickness INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature !! 0.2 Output variables !! 0.3 Modified variables REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt ice for ice ablation case (K) REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Ice temperature !! 0.4 Local variables INTEGER(i_std) :: ji, jg !_ !================================================================================================================================ !! 1. Computes the snow temperatures ptn. DO ji = 1,kjpindex IF (njsc(ji).EQ.13.) THEN ! ice grid point IF (sum(snowdz(ji,:)).GT.0.) THEN ! grid point with snow icetemp(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * snowtemp(ji,nsnow) temp_sol_add(ji) = zero DO jg = 1,nice-1 icetemp(ji,jg+1) = cgrnd_ice(ji,jg) + dgrnd_ice(ji,jg) * icetemp(ji,jg) ENDDO ELSE ! no snow icetemp(ji,1) = (lambda_ice(ji) * cgrnd_ice(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / & (lambda_ice(ji) * (un - dgrnd_ice(ji,1)) + un) temp_sol_add(ji) = zero DO jg = 1,nice-1 icetemp(ji,jg+1) = cgrnd_ice(ji,jg) + dgrnd_ice(ji,jg) * icetemp(ji,jg) ENDDO ENDIF ELSE icetemp(ji,:) = tp_00 ! ice temperature on non ice-sheet area (could be an other value, to be tested) ENDIF ENDDO IF (printlev>=3) WRITE (numout,*) ' explicitsnow_iceprofile done ' END SUBROUTINE explicitsnow_iceprofile !! ================================================================================================================================ !! SUBROUTINE : explicitsnow_read_reftempicefile !! !>\BRIEF !! !! DESCRIPTION : Read file with longterm ice temperature !! !! !! RECENT CHANGE(S) : None !! !! MAIN OUTPUT VARIABLE(S): reftempice : Reference temperature for ice !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE explicitsnow_read_reftempicefile(kjpindex,lalo,reftempice) USE interpweight IMPLICIT NONE !! 0. Variables and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo !! 0.2 Output variables REAL(r_std), DIMENSION(kjpindex, nice), INTENT(out) :: reftempice !! 0.3 Local variables INTEGER(i_std) :: ib CHARACTER(LEN=80) :: filename REAL(r_std),DIMENSION(kjpindex) :: reftempice_file !! Horizontal temperature field interpolated from file [K] INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours REAL(r_std) :: vmin, vmax !! min/max values to use for the !! renormalization REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation CHARACTER(LEN=80) :: variablename !! Variable to interpolate !! the file CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable !! (variabletypevals(1) = -un, not used) CHARACTER(LEN=50) :: fractype !! method of calculation of fraction !! 'XYKindTime': Input values are kinds !! of something with a temporal !! evolution on the dx*dy matrix' LOGICAL :: nonegative !! whether negative values should be removed CHARACTER(LEN=50) :: maskingtype !! Type of masking !! 'nomask': no-mask is applied !! 'mbelow': take values below maskvals(1) !! 'mabove': take values above maskvals(1) !! 'msumrange': take values within 2 ranges; !! maskvals(2) <= SUM(vals(k)) <= maskvals(1) !! maskvals(1) < SUM(vals(k)) <= maskvals(3) !! (normalized by maskvals(3)) !! 'var': mask values are taken from a !! variable inside the file (>0) REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to !! `maskingtype') CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask REAL(r_std) :: reftemp_norefinf REAL(r_std) :: reftemp_default !! Default value !Config Key = SOIL_REFTEMP_FILE !Config Desc = File with climatological soil temperature !Config If = READ_REFTEMPICE !Config Def = reftempice.nc !Config Help = !Config Units = [FILE] filename = 'reftempice.nc' CALL getin_p('REFTEMPICE_FILE',filename) variablename = 'icetemp' IF (printlev >= 3) WRITE(numout,*) " in thermosoil_read_reftempicefile filename '" // TRIM(filename) // & "' variable name: '" //TRIM(variablename) // "'" ! For this case there are not types/categories. We have 'only' a continuos field ! Assigning values to vmin, vmax vmin = 0. vmax = 9999. ! For this file we do not need neightbours! neighbours = 0 !! Variables for interpweight ! Type of calculation of cell fractions fractype = 'default' ! Name of the longitude and latitude in the input file lonname = 'nav_lon' latname = 'nav_lat' ! Default value when no value is get from input file reftemp_default = 1. ! Reference value when no value is get from input file reftemp_norefinf = 1. ! Should negative values be set to zero from input file? nonegative = .FALSE. ! Type of mask to apply to the input data (see header for more details) maskingtype = 'nomask' ! Values to use for the masking (here not used) maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /) ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used) namemaskvar = '' CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, & contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, & maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, & reftempice_file, areftemp) IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempicefile after interpweight_2Dcont' ! Copy reftempice_file temperature to all ground levels DO ib=1, kjpindex reftempice(ib, :) = reftempice_file(ib) END DO END SUBROUTINE explicitsnow_read_reftempicefile !! ================================================================================================================================ !! SUBROUTINE : explicitsnow_maxmass !! !>\BRIEF !! !! DESCRIPTION : limits the snow mass below a threshold (maxmass_snow) !! !! !! RECENT CHANGE(S) : adapted for nsnow=12 !! !! MAIN OUTPUT VARIABLE(S): snow, snowdz & snowmelt_from_maxmass !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE explicitsnow_maxmass(kjpindex,snowrho,soilcap,snow,snowdz,snowmelt_from_maxmass) IMPLICIT NONE ! variables pas necessaires a transmettre ? : maxmass_snow, chalfu0 !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density (Kg/m^3) REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity !! 0.2 Output variables REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: snowmelt_from_maxmass !! snowmelt to limit snow mass under maxmass_snow !! 0.3 Local variables INTEGER(i_std) :: ji, jj INTEGER(i_std) :: nsnowstart, nsnowstop !! start and stop index layer INTEGER(i_std) :: locjj REAL(r_std) :: snow_d1k REAL(r_std),DIMENSION(kjpindex,nsnow) :: WSNOWDZ !! snow mass (kg/m2) REAL(r_std),DIMENSION(kjpindex) :: SMASSC !! snow mass cumulated REAL(r_std),DIMENSION(kjpindex,nsnow) :: SMASS !! snow mass cumulated starting from bottom layer REAL(r_std),DIMENSION(kjpindex) :: ZSNOWMELT_XS !! thickness of the layer that has been partially removed REAL(r_std),DIMENSION(kjpindex) :: ZSNOWDZ !! new thickness of the snow layer REAL(r_std),DIMENSION(kjpindex) :: snow_remove !! snow mass of layers that completely disapears ! initialisation ! snowmelt_from_maxmass(:) = zero IF (nsnow.EQ.3) THEN ! snowmelt_from_maxmass can be applied to layers 3 to 2 nsnowstart=3 nsnowstop=2 ELSEIF (nsnow.EQ.12) THEN ! snowmelt_from_maxmass can be applied to layers 9 to 6 because 10, 11 and 12 are thin nsnowstart=9 nsnowstop=6 ELSE WRITE(numout,*) 'explicitsnow_maxmass : nsnow has a non implemented value : ', nsnow STOP 'explicitsnow_maxmass' ENDIF DO ji=1,kjpindex !domain size snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:)) IF(snow(ji).GT.maxmass_snow)THEN IF (snow(ji).GT.(2.5*maxmass_snow)) THEN ! melt much faster for points where accumulation is very very high snow_d1k = 10. * soilcap(ji) / chalfu0 ELSE IF (snow(ji).GT.(1.5*maxmass_snow)) THEN ! melt much faster for points where accumulation is very high snow_d1k = six * soilcap(ji) / chalfu0 ELSE IF (snow(ji).GT.(1.2*maxmass_snow)) THEN ! melt faster for points where accumulation is too high snow_d1k = trois * soilcap(ji) / chalfu0 ELSE ! standard melting snow_d1k = un * soilcap(ji) / chalfu0 ENDIF snowmelt_from_maxmass(ji) = MIN((snow(ji) - maxmass_snow),snow_d1k) ! Calculating the snow accumulation WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) ! WSNOWDZ in kg/m2 SMASSC(ji)= 0.0 DO jj=nsnowstart,nsnowstop,-1 SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj) SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj) ENDDO ! Finding the layer locjj = nsnowstart ! search the layer starting from nsnowstart DO jj=nsnowstart,nsnowstop,-1 IF ((SMASS(ji,jj) .LE. snowmelt_from_maxmass(ji)) .AND. (SMASS(ji,jj-1) .GE. snowmelt_from_maxmass(ji)) ) THEN locjj=jj-1 ENDIF ENDDO ! Calculating the removal of snow depth IF (locjj .EQ. nsnowstart) THEN ! ZSNOWMELT_XS : Epaisseur de la couche qui a disparu partiellement ZSNOWMELT_XS(ji) = snowmelt_from_maxmass(ji)/snowrho(ji,nsnowstart) ZSNOWDZ(ji) = snowdz(ji,nsnowstart) - ZSNOWMELT_XS(ji) snowdz(ji,nsnowstart) = MAX(0.0, ZSNOWDZ(ji)) ELSE IF (locjj .LT. nsnowstart) THEN snow_remove(ji)=0 ! masse de neige des couches qui disparaissent totalement DO jj=nsnowstart,locjj+1,-1 snow_remove(ji)=snow_remove(ji)+snowdz(ji,jj)*snowrho(ji,jj) snowdz(ji,jj)=0 ENDDO ZSNOWMELT_XS(ji) = (snowmelt_from_maxmass(ji) - snow_remove(ji))/snowrho(ji,locjj) ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWMELT_XS(ji) snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji)) ENDIF ENDIF ENDDO END SUBROUTINE explicitsnow_maxmass !! ================================================================================================================================ !! SUBROUTINE : explicitsnow_subli !! !>\BRIEF !! !! DESCRIPTION : sublimation on snow !! !! !! RECENT CHANGE(S) : !! !! MAIN OUTPUT VARIABLE(S): snow, snowdz, subsnownobio, subsinksoil, snowliq, snowtemp, vevapsno !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE explicitsnow_subli(kjpindex,snowrho,snowdz,vevapsno, & snowliq,snowtemp,snow,subsnownobio,subsinksoil) IMPLICIT NONE !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density (Kg/m^3) !! 0.2 Output variables REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Snow liquid content (m) REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...) REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil !! 0.3 Local variables INTEGER(i_std) :: ji, jj INTEGER(i_std) :: nsnowstart, nsnowstop !! start and stop index layer INTEGER(i_std) :: locjj REAL(r_std),DIMENSION(kjpindex,nsnow) :: WSNOWDZ !! snow mass (kg/m2) REAL(r_std),DIMENSION(kjpindex) :: SMASSC !! snow mass cumulated REAL(r_std),DIMENSION(kjpindex,nsnow) :: SMASS !! snow mass cumulated starting from bottom layer REAL(r_std),DIMENSION(kjpindex) :: ZSNOWEVAPS !! thickness of the layer that has been removed by sublimation REAL(r_std),DIMENSION(kjpindex) :: ZSNOWDZ !! new thickness of the snow layer REAL(r_std), DIMENSION (kjpindex) :: snowacc snow(:) = 0.0 DO ji=1,kjpindex !domain size snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:)) ENDDO subsnownobio(:,:) = zero subsinksoil(:) = zero DO ji=1, kjpindex ! domain size IF (vevapsno(ji) .GT. snow(ji)) THEN subsinksoil (ji) = vevapsno(ji) - snow(ji) vevapsno(ji) = snow(ji) snow(ji) = zero snowdz(ji,:) = 0 snowliq(ji,:) = 0 snowtemp(ji,:) = tp_00 ELSE ! Calculating the snow accumulation WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) SMASSC (ji)= 0.0 DO jj=1,nsnow SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj) SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj) ENDDO ! Finding the layer locjj=1 DO jj=1,nsnow-1 IF ((SMASS(ji,jj) .LE. vevapsno(ji)) .AND. (SMASS(ji,jj+1) .GE. vevapsno(ji)) ) THEN locjj=jj+1 ENDIF ENDDO ! Calculating the removal of snow depth IF (locjj .EQ. 1) THEN ZSNOWEVAPS(ji) = vevapsno(ji)/snowrho(ji,1) ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji) snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji)) ELSE IF (locjj .GT. 1) THEN snowacc(ji)=0 DO jj=1,locjj-1 snowacc(ji)=snowacc(ji)+snowdz(ji,jj)*snowrho(ji,jj) snowdz(ji,jj)=0 ENDDO ZSNOWEVAPS(ji) = (vevapsno(ji)-snowacc(ji))/snowrho(ji,locjj) ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWEVAPS(ji) snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji)) ! ELSE ! ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1) ! ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji) ! snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji)) ENDIF ENDIF ENDDO END SUBROUTINE explicitsnow_subli !! ================================================================================================================================ !! SUBROUTINE : explicitsnow_age !! !>\BRIEF !! !! DESCRIPTION : compute snow age for albedo !! !! !! RECENT CHANGE(S) : !! !! MAIN OUTPUT VARIABLE(S): snowage, snownobioage !! !! REFERENCE(S) : !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE explicitsnow_age(kjpindex,snow,precip_snow,precip_rain,frac_snow_nobio, & temp_sol_new,snow_age,snow_nobio_age) IMPLICIT NONE !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjpindex REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature !! 0.2 Output variables REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ... !! 0.3 Local variables INTEGER(i_std) :: ji REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change REAL(r_std), DIMENSION (kjpindex) :: xx !! Temporary DO ji = 1, kjpindex !! 5.1. Snow age on land IF (snow(ji) .LE. zero) THEN snow_age(ji) = zero ELSE snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) & & * EXP(-precip_snow(ji) / snow_trans) ENDIF !! 5.2. Snow age on land ice !! Age of snow on ice: a little bit different because in cold regions, we really !! cannot negect the effect of cold temperatures on snow metamorphism any more. IF ( (frac_snow_nobio(ji,iice) .LE. zero) .OR. (snow(ji) .LE. zero) ) THEN snow_nobio_age(ji,iice) = zero ELSE d_age(ji) = ( snow_nobio_age(ji,iice) + & & (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * & & EXP(-precip_snow(ji) / snow_trans_nobio) - snow_nobio_age(ji,iice) IF (d_age(ji) .GT. 0. ) THEN xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero ) !!-SC xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std ! paramétrisation standard !!-SC Pour ralentir le vieillissemnt (surtout pour les basses températures) : xx(ji) = ( xx(ji) / omg1 ) ** omg2 !! omg1,omg2 values in run.def d_age(ji) = d_age(ji) / (un+xx(ji)) !cdc vieillissement accelere par 2 si il pleut : IF (precip_rain(ji) .GT.0.) THEN d_age(ji) = d_age(ji)*2. ENDIF ENDIF snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero ) ENDIF ENDDO END SUBROUTINE explicitsnow_age END MODULE explicitsnow