MODULE surf_inlandsis_mod IMPLICIT NONE CONTAINS SUBROUTINE surf_inlandsis(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, & rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, & precip_rain, precip_snow, precip_snow_adv, snow_adv, & zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & rugos, snow_cont_air, alb_soil, slope, cloudf, & radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, & tsurf_new, alb1, alb2, alb3, & 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 | ! | Update Etienne Vignon, LMD, Novembre 2020 | ! | | ! +------------------------------------------------------------------------+ ! | ! | 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 | ! | SnoMod: Snow Pack is set up when .T. | ! | reaLBC: Update Bound.Condit.when .T. | ! | | ! | INPUT (via MODULES VARxSV, VARySV, VARtSV) | ! | ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables | ! | | ! | Preprocessing Option: SISVAT PHYSICS | ! | ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^ | ! | # #HY | ! | # #SN: Snow Model | ! | # #BS: Blowing Snow Parameterization | ! +------------------------------------------------------------------------+ 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 IMPLICIT NONE ! +--INTERFACE Variables ! + =================== ! include "dimsoil.h" ! +--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) :: precip_snow_adv, snow_adv !Snow Drift 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,snow_cont_air REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope 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), INTENT(OUT) :: emis_new ! Surface Emissivity REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff 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) :: 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 INTEGER :: i, ig, ikl, isl, isn, nt INTEGER :: gp_outfor, un_outfor REAL, PARAMETER :: f1=0.5 REAL, PARAMETER :: sn_upp=5000.,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 REAL :: lambda ! Par. soil discret. REAL, DIMENSION(nsoilmx), SAVE :: dz1,dz2 ! Soil layer thicknesses !$OMP THREADPRIVATE(dz1,dz2) LOGICAL, SAVE :: firstcall !$OMP THREADPRIVATE(firstcall) ! +--Internal Variables ! + ================== INTEGER :: iso LOGICAL :: file_exists CHARACTER(len=20) :: fichnom !======================================================================== PRINT*, 'je rentre dans inlandsis' zsigma=1000. dt__SV=dtime ! write(*,*)'Start of simulation? ',debut !hj 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 fn_outfor='outfor_SV.dat' END IF ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + INITIALISATION: BEGIN +++ ! + ------------------------- ! + ! + Compute soil discretization (as for LMDZ) ! + ----------------------------------------- IF (firstcall) THEN ! +--Array size klonv=klon knonv=knon write(*,*)'klon',klon,'klonv',klonv,'knon',knon,'nsol',nsol,'nsno',nsno CALL INIT_VARtSV CALL INIT_VARxSV CALL INIT_VARySV eps0SL(:)=0. ! +--Soil layer thickness ! + ----------------------- ! 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 DO ikl=1,knon ! Initialise variables ispiSV(ikl) = 0 iiceSV(ikl) = 0 rusnSV(ikl) = 0. toicSV(ikl) = 0. isnoSV(ikl) = 0. ! # snow layers istoSV(ikl,:) = 0. eta_SV(ikl,:) = 0. TsisSV(ikl,:) = 0. ro__SV(ikl,:) = 0. G1snSV(ikl,:) = 0. G2snSV(ikl,:) = 0. agsnSV(ikl,:) = 0. dzsnSV(ikl,:) = 0. zzsnsv(ikl,:) = 0. BufsSV(ikl) = 0. qsnoSV(ikl) = 0. ! BL snow content zWEcSV(ikl) = 0. dbs_SV(ikl) = 0. dsnbSV(ikl) = 0. esnbSV(ikl) = 0. BrosSV(ikl) = 0. BG1sSV(ikl) = 0. BG2sSV(ikl) = 0. SWS_SV(ikl) = 0. RnofSV(ikl) = 0. ! RunOFF Intensity RRs_SV(ikl) = 0. DDs_SV(ikl) = 0. VVs_SV(ikl) = 0. cld_SV(ikl) = 0. uts_SV(ikl) = 0. ! u*T* arbitrary uqs_SV(ikl) = 0. ! u*q* " uss_SV(ikl) = 0. ! u*s* " LMO_SV(ikl) = 0. ! Set variables LSmask(ikl) = 1 ! Land/Sea Mask isotSV(ikl) = 12 ! Soil Type -> 12= ice iWaFSV(ikl) = 1 ! Soil Drainage eps0SL(ikl )= 1. alb0SV(ikl) = alb_soil(ikl) ! Soil Albedo Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L. Z0h_SV(ikl) = z0h(ikl) ! heat Roughn.L. ! + Soil Upward IR Flux, Water Fluxes, roughness length IRs_SV(ikl) = & -eps0SL(ikl)* StefBo*(temp_air(ikl)**4) ! Upward IR Flux Tsf_SV(ikl) = min(temp_air(ikl),TfSnow) ! + 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 !! Initialise with snow ! G1snSV(ikl,0) = 0. ! [-] ! G2snSV(ikl,0) = 1.6 ! [-] [0.0001 m] ! dzsnSV(ikl,0) = dz_dSV(0) ! [m] ! if (snow(ikl) .GT. 0.) then ! isnoSV(ikl) = 1 ! snow layers ! istoSV(ikl,1:nsno) = 0 ! 0,...,5 : Snow History (see istdSV data) ! eta_SV(ikl,1:nsno) = epsi ! TsisSV(ikl,1:nsno) = tsoil(ikl,1) ! ro__SV(ikl,1:nsno) = 350.0 ! G1snSV(ikl,1:nsno) = 0. ! [-] ! G2snSV(ikl,1:nsno) = 1.6 ! [-] [0.0001 m] ! agsnSV(ikl,1:nsno) = 50. ! [day] ! dzsnSV(ikl,1) = snow(ikl)/max(ro__SV(ikl,1),epsi) ![m] ! ! ecrete si trop de neige: ! IF (snow(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 ! zzsnsv(ikl,1) = dzsnSV(ikl,1) ! Total snow pack thickness ! endif ! Initialise la neige avec un profil de densité prochde des conditions de Dôme C (~10m de neige avec 19 niveaux) (Etienne): 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)) = 0 G2snSV(ikl,1:isnoSV(ikl)) = 1.6 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) =0.100 END DO ! +--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 ! + SISVAT_ini (as for use with MAR, but not computing soil layers) ! + ------------------------------------------------------------- ! write(*,'(/a)') 'Start SISVAT initialization: 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 write(un_outfor,*) 'nsnow - albedo - z0m - z0h , dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46] & & G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]' 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) ! + ^^^^^^^^^^^^^^^^^^^ ^^^^^ zSBLSV = 1000. ! [m] za__SV(ikl) = zsl_height(ikl) ! surface layer height (fisr model level) [m] Ua_min = epsi ! 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] !c #BS dbsnow = -SLussl(i,j,n) ! Erosion !c #BS. *dtPhys *rhT_SV(ikl) /ro_Wat !c #BS dsnbSV(ikl) = snow_adv(ikl) ! min(max(zero,dbsnow) !c #BS. / max(epsi,d_snow),unun) !c #BS dbs_SV(ikl) = snow_cont_air(ikl) !c #BS blowSN(i,j,n) ! [kg/m2] ! +--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) = cdragh(ikl) cdM_SV(ikl) = cdragm(ikl) Us_min = 0.01 us__SV(ikl) = max(Us_min, ustar(ikl)) ram_sv(ikl) = 1./(cdragm(ikl)*max(VV__SV(ikl),eps6)) rah_sv(ikl) = 1./(cdragh(ikl)*max(VV__SV(ikl),eps6)) ! +--Energy Fluxes (INPUT/OUTPUT) ! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^ IF (.not.firstcall) THEN Tsf_SV(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..) runoff_lic(ikl) = RnofSV(ikl)*dtime ! RunOFF: intensity* time step 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) = HLs_sv(ikl)/LHvH2O ! Evaporation snow(ikl) = 0. snowhgt(ikl) = 0. qsnow(ikl) = 0. qsol(ikl) = 0. ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! + ! + Check snow thickness, substract if too thick (commended by etienne: 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. !() DO isn = 1,isnoSV(ikl) sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn) snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn) qsnow(ikl) = qsnow(ikl)+1e03*eta_SV(ikl,isn)*dzsnSV(ikl,isn) END DO ! Etienne: pourquoi ajouter toicSV ici? Pour bilan d'eau? snow(ikl) = sissnow(ikl)+toicSV(ikl) to_ice(ikl) = toicSV(ikl) DO isl = -nsol,0 tsoil(ikl,1-isl) = TsisSV(ikl,isl) ! Soil Temperature 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 ! Albedo NIR alb3(ikl) = alb3sv(ikl) ! Albedo FIR tsurf_new(ikl) =Tsrfsv(ikl) zfra(ikl) = max(min(isnoSV(ikl)-iiceSV(ikl),1),0) qsurf(ikl) = QaT_SV(ikl) emis_new(ikl) = eps0SL(ikl) z0m(ikl) = Z0m_SV(ikl) z0h(ikl) = Z0h_SV(ikl) END DO ! ikl ! write variables in output file IF (ok_outfor) THEN ikl=gp_outfor ! write(un_outfor,*) 'nsnow [-,1], dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46]' ! write(un_outfor,*) 'G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]' write(un_outfor,*) '+++++++++++++++++++++++++++++++++++++++++++++++' write(un_outfor,*) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(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,:) 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 ! INCLUDE "dimsoil.h" 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) PRINT*,'full layers, intermediate layers (seconds)' DO jk=1,nsoilmx rk=jk rk1=jk+.5 rk2=jk-.5 PRINT *,'fz=', & fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 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 ! Piecewis.https://www.lequipe.fr/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 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 IMPLICIT none !====================================================================== ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 ! Objet: Lecture du fichier de conditions initiales pour SISVAT !====================================================================== include "netcdf.inc" ! include "indicesol.h" ! include "dimsoil.h" include "clesphys.h" include "thermcell.h" include "compbl.h" !====================================================================== 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 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) ! IF (.NOT. found) THEN ! PRINT*, 'phyetat0: Le champ est absent' ! rusn(:)=0. ! ENDIF CALL get_field("to_ice",toic,found) IF (.NOT. found) THEN PRINT*, 'phyetat0: Le champ 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 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 IMPLICIT none include "netcdf.inc" ! include "indicesol.h" ! include "dimsoil.h" include "clesphys.h" include "thermcell.h" include "compbl.h" !====================================================================== 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 print*, 'je rentre dans restart inlandsis' 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 print*, 'je call open_restart' CALL open_restartphy(fichnom) print*, 'je sors open_restart' 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