MODULE surf_inlandsis_mod IMPLICIT NONE; PRIVATE PUBLIC surf_inlandsis, get_soil_levels, SISVAT_ini, sisvatetat0, sisvatredem CONTAINS SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, & rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, & precip_rain, precip_snow, & zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & rugos, snow_cont_air, alb_soil, alt, slope, cloudf, & radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, & tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf) ! | | ! | SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow| ! | (INLANDSIS) | ! | SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme) | ! | surface scheme of the Modele Atmospherique Regional (MAR) | ! | Author: Heinz Juergen Punge, LSCE June 2009 | ! | based on the MAR-SISVAT interface by Hubert Gallee | ! | Updated by Etienne Vignon, Cecile Agosta | ! | | ! +------------------------------------------------------------------------+ ! | ! | In the current setup, SISVAT is used only to model the land ice | ! | part of the surface; hence it is called with the compressed variables| ! | from pbl_surface, and only by the surf_landice routine. | ! | | ! | In this interface it is assumed that the partitioning of the soil, | ! | and hence the number of grid points is constant during a simulation, | ! | hence eg. snow properties remain stored in the global SISVAT | ! | variables between the calls and don't need to be handed over as | ! | arguments. When the partitioning is supposed to change, make sure to | ! | update the variables. | ! | | ! | INPUT (via MODULES VARxSV, VARySV, VARtSV ...) | ! | ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables | ! | | ! +------------------------------------------------------------------------+ USE dimphy USE VAR_SV USE VARdSV USE VARxSV USE VARySV USE VARtSV USE VARphy USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx IMPLICIT NONE ! +--Global Variables ! + ================ ! Input Variables for SISVAT INTEGER, INTENT(IN) :: knon INTEGER, INTENT(IN) :: itime REAL, INTENT(IN) :: dtime LOGICAL, INTENT(IN) :: debut ! true if first step LOGICAL, INTENT(IN) :: lafin ! true if last step INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i ! Index Decompression REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cos sol. zenith angle REAL, DIMENSION(klon), INTENT(IN) :: swdown ! REAL, DIMENSION(klon), INTENT(IN) :: lwdown ! REAL, DIMENSION(klon), INTENT(IN) :: albedo_old REAL, DIMENSION(klon), INTENT(IN) :: pexner ! Exner potential REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf REAL, DIMENSION(klon), INTENT(IN) :: rugos REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope REAL, DIMENSION(klon), INTENT(IN) :: alt ! surface elevation REAL, DIMENSION(klon), INTENT(IN) :: cloudf REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh REAL, DIMENSION(klon), INTENT(IN) :: ustar ! friction velocity ! Variables exchanged between LMDZ and SISVAT REAL, DIMENSION(klon), INTENT(IN) :: radsol ! Surface absorbed rad. REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! Tot snow mass [kg/m2] REAL, DIMENSION(klon), INTENT(INOUT) :: zfra ! snwo surface fraction [0-1] REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! Soil Water Content REAL, DIMENSION(klon), INTENT(INOUT) :: z0m ! Momentum Roughn Lgt REAL, DIMENSION(klon), INTENT(INOUT) :: z0h ! Momentum Roughn Lgt ! Output Variables for LMDZ REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! Albedo SW REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo REAL, DIMENSION(klon), INTENT(OUT) :: emis_new ! Surface Emissivity REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff REAL, DIMENSION(klon), INTENT(OUT) :: ffonte ! enthalpy flux due to surface melting REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte ! water flux due to surface melting REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s ! d/dT sens. ht flux REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l ! d/dT latent ht flux REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens ! Sensible ht flux REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat ! Latent heat flux REAL, DIMENSION(klon), INTENT(OUT) :: evap ! Evaporation REAL, DIMENSION(klon), INTENT(OUT) :: erod ! Erosion of surface snow (flux) REAL, DIMENSION(klon), INTENT(OUT) :: agesno ! Snow age (top layer) REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! Surface Humidity ! Specific INLANDIS outputs REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! Total H2O snow[kg/m2] REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! Snow height (m) REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! Snow passed to ice REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! Snow in model (kg/m2) ! +--Internal Variables ! + =================== CHARACTER(len = 20) :: fn_outfor ! Name for output file CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'surf_inlandsis_mod' INTEGER :: i, ig, ikl, isl, isn, nt INTEGER :: gp_outfor, un_outfor REAL, PARAMETER :: f1 = 0.5 REAL, PARAMETER :: sn_upp = 10000., sn_low = 500. REAL, PARAMETER :: sn_add = 400., sn_div = 2. ! snow mass upper,lower limit, ! added mass/division lowest layer REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6 REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3 ! Parameters for drainage ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ ! Tuning ! +... Run Off Parameters ! + 86400*1.5 day ...*25 days (Modif. ETH Camp: 86400*0.3day) ! + (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317) REAL, DIMENSION(klon) :: eps0SL ! surface Emissivity REAL :: zsigma, Ua_min, Us_min, lati REAL, PARAMETER :: cdmax=0.05 REAL :: lambda ! Par. soil discret. REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2 ! Soil layer thicknesses !$OMP THREADPRIVATE(dz1,dz2) LOGICAL, SAVE :: firstcall !$OMP THREADPRIVATE(firstcall) INTEGER :: iso LOGICAL :: file_exists CHARACTER(len = 20) :: fichnom LOGICAL :: is_init_domec ! CA initialization ! dz_profil_15 : 1 m in 15 layers [m] real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, & 0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/) ! mean_temp : mean annual surface temperature [K] real, dimension(klon) :: mean_temp ! mean_dens : mean surface density [kg/m3] real, dimension(klon) :: mean_dens ! lat_scale : temperature lapse rate against latitude [K degree-1] real :: lat_scale ! sh_scale : temperature lapse rate against altitude [K km-1] real :: sh_scale ! variables for density profile ! E0, E1 : exponent real :: E0, E1 ! depth at which 550 kg m-3 is reached [m] real :: z550 ! depths of snow layers real :: depth, snow_depth, distup ! number of initial snow layers integer :: nb_snow_layer ! For density calc. real :: alpha0, alpha1, ln_smb ! theoritical densities [kg m-3] real :: rho0, rho1, rho1_550 ! constants for density profile ! C0, C1 : constant, 0.07 for z <= 550 kg m-3 real, parameter :: C0 = 0.07 real, parameter :: C1 = 0.03 ! rho_i : ice density [kg m-3] real, parameter :: rho_ice = 917. ! E_c : activation energy [J mol-1] real, parameter :: E_c = 60000. ! E_g : activation energy [J mol-1] real, parameter :: E_g = 42400. ! R : gas constant [J mol-1 K-1] real, parameter :: R = 8.3144621 ! + PROGRAM START ! + ----------------------------------------- zsigma = 1000. dt__SV = dtime IF (debut) THEN firstcall = .TRUE. INI_SV = .false. ELSE firstcall = .false. INI_SV = .true. END IF IF (ok_outfor) THEN un_outfor = 51 ! unit number for point output file gp_outfor = 1 ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC fn_outfor = 'outfor_SV.dat' END IF ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + INITIALISATION: BEGIN +++ ! + ----------------------------------------- IF (firstcall) THEN ! +--Array size ! + ----------------------- klonv = klon knonv = knon write(*, *) 'ikl, lon and lat in INLANDSIS' DO ikl = 1, knon i=ikl2i(ikl) write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i) END DO ! +--Variables initizialisation ! + --------------------------- CALL INIT_VARtSV CALL INIT_VARxSV CALL INIT_VARySV ! +--Surface Fall Line Slope ! + ----------------------- IF (SnoMod) THEN DO ikl = 1, knon slopSV(ikl) = slope(ikl) SWf_SV(ikl) = & ! Normalized Decay of the exp(-dt__SV & ! Surficial Water Content / (c1_zuo & !(Zuo and Oerlemans 1996, + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl))))) ! J.Glacio. 42, 305--317) END DO END IF ! +--Soil layer thickness . Compute soil discretization (as for LMDZ) ! + ---------------------------------------------------------------- ! write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx CALL get_soil_levels(dz1, dz2, lambda) lambSV = lambda dz1_SV(1:knon, 1:) = 0. dz2_SV(1:knon, 1:) = 0. DO isl = -nsol, 0 dz_dSV(isl) = 0.5e-3 * dz2(1 - isl) ! Soil layer thickness DO ikl = 1, knon dz1_SV(ikl, isl) = dz1(1 - isl) !1.e-3* dz2_SV(ikl, isl) = dz2(1 - isl) !1.e-3* END DO END DO ! Set variables ! ============= DO ikl = 1, knon ! LSmask : Land/Sea Mask LSmask(ikl) = 1 ! isotSV : Soil Type -> 12 = ice isotSV(ikl) = 12 ! iWaFSV : Soil Drainage (1,0)=(y,n) iWaFSV(ikl) = 1 ! eps0SL : Surface Emissivity eps0SL(ikl) = 1. ! alb0SV : Soil Albedo alb0SV(ikl) = alb_soil(ikl) ! Tsf_SV : Surface Temperature, must be bellow freezing Tsf_SV(ikl) = min(temp_air(ikl), TfSnow) END DO ! +--Initialization of soil and snow variables in case startsis is not read ! + ---------------------------------------------------------------------- is_init_domec=.FALSE. IF (is_init_domec) THEN ! Coarse initilization inspired from vertcical profiles at Dome C, ! Antarctic Plateaui (10m of snow, 19 levels) DO ikl = 1,knon ! + Soil DO isl = -nsol,0 TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2) !temp_air(ikl) !tsoil(ikl,1-isl) Soil Temperature !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2) eta_SV(ikl,isl) = epsi !etasoil(ikl,1-isl) Soil Water[m3/m3] ro__SV(ikl,isl) = rhoIce !rosoil(ikl,1-isl) volumic mass END DO ! Snow isnoSV(ikl) = 19 istoSV(ikl, 1:isnoSV(ikl)) = 100 ro__SV(ikl, 1:isnoSV(ikl)) = 350. eta_SV(ikl, 1:isnoSV(ikl)) = epsi TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2) G1snSV(ikl, 1:isnoSV(ikl)) = 99. G2snSV(ikl, 1:isnoSV(ikl)) = 2. agsnSV(ikl, 1:isnoSV(ikl)) = 50. dzsnSV(ikl, 19) = 0.015 dzsnSV(ikl, 18) = 0.015 dzsnSV(ikl, 17) = 0.020 dzsnSV(ikl, 16) = 0.030 dzsnSV(ikl, 15) = 0.040 dzsnSV(ikl, 14) = 0.060 dzsnSV(ikl, 13) = 0.080 dzsnSV(ikl, 12) = 0.110 dzsnSV(ikl, 11) = 0.150 dzsnSV(ikl, 10) = 0.200 dzsnSV(ikl, 9) = 0.300 dzsnSV(ikl, 8) = 0.420 dzsnSV(ikl, 7) = 0.780 dzsnSV(ikl, 6) = 1.020 dzsnSV(ikl, 5) = 0.980 dzsnSV(ikl, 4) = 1.020 dzsnSV(ikl, 3) = 3.970 dzsnSV(ikl, 2) = 1.020 dzsnSV(ikl, 1) = 1.020 END DO ELSE ! Initilialisation with climatological temperature and density ! profiles as in MAR. Methodology developed by Cecile Agosta ! initialize with 0., for unused snow layers dzsnSV = 0. G1snSV = 0. G2snSV = 0. istoSV = 0 TsisSV = 0. ! initialize mean variables (unrealistic) mean_temp = TfSnow mean_dens = 300. ! loop on grid cells DO ikl = 1, knon lati=rlat(ikl2i(ikl)) ! approximations for mean_temp and mean_dens ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1) ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1 ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica, ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed. if (lati > 60.) then ! CA todo : add longitude bounds ! Greenland mean temperature : function of altitude and latitude ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1 lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9 lat_scale = max(min(lat_scale, 0.9), 0.75) ! sh_scale equals the environmental lapse rate : 6.5 °C km-1 sh_scale = 6.5 mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.) ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051 mean_dens(ikl) = 315. else if (lati < -60.) then ! Antarctica mean temperature : function of altitude and latitude ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1 lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3 lat_scale = max(min(lat_scale, 1.3), 0.6) ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1 sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5 sh_scale = max(min(sh_scale, 9.8), 6.5) mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.) ! Antarctica surface density : function of mean annual temperature ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013) ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m. ! Weinhart et al 2020 https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201 ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l. ! Dome C : st_ant_param(3233, -75.1) = -47.7 ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7 mean_dens(ikl) = (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450. mean_dens(ikl) = min(450., max(320., mean_dens(ikl))) else ! write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica' mean_dens(ikl) =350. mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2) !abort_message='temperature initialization is only defined for Greenland and Antarctica' !CALL abort_physic(modname,abort_message,1) end if ! mean_temp is defined for ice ground only mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2) ! Soil layers ! =========== DO isl = -nsol, 0 ! TsisSV : Temperature [K] TsisSV(ikl, isl) = mean_temp(ikl) ! eta_SV : Soil Water [m3/m3] eta_SV(ikl, isl) = epsi ! ro__SV : Volumic Mass [kg/m3] ro__SV(ikl, isl) = rhoIce END DO ! Snow layers ! =========== ! snow_depth : initial snow depth snow_depth = 20. ! nb_snow_layer : initial nb of snow layers nb_snow_layer = 15 ! isnoSV : total nb of snow layers isnoSV(ikl) = nb_snow_layer ! depth : depth of each layer depth = snow_depth do isl = 1, nb_snow_layer ! dzsnSV : snow layer thickness dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1)) ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical G1snSV(ikl, isl) = 99. ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size G2snSV(ikl, isl) = 3. agsnSV(ikl, isl) = 0. istoSV(ikl, isl) = 0 ! eta_SV : Liquid Water Content [m3/m3] eta_SV(ikl, isl) = 0. ! distance to surface depth = depth - dzsnSV(ikl,isl) / 2. distup = min(1., max(0., depth / snow_depth)) ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom) TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2 ! firn density : densification formulas from : ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/) ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306) ! Integration of the steady state equation ! ln_smb approximated as a function of temperature ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.) ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1 alpha0 = max(1.435 - 0.151 * ln_smb, 0.25) alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25) E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0 E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1 z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0 rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice if (depth <= z550) then ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice else ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice end if depth = depth - dzsnSV(ikl,isl) / 2. end do END DO END IF ! + Numerics paramaters, SISVAT_ini ! + ---------------------- CALL SISVAT_ini(knon) ! +--Read restart file ! + ================================================= INQUIRE(FILE = "startsis.nc", EXIST = file_exists) IF (file_exists) THEN CALL sisvatetat0("startsis.nc", ikl2i) END IF ! +--Output ascii file ! + ================================================= ! open output file IF (ok_outfor) THEN open(unit = un_outfor, status = 'replace', file = fn_outfor) ikl = gp_outfor ! index sur la grille land ice write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl)) write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] & & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]' END IF END IF ! firstcall ! + ! + +++ INITIALISATION: END +++ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + READ FORCINGS ! + ------------------------ ! + Update Forcings for SISVAT given by the LMDZ model. ! + DO ikl = 1, knon ! +--Atmospheric Forcing (INPUT) ! + ^^^^^^^^^^^^^^^^^^^ ^^^^^ za__SV(ikl) = zsl_height(ikl) ! surface layer height (fisr model level) [m] Ua_min = 0.2 * sqrt(za__SV(ikl)) ! VV__SV(ikl) = max(Ua_min, wind_velo(ikl)) ! Wind velocity [m/s] TaT_SV(ikl) = temp_air(ikl) ! BL top Temperature [K] ExnrSV(ikl) = pexner(ikl) ! Exner potential rhT_SV(ikl) = dens_air(ikl) ! Air density QaT_SV(ikl) = spechum(ikl) ! Specific humidity ps__SV(ikl) = ps(ikl) ! surface pressure [Pa] p1l_SV(ikl) = p1lay(ikl) ! lowest atm. layer press[Pa] ! +--Surface properties ! + ^^^^^^^^^^^^^^^^^^ Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L. Z0h_SV(ikl) = z0h(ikl) ! Moment.Roughn.L. ! +--Energy Fluxes (INPUT) ! + ^^^^^^^^^^^^^ ^^^^^ coszSV(ikl) = max(czemin, rmu0(ikl)) ! cos(zenith.Dist.) sol_SV(ikl) = swdown(ikl) ! downward Solar IRd_SV(ikl) = lwdown(ikl) ! downward IR rsolSV(ikl) = radsol(ikl) ! surface absorbed rad. ! +--Water Fluxes (INPUT) ! + ^^^^^^^^^^^^^ ^^^^^ drr_SV(ikl) = precip_rain(ikl) ! Rain fall rate [kg/m2/s] dsn_SV(ikl) = precip_snow(ikl) ! Snow fall rate [kg/m2/s] ! #BS dbs_SV(ikl) = blowSN(i,j,n) ! dbs_SV = Maximum potential erosion amount [kg/m2] ! => Upper bound for eroded snow mass ! uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f) ! #BS if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then ! #BS dsnbSV(ikl) =1.0-min(qsHY(i,j,kB) !BS neglib. at kb ~100 magl) ! #BS. /max(qshy(i,j,mz),eps9),unun) ! #BS dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl)) ! #BS dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl))) ! #BS else ! #BS dsnbSV(ikl) = 0. ! #BS endif ! dsnbSV is the drift fraction of deposited snow updated in sisvat.f ! will be used for characterizing the Buffer Layer ! (see update of Bros_N, G1same, G2same, zroOLD, zroNEW) ! #BS if(n==1) qbs_HY(i,j) = dsnbSV(ikl) qsnoSV(ikl) = snow_cont_air(ikl) ! +--Soil/BL (INPUT) ! + ^^^^^^^ ^^^^^ alb0SV(ikl) = alb_soil(ikl) ! Soil background Albedo AcoHSV(ikl) = AcoefH(ikl) BcoHSV(ikl) = BcoefH(ikl) AcoQSV(ikl) = AcoefQ(ikl) BcoQSV(ikl) = BcoefQ(ikl) cdH_SV(ikl) = min(cdragh(ikl),cdmax) cdM_SV(ikl) = min(cdragm(ikl),cdmax) rcdmSV(ikl) = sqrt(cdM_SV(ikl)) Us_min = 0.01 us__SV(ikl) = max(Us_min, ustar(ikl)) ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6)) rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6)) ! +--Energy Fluxes (INPUT/OUTPUT) ! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^ !IF (.not.firstcall) THEN Tsrfsv(ikl) = tsurf(ikl) !hj 12 03 2010 cld_SV(ikl) = cloudf(ikl) ! Cloudiness !END IF END DO ! ! + +++ READ FORCINGS: END +++ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! +--SISVAT EXECUTION ! + ---------------- call INLANDSIS(SnoMod, BloMod, 1) ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + RETURN RESULTS ! + -------------- ! + Return (compressed) SISVAT variables to LMDZ ! + DO ikl = 1, knon ! use only 1:knon (actual ice sheet..) dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der. dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der. fluxsens(ikl) = HSs_sv(ikl) ! HS fluxlat(ikl) = HLs_sv(ikl) ! HL evap(ikl) = -1*HLs_sv(ikl) / LHvH2O ! Evaporation erod(ikl) = 0. IF (BloMod) THEN ! + Blowing snow ! SLussl(i,j,n)= 0. ! #BS SLussl(i,j,n)= !Effective erosion ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl)) !~u*qs* from previous time step ! #BS blowSN(i,j,n)= dt*uss_SV(ikl) !New max. pot. Erosion [kg/m2] ! #BS. *rhT_SV(ikl) !(further bounded in sisvat_bsn.f) ! #BS erprev(i,j,n) = dbs_Er(ikl)/dt__SV erod(ikl) = dbs_Er(ikl) / dt__SV ENDIF ! + Check snow thickness, substract if too thick, add if too thin sissnow(ikl) = 0. !() DO isn = 1, isnoSV(ikl) sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn) END DO IF (sissnow(ikl) .LE. sn_low) THEN !add snow IF (isnoSV(ikl).GE.1) THEN dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi) toicSV(ikl) = toicSV(ikl) - sn_add ELSE write(*, *) 'Attention, bare ice... point ', ikl isnoSV(ikl) = 1 istoSV(ikl, 1) = 0 ro__SV(ikl, 1) = 350. dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi) ! 1. eta_SV(ikl, 1) = epsi TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2) G1snSV(ikl, 1) = 0. G2snSV(ikl, 1) = 0.3 agsnSV(ikl, 1) = 10. toicSV(ikl) = toicSV(ikl) - sn_add END IF END IF IF (sissnow(ikl) .ge. sn_upp) THEN !thinnen snow layer below dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div END IF sissnow(ikl) = 0. qsnow(ikl) = 0. snow(ikl) = 0. snowhgt(ikl) = 0. DO isn = 1, isnoSV(ikl) sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn) snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn) ! Etienne: check calc qsnow qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn) END DO zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0) ! Etienne: comment following line ! snow(ikl) = sissnow(ikl)+toicSV(ikl) snow(ikl) = sissnow(ikl) to_ice(ikl) = toicSV(ikl) runoff_lic(ikl) = RnofSV(ikl) ! RunOFF: intensity (flux due to melting + liquid precip) fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing ffonte(ikl)=fqfonte(ikl)*Lf_H2O qsol(ikl) = 0. DO isl = -nsol, 0 tsoil(ikl, 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature ! Etienne: check calc qsol qsol(ikl) = qsol(ikl) & + eta_SV(ikl, isl) * dz_dSV(isl) END DO agesno(ikl) = agsnSV(ikl, isnoSV(ikl)) ! [day] alb1(ikl) = alb1sv(ikl) ! Albedo VIS ! alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl) & ! & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1 alb2(ikl)=alb2sv(ikl) ! Albedo NIR alb3(ikl) = alb3sv(ikl) ! Albedo FIR ! 6 band Albedo alb6(ikl,:)=alb6sv(ikl,:) tsurf_new(ikl) = Tsrfsv(ikl) qsurf(ikl) = QsT_SV(ikl) emis_new(ikl) = eps0SL(ikl) z0m(ikl) = Z0m_SV(ikl) z0h(ikl) = Z0h_SV(ikl) END DO IF (ok_outfor) THEN ikl= gp_outfor write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++' write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl) write(un_outfor, *) dzsnSV(ikl, :) write(un_outfor, *) TsisSV(ikl, :) write(un_outfor, *) ro__SV(ikl, :) write(un_outfor, *) eta_SV(ikl, :) write(un_outfor, *) G1snSV(ikl, :) write(un_outfor, *) G2snSV(ikl, :) write(un_outfor, *) agsnSV(ikl, :) write(un_outfor, *) istoSV(ikl, :) write(un_outfor, *) DOPsnSV(ikl, :) ENDIF ! + ----------------------------- ! + END --- RETURN RESULTS ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (lafin) THEN fichnom = "restartsis.nc" CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat) IF (ok_outfor) THEN close(unit = un_outfor) END IF END IF END SUBROUTINE surf_inlandsis !======================================================================= SUBROUTINE get_soil_levels(dz1, dz2, lambda) ! ====================================================================== ! Routine to compute the vertical discretization of the soil in analogy ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case ! of SISVAT, therefore it's needed here. ! USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root USE mod_phys_lmdz_para USE VAR_SV USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1 REAL, INTENT(OUT) :: lambda !----------------------------------------------------------------------- ! Depthts: ! -------- REAL fz, rk, fz1, rk1, rk2 REAL min_period, dalph_soil INTEGER ierr, jk fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.) ! write(*,*)'Start soil level computation' !----------------------------------------------------------------------- ! Calculation of some constants ! NB! These constants do not depend on the sub-surfaces !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ground levels ! grnd=z/l where l is the skin depth of the diurnal cycle: !----------------------------------------------------------------------- min_period = 1800. ! en secondes dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ. ! !$OMP MASTER ! IF (is_mpi_root) THEN ! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr) ! IF (ierr == 0) THEN ! Read file only if it exists ! READ(99,*) min_period ! READ(99,*) dalph_soil ! PRINT*,'Discretization for the soil model' ! PRINT*,'First level e-folding depth',min_period, & ! ' dalph',dalph_soil ! CLOSE(99) ! END IF ! ENDIF ! !$OMP END MASTER ! CALL bcast(min_period) ! CALL bcast(dalph_soil) ! la premiere couche represente un dixieme de cycle diurne fz1 = SQRT(min_period / 3.14) DO jk = 1, nsoilmx rk1 = jk rk2 = jk - 1 dz2(jk) = fz(rk1) - fz(rk2) ENDDO DO jk = 1, nsoilmx - 1 rk1 = jk + .5 rk2 = jk - .5 dz1(jk) = 1. / (fz(rk1) - fz(rk2)) ENDDO lambda = fz(.5) * dz1(1) DO jk = 1, nsoilmx rk = jk rk1 = jk + .5 rk2 = jk - .5 ENDDO END SUBROUTINE get_soil_levels !=========================================================================== SUBROUTINE SISVAT_ini(knon) !C +------------------------------------------------------------------------+ !C | MAR SISVAT_ini Jd 11-10-2007 MAR | !C | SubRoutine SISVAT_ini generates non time dependant SISVAT parameters | !C +------------------------------------------------------------------------+ !C | PARAMETERS: klonv: Total Number of columns = | !C | ^^^^^^^^^^ = Total Number of continental grid boxes | !C | X Number of Mosaic Cell per grid box | !C | | !C | INPUT: dt__SV : Time Step [s] | !C | ^^^^^ dz_dSV : Layer Thickness [m] | !C | | !C | OUTPUT: [-] | !C | ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] | !C | etamSV : Soil Minimum Humidity [m3/m3] | !C | (based on a prescribed Soil Relative Humidity) | !C | s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. | !C | s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. | !C | aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] | !C | bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] | !C | dzsnSV(0): Soil first Layer Thickness [m] | !C | dzmiSV : Distance between two contiguous levels [m] | !C | dz78SV : 7/8 (Layer Thickness) [m] | !C | dz34SV : 3/4 (Layer Thickness) [m] | !C | dz_8SV : 1/8 (Layer Thickness) [m] | !C | dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] | !C | dtz_SV : dt/dz [s/m] | !C | OcndSV : Swab Ocean / Soil Ratio [-] | !C | Implic : Implicit Parameter (0.5: Crank-Nicholson) | !C | Explic : Explicit Parameter = 1.0 - Implic | !C | | !C | # OPTIONS: #ER: Richards Equation is not smoothed | !C | # ^^^^^^^ #kd: De Ridder Discretization | !C | # #SH: Hapex-Sahel Values ! !C | | !C +------------------------------------------------------------------------+ ! ! !C +--Global Variables !C + ================ USE dimphy USE VARphy USE VAR_SV USE VARdSV USE VAR0SV USE VARxSV USE VARtSV USE VARxSV USE VARySV IMPLICIT NONE !C +--Arguments !C + ================== INTEGER, INTENT(IN) :: knon !C +--Internal Variables !C + ================== INTEGER :: ivt, ist, ikl, isl, isn, ikh INTEGER :: misl_2, nisl_2 REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2 REAL, PARAMETER :: RHsMin = 0.001 ! Min.Soil Relative Humidity REAL :: PsiMax ! Max.Soil Water Potential REAL :: a_Khyd, b_Khyd ! Water conductivity !c #WR REAL :: Khyd_x,Khyd_y !C +--Non Time Dependant SISVAT parameters !C + ==================================== !C +--Soil Discretization !C + ------------------- !C +--Numerical Scheme Parameters !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Implic = 0.75 ! 0.5 <==> Crank-Nicholson Explic = 1.00 - Implic ! !C +--Soil/Snow Layers Indices !C + ^^^^^^^^^^^^^^^^^^^^^^^^ DO isl = -nsol, 0 islpSV(isl) = isl + 1 islpSV(isl) = min(islpSV(isl), 0) islmSV(isl) = isl - 1 islmSV(isl) = max(-nsol, islmSV(isl)) END DO DO isn = 1, nsno isnpSV(isn) = isn + 1 isnpSV(isn) = min(isnpSV(isn), nsno) END DO !C +--Soil Layers Thicknesses !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels! !c #kd IF (nsol.gt.4) THEN !c #kd DO isl=-5,-nsol,-1 !c #kd dz_dSV(isl)= 1. !c #kd END DO !c #kd END IF ! ! IF (nsol.ne.4) THEN ! DO isl= 0,-nsol,-1 ! misl_2 = -mod(isl,2) ! nisl_2 = -isl/2 ! dz_dSV(isl)=(((1-misl_2) * 0.001 ! . + misl_2 * 0.003) * 10**(nisl_2)) * 4. !C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm ! !c +SH dz_dSV(isl)=(((1-misl_2) * 0.001 !c +SH. + misl_2 * 0.003) * 10**(nisl_2)) * 1. ! !c #05 dz_dSV(isl)=(((1-misl_2) * 0.001 !c #05. + misl_2 * 0.008) * 10**(nisl_2)) * 0.5 ! END DO ! dz_dSV(0) = 0.001 ! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004 ! END IF zz_dSV = 0. DO isl = -nsol, 0 dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl))) dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl) dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl) dtz_SV(isl) = dt__SV / dz_dSV(isl) dtz_SV2(isl) = 1. / dz_dSV(isl) dz78SV(isl) = 0.875 * dz_dSV(isl) dz34SV(isl) = 0.750 * dz_dSV(isl) dz_8SV(isl) = 0.125 * dz_dSV(isl) dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl)) & & + 0.750 * dz_dSV(isl) & & + 0.125 * dz_dSV(islpSV(isl)) zz_dSV = zz_dSV + dz_dSV(isl) END DO DO ikl = 1, knon !v dzsnSV(ikl, 0) = dz_dSV(0) END DO !C +--Conversion to a 50 m Swab Ocean Discretization !C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ OcndSV = 0. DO isl = -nsol, 0 OcndSV = OcndSV + dz_dSV(isl) END DO OcndSV = 50. / OcndSV !C +--Secondary Soil Parameters !C + ------------------------------- DO ist = 0, nsot rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6 ! Soil Contrib. to (ro c)_s s1__SV(ist) = bCHdSV(ist) & ! Factor of (eta)**(b+2) & * psidSV(ist) * Ks_dSV(ist) & ! in DR97, Eqn.(3.36) & / (etadSV(ist)**(bCHdSV(ist) + 3.)) ! s2__SV(ist) = Ks_dSV(ist) & ! Factor of (eta)**(2b+3) & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.)) ! in DR97, Eqn.(3.35) !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity) !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Psimax = -(log(RHsMin)) / 7.2E-5 ! DR97, Eqn 3.15 Inversion etamSV(ist) = etadSV(ist) & & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist))) END DO etamSV(12) = 0. !C +--Piecewise Hydraulic Conductivity Profiles !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ist = 0, nsot d__eta = etadSV(ist) / nkhy eta__1 = 0. eta__2 = d__eta DO ikh = 0, nkhy Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35) & * (eta__1 **(2. * bCHdSV(ist) + 3.)) ! Khyd_2 = s2__SV(ist) &! & * (eta__2 **(2. * bCHdSV(ist) + 3.)) ! a_Khyd = (Khyd_2 - Khyd_1) / d__eta ! b_Khyd = Khyd_1 - a_Khyd * eta__1 ! !c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd ! !c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd ! aKdtSV(ist, ikh) = a_Khyd * dt__SV ! bKdtSV(ist, ikh) = b_Khyd * dt__SV ! eta__1 = eta__1 + d__eta eta__2 = eta__2 + d__eta END DO END DO return END SUBROUTINE SISVAT_ini !*************************************************************************** SUBROUTINE sisvatetat0 (fichnom, ikl2i) USE clesphys_mod_h USE dimphy USE mod_grid_phy_lmdz USE mod_phys_lmdz_para USE iostart USE VAR_SV USE VARdSV USE VARxSV USE VARtSV USE indice_sol_mod USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx USE clesphys_mod_h USE compbl_mod_h IMPLICIT none !====================================================================== ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 ! Objet: Lecture du fichier de conditions initiales pour SISVAT !====================================================================== CHARACTER(LEN = *) :: fichnom INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i REAL, DIMENSION(klon) :: rlon REAL, DIMENSION(klon) :: rlat ! les variables globales ecrites dans le fichier restart REAL, DIMENSION(klon) :: isno REAL, DIMENSION(klon) :: ispi REAL, DIMENSION(klon) :: iice REAL, DIMENSION(klon) :: rusn REAL, DIMENSION(klon, nsno) :: isto REAL, DIMENSION(klon, nsismx) :: Tsis REAL, DIMENSION(klon, nsismx) :: eta REAL, DIMENSION(klon, nsismx) :: ro REAL, DIMENSION(klon, nsno) :: dzsn REAL, DIMENSION(klon, nsno) :: G1sn REAL, DIMENSION(klon, nsno) :: G2sn REAL, DIMENSION(klon, nsno) :: agsn REAL, DIMENSION(klon) :: toic INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts CHARACTER (len = 2) :: str2 LOGICAL :: found errT = 0 errro = 0 erreta = 0 errdz = 0 snopts = 0 ! Ouvrir le fichier contenant l'etat initial: CALL open_startphy(fichnom) ! Lecture des latitudes, longitudes (coordonnees): CALL get_field("latitude", rlat, found) CALL get_field("longitude", rlon, found) CALL get_field("n_snows", isno, found) IF (.NOT. found) THEN PRINT*, 'phyetat0: Le champ <n_snows> est absent' PRINT *, 'fichier startsisvat non compatible avec sisvatetat0' ENDIF CALL get_field("n_ice_top", ispi, found) CALL get_field("n_ice", iice, found) CALL get_field("surf_water", rusn, found) CALL get_field("to_ice", toic, found) IF (.NOT. found) THEN PRINT*, 'phyetat0: Le champ <to_ice> est absent' toic(:) = 0. ENDIF DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("AGESNOW" // str2, & agsn(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("DZSNOW" // str2, & dzsn(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("G2SNOW" // str2, & G2sn(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("G1SNOW" // str2, & G1sn(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("ETA" // str2, & eta(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("RO" // str2, & ro(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("TSS" // str2, & Tsis(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL get_field("HISTORY" // str2, & isto(:, isn), found) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO write(*, *)'Read ', fichnom, ' finished!!' !********************************************************************************* ! Compress restart file variables for SISVAT DO ikl = 1, klon i = ikl2i(ikl) IF (i > 0) THEN isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay. ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay. iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay. DO isl = -nsol, 0 ro__SV(ikl, isl) = ro(i, nsno + 1 - isl) ! eta_SV(ikl, isl) = eta(i, nsno + 1 - isl) ! Soil Humidity !hjp 15/10/2010 IF (eta_SV(ikl, isl) <= 1.e-6) THEN !hj check eta_SV(ikl, isl) = 1.e-6 ENDIF TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl) ! Soil Temperature IF (TsisSV(ikl, isl) <= 1.) THEN !hj check ! errT=errT+1 TsisSV(ikl, isl) = 273.15 - 0.2 ! Etienne: negative temperature since soil is ice ENDIF END DO write(*, *)'Copy histo', ikl DO isn = 1, isnoSV(ikl) !nsno snopts = snopts + 1 IF (isto(i, isn) > 10.) THEN !hj check write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn) isto(i, isn) = 1. ENDIF istoSV(ikl, isn) = INT(isto(i, isn)) ! Snow History ro__SV(ikl, isn) = ro(i, isn) ! [kg/m3] eta_SV(ikl, isn) = eta(i, isn) ! [m3/m3] TsisSV(ikl, isn) = Tsis(i, isn) ! [K] IF (TsisSV(ikl, isn) <= 1.) THEN !hj check errT = errT + 1 TsisSV(ikl, isn) = TsisSV(ikl, 0) ENDIF IF (TsisSV(ikl, isn) <= 1.) THEN !hj check TsisSV(ikl, isn) = 263.15 ENDIF IF (eta_SV(ikl, isn) < 1.e-9) THEN !hj check eta_SV(ikl, isn) = 1.e-6 erreta = erreta + 1 ENDIF IF (ro__SV(ikl, isn) <= 10.) THEN !hj check ro__SV(ikl, isn) = 11. errro = errro + 1 ENDIF write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn) G1snSV(ikl, isn) = G1sn(i, isn) ! [-] [-] G2snSV(ikl, isn) = G2sn(i, isn) ! [-] [0.0001 m] dzsnSV(ikl, isn) = dzsn(i, isn) ! [m] agsnSV(ikl, isn) = agsn(i, isn) ! [day] END DO rusnSV(ikl) = rusn(i) ! Surficial Water toicSV(ikl) = toic(i) ! bilan snow to ice END IF END DO END SUBROUTINE sisvatetat0 !====================================================================== SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat) !====================================================================== ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 ! Objet: Ecriture de l'etat de redemarrage pour SISVAT !====================================================================== USE compbl_mod_h USE mod_grid_phy_lmdz USE mod_phys_lmdz_para USE iostart USE VAR_SV USE VARxSV USE VARySV !hj tmp 12 03 2010 USE VARtSV USE indice_sol_mod USE dimphy USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx IMPLICIT none !====================================================================== CHARACTER(LEN = *) :: fichnom INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i REAL, DIMENSION(klon), INTENT(IN) :: rlon REAL, DIMENSION(klon), INTENT(IN) :: rlat ! les variables globales ecrites dans le fichier restart REAL, DIMENSION(klon) :: isno REAL, DIMENSION(klon) :: ispi REAL, DIMENSION(klon) :: iice REAL, DIMENSION(klon, nsnowmx) :: isto REAL, DIMENSION(klon, nsismx) :: Tsis REAL, DIMENSION(klon, nsismx) :: eta REAL, DIMENSION(klon, nsnowmx) :: dzsn REAL, DIMENSION(klon, nsismx) :: ro REAL, DIMENSION(klon, nsnowmx) :: G1sn REAL, DIMENSION(klon, nsnowmx) :: G2sn REAL, DIMENSION(klon, nsnowmx) :: agsn REAL, DIMENSION(klon) :: IRs REAL, DIMENSION(klon) :: LMO REAL, DIMENSION(klon) :: rusn REAL, DIMENSION(klon) :: toic REAL, DIMENSION(klon) :: Bufs REAL, DIMENSION(klon) :: alb1, alb2, alb3 INTEGER isl, ikl, i, isn, ierr CHARACTER (len = 2) :: str2 INTEGER :: pass isno(:) = 0 ispi(:) = 0 iice(:) = 0 IRs(:) = 0. LMO(:) = 0. eta(:, :) = 0. Tsis(:, :) = 0. isto(:, :) = 0 ro(:, :) = 0. G1sn(:, :) = 0. G2sn(:, :) = 0. dzsn(:, :) = 0. agsn(:, :) = 0. rusn(:) = 0. toic(:) = 0. Bufs(:) = 0. alb1(:) = 0. alb2(:) = 0. alb3(:) = 0. !*************************************************************************** ! Uncompress SISVAT output variables for storage DO ikl = 1, klon i = ikl2i(ikl) IF (i > 0) THEN isno(i) = 1. * isnoSV(ikl) ! Nb Snow/Ice Lay. ispi(i) = 1. * ispiSV(ikl) ! Nb Supr.Ice Lay. iice(i) = 1. * iiceSV(ikl) ! Nb Ice Lay. ! IRs(i) = IRs_SV(ikl) ! LMO(i) = LMO_SV(ikl) DO isl = -nsol, 0 ! eta(i, nsno + 1 - isl) = eta_SV(ikl, isl) ! Soil Humidity Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature ro(i, nsno + 1 - isl) = ro__SV(ikl, isl) ! [kg/m3] END DO DO isn = 1, nsno isto(i, isn) = 1. * istoSV(ikl, isn) ! Snow History ro(i, isn) = ro__SV(ikl, isn) ! [kg/m3] eta(i, isn) = eta_SV(ikl, isn) ! [m3/m3] Tsis(i, isn) = TsisSV(ikl, isn) ! [K] G1sn(i, isn) = G1snSV(ikl, isn) ! [-] [-] G2sn(i, isn) = G2snSV(ikl, isn) ! [-] [0.0001 m] dzsn(i, isn) = dzsnSV(ikl, isn) ! [m] agsn(i, isn) = agsnSV(ikl, isn) ! [day] END DO rusn(i) = rusnSV(ikl) ! Surficial Water toic(i) = toicSV(ikl) ! to ice alb1(i) = alb1sv(ikl) alb2(i) = alb2sv(ikl) alb3(i) = alb3sv(ikl) ! Bufs(i) = BufsSV(ikl) END IF END DO CALL open_restartphy(fichnom) DO pass = 1, 2 CALL put_field(pass, "longitude", & "Longitudes de la grille physique", rlon) CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat) CALL put_field(pass, "n_snows", "number of snow/ice layers", isno) CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi) CALL put_field(pass, "n_ice", "number of ice layers", iice) CALL put_field(pass, "IR_soil", "Soil IR flux", IRs) CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO) CALL put_field(pass, "surf_water", "Surficial water", rusn) CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs) CALL put_field(pass, "alb_1", "albedo sw", alb1) CALL put_field(pass, "alb_2", "albedo nIR", alb2) CALL put_field(pass, "alb_3", "albedo fIR", alb3) CALL put_field(pass, "to_ice", "Snow passed to ice", toic) DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "AGESNOW" // str2, & "Age de la neige layer No." // str2, & agsn(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "DZSNOW" // str2, & "Snow/ice thickness layer No." // str2, & dzsn(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "G2SNOW" // str2, & "Snow Property 2, layer No." // str2, & G2sn(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "G1SNOW" // str2, & "Snow Property 1, layer No." // str2, & G1sn(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "ETA" // str2, & "Soil/snow water content layer No." // str2, & eta(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx !nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "RO" // str2, & "Snow density layer No." // str2, & ro(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsismx IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "TSS" // str2, & "Soil/snow temperature layer No." // str2, & Tsis(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO DO isn = 1, nsno IF (isn.LE.99) THEN WRITE(str2, '(i2.2)') isn CALL put_field(pass, "HISTORY" // str2, & "Snow history layer No." // str2, & isto(:, isn)) ELSE PRINT*, "Trop de couches" CALL abort ENDIF ENDDO CALL enddef_restartphy ENDDO CALL close_restartphy END SUBROUTINE sisvatredem END MODULE surf_inlandsis_mod