Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 4082) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 4119) @@ -23,5 +23,5 @@ use aerosol_mod, only: i_haze, haze_prof use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & - dryness + dryness, qsurfyear, phisfibed use comdiurn_h, only: coslat, sinlat, coslon, sinlon use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen @@ -29,4 +29,5 @@ use geometry_mod, only: latitude, longitude, cell_area, & cell_area_for_lonlat_outputs + use control_mod, only: iphysiq USE comgeomfi_h, only: totarea, totarea_planet USE tracer_h, only: noms, mmol, radius, rho_q, qext, & @@ -254,8 +255,9 @@ ! ---------------------- integer,save :: day_ini ! Initial date of the run (sol since Ls=0). + real,save :: time_phys ! Initial time of day of the run integer,save :: icount ! Counter of calls to physiq during the run. !$OMP THREADPRIVATE(day_ini,icount) - !Pluto specific + ! Pluto specific REAL,save :: acond,bcond REAL,save :: tcond1p4Pa @@ -264,9 +266,7 @@ ! Local variables : ! ----------------- - ! Tendencies for the paleoclimate mode - REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run - REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g) - REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi - REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi + ! Variables for the paleoclimate mode + REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi + REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi REAL oblipal ! change of obliquity REAL peri_daypal ! new periday @@ -276,20 +276,21 @@ REAL pdaypal ! new pday = day_ini + step REAL zdt_tot ! time range corresponding to the flux of qsurfyear - REAL massacc(nq) ! accumulated mass - REAL masslost(nq) ! accumulated mass - - REAL globave ! globalaverage 2D ps - REAL globaveice(nq) ! globalaverage 2D ice - REAL globavenewice(nq) ! globalaverage 2D ice - INTEGER lecttsoil ! lecture of tsoil from proftsoil - REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2 - REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1 - REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1 - REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) - REAL vecnull(ngrid,nq) ! null vector used to check conservation of tracer + REAL massacc(nq) ! accumulated mass + REAL masslost(nq) ! accumulated mass + REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2 + REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1 + REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1 + REAL globaveice(nq) ! globalaverage 2D ice + REAL globavenewice(nq) ! globalaverage 2D ice REAL,SAVE :: ptime0 ! store the first time REAL dstep REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier + + ! Other variables + REAL globave ! globalaverage 2D ps + INTEGER lecttsoil ! lecture of tsoil from proftsoil + REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) ! pressure density levels + REAL vecnull(ngrid,nq) ! null vector used to check conservation of tracer ! Aerosol (dust or ice) extinction optical depth at reference wavelength @@ -343,7 +344,4 @@ real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. - REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine - REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine -!$OMP THREADPRIVATE(zdqchim,zdqschim) !! PLUTO variables @@ -443,5 +441,4 @@ real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T). real vmr(ngrid,nlayer) ! volume mixing ratio - real time_phys real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. @@ -667,4 +664,16 @@ endif +! Initialize paleo variables +! ~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (paleo) then + IF (.not. ALLOCATED(qsurfyear)) ALLOCATE(qsurfyear(ngrid,nq)) + IF (.not. ALLOCATED(phisfibed)) ALLOCATE(phisfibed(ngrid)) + write(*,*) 'Paleo time tpal = ',tpal + qsurfyear(:,:)=0. + DO ig=1,ngrid + phisfibed(ig)=phisfi(ig)-qsurf(ig,igcm_n2)*g/rho_q(igcm_n2) ! topo of bedrock below the ice + ENDDO + endif + ! Initialize correlated-k. ! ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -783,21 +792,24 @@ glat(:) = g !AF24: removed oblateness - !zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) - !do l=1,nlayer - ! zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) - !enddo - !zzlev(1:ngrid,1)=0. - !do l=2,nlayer - ! do ig=1,ngrid - ! z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) - ! z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) - ! zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) - ! enddo - !enddo - !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 - !zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) + if (fast) then + zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) + do l=1,nlayer + zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) + enddo + zzlev(1:ngrid,1)=0. + do l=2,nlayer + do ig=1,ngrid + z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) + z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) + zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) + enddo + enddo + !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 + zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) + + else - ! Calculation zzlev & zzlay with g variable (from Mars PCM) - do ig = 1, ngrid + ! Calculation zzlev & zzlay with g variable (from Mars PCM) + do ig = 1, ngrid ! First layer zzlay(ig,1) = -(log(pplay(ig,1)/pplev(ig,1)))*r*pt(ig,1)/g @@ -817,8 +829,11 @@ z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) + print*,ig,l,z1,z2,zzlay(ig,l-1),zzlay(ig,l) zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) enddo zzlev(ig,nlayer+1) = 2*zzlev(ig,nlayer)-zzlev(ig,nlayer-1) - enddo + enddo + + endif !fast ! Compute potential temperature @@ -1981,5 +1996,5 @@ IF (paleo) then call spreadglacier_paleo(ngrid,nq,qsurf, & - phisfinew,dstep,tsurf) + phisfibed,dstep,tsurf) else call spreadglacier_simple(ngrid,nq,qsurf,dstep) @@ -2283,5 +2298,5 @@ ! update new geopotential depending on the ice reservoir - phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000. + phisfipal(:)=phisfibed(:)+qsurfpal(:,igcm_n2)*g/rho_q(igcm_n2) !phisfipal(ig)=phisfi(ig) @@ -2313,18 +2328,10 @@ ! create restartfi if (ngrid.ne.1) then + ztime_restart = ptime + ptimestep/(iphysiq*daysec) call physdem0pal("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & - ptimestep,pdaypal,time_phys,cell_area, & + ptimestep,pdaypal,ztime_restart,cell_area, & albedo_bareground,zmea,zstd,zsig,zgam,zthe, & oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal) - - !call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, & - ! ptimestep,pdaypal, & - ! ztime_restart,tsurf,tsoil,emis,q2,qsurfpal, & - ! cell_area,albedodat,therm_inertia,zmea,zstd,zsig, & - ! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, & - ! peri_daypal) endif - else ! 'paleo' - endif ! end of 'paleo' @@ -2345,5 +2352,5 @@ ! in 'restart'. Between now and the writing of 'restart', ! there will have been the itau=itau+1 instruction and -! a reset of 'time' (lastacll = .true. when itau+1= itaufin) +! a reset of 'time' (lastcall = .true. when itau+1= itaufin) ! thus we store for time=time+dtvr @@ -2364,5 +2371,10 @@ endif ! ngrid endif ! is_omp_master - + if (lastcall.and.paleo) then + call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & + ptimestep,ztime_restart,tsurf, & + tsoil,therm_inertia,emis,albedo,q2,qsurfpal,n2frac) + if (is_master) write(*,*)'PHYSIQ: writing restartfi at time =',ztime_restart + endif @@ -2479,18 +2491,42 @@ ! Diagnostics of optical thickness (dtau = dtau_gas + dtau_rayaer + dtau_cont). ! Warning this is exp(-dtau), I let you postproc with -log to have tau and k itself - ! VI diagnostics for ALICE (/!\ for 28+3 VI bands) - call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,31)) ! 0.185 um - ! IR diagnostics for JWST (/!\ for 20+6 IR bands) - call write_output('dtaui_25250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,3)) ! 25.250 um - call write_output('dtaui_20800nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,5)) ! 20.800 um - call write_output('dtaui_18000nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,7)) ! 18.000 um - call write_output('dtaui_15050nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,9)) ! 15.050 um + !! VI + !call write_output('dtauv_4656nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,2)) ! 4.656 um (28 VIS Bands) + !call write_output('dtauv_1181nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,21)) ! 1.181 um (28 VIS Bands) + !call write_output('dtauv_700nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,24)) ! 0.700 um (28 VIS Bands) + !call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,27)) ! 0.185 um (28 VIS Bands) + !call write_output('dtauv_118nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,28)) ! 0.118 um (28 VIS Bands) + !! IR + !call write_output('dtaui_81250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,2)) ! 81.250 um (17 IR Bands) + !call write_output('dtaui_3859nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,16)) ! 3.859 um (17 IR Bands) + + !call write_output('dtaui_25250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,4)) ! 25.250 um (25 IR Bands) + !call write_output('dtaui_20800nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,6)) ! 20.800 um (25 IR Bands) + !call write_output('dtaui_18000nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,8)) ! 18.000 um (25 IR Bands) + !call write_output('dtaui_15050nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,10)) ! 15.050 um (25 IR Bands) + !if (callmufi) then - ! ! Aerosol optical thickness - ! call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,31,1)) - ! call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,31,2)) - ! ! Aerosols single scattering albedo - ! call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,31,1)) - ! call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,31,2)) + ! Aerosol optical thickness + !call write_output('dtauv_aers_4656nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,1)) + !call write_output('dtauv_aerf_4656nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,2)) + !call write_output('dtauv_aers_1181nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,1)) + !call write_output('dtauv_aerf_1181nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,2)) + !call write_output('dtauv_aers_700nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,1)) + !call write_output('dtauv_aerf_700nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,2)) + !call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,1)) + !call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,2)) + !call write_output('dtauv_aers_118nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,1)) + !call write_output('dtauv_aerf_118nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,2)) + !! Aerosols single scattering albedo + !call write_output('wbarv_aers_4656nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,1)) + !call write_output('wbarv_aerf_4656nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,2)) + !call write_output('wbarv_aers_1181nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,1)) + !call write_output('wbarv_aerf_1181nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,2)) + !call write_output('wbarv_aers_700nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,1)) + !call write_output('wbarv_aerf_700nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,2)) + !call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,1)) + !call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,2)) + !call write_output('wbarv_aers_118nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,1)) + !call write_output('wbarv_aerf_118nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,2)) !endif ! end callmufi Index: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 (revision 4082) +++ trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 (revision 4119) @@ -48,7 +48,8 @@ REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf REAL :: zdqsurf(ngrid) ! tendency of qsurf + REAL :: diff(ngrid) ! height difference (topography+ice layer) between two adjacent points REAL massflow ! function - REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp - INTEGER compt !compteur + REAL hlim,qlim,difflim,stock,H0,totmasstmp + INTEGER compt,nbedge INTEGER slid ! option slid 0 or 1 INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W @@ -63,4 +64,12 @@ difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier zdqsurf(:)=0. + + ! We set the number of edges to make it applicable for 1D case + IF (iim.eq.1) THEN + nbedge=2 ! 1D Case: Only check North/South + ELSE + nbedge=4 ! 3D Cse: Check N,S,E,W + ENDIF + !--------------- Dimensions ------------------------------------- @@ -73,17 +82,15 @@ qlim=hlim*1000. ! kg !! Security for not depleted all ice too fast in one timestep - secu=4 - + secu=2 + !************************************* ! Loop over all points !************************************* DO ig=1,ngrid - !if (ig.eq.ngrid) then - ! print*, 'qpole=',pqsurf(ig,igcm_n2),qlim - !endif !************************************* ! if significant amount of ice, where qsurf > qlim !************************************* + IF (pqsurf(ig,igcm_n2).gt.qlim) THEN @@ -94,6 +101,10 @@ tgla(ig)=ptsurf(ig) ! temperature at the base (we neglect convection) - tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. - if (tbase(ig).gt.63.) then + ! tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. + + ! DEBUG test: paleo simulation without basal melting + tbase(ig) = tgla(ig) + + if (tbase(ig).gt.63.) then slid=1 ! activate slide else @@ -150,4 +161,5 @@ masstmp(:)=0. ! mass to be transferred + diff(:)=0. ! height difference between adjacent points (m) totmasstmp=0. ! total mass to be transferred H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m) @@ -161,7 +173,7 @@ !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point DO i=inddeb,indfin - diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) - if (diff.gt.difflim) then - masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) + diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) + if (diff(i).gt.difflim) then + masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) else masstmp(i)=0. @@ -171,5 +183,5 @@ if (totmasstmp.gt.0.) then !! 2) Available mass to be transferred - stock=maxval(masstmp(:))/secu ! kg/m2/s + stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s !! 3) Real amounts of ice to be transferred : masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb) !kg/m2/s all area around the pole are the same @@ -217,9 +229,13 @@ !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point - DO compt=1,4 + DO compt=1,nbedge i=edge(compt) - diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) - if (diff.gt.difflim) then - masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) + diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) + + + if (diff(i).gt.difflim) then + masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid) + + else masstmp(i)=0. @@ -229,7 +245,7 @@ if (totmasstmp.gt.0) then !! 2) Available mass to be transferred - stock=maxval(masstmp(:))/secu ! kg/m2/s + stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s !! 3) Real amounts of ice to be transferred : - DO compt=1,4 + DO compt=1,nbedge i=edge(compt) masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i) !kg/m2/s @@ -246,5 +262,5 @@ pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 !! add CH4 - DO compt=1,4 + DO compt=1,nbedge i=edge(compt) pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 @@ -257,5 +273,5 @@ pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco !! add CO - DO compt=1,4 + DO compt=1,nbedge i=edge(compt) pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco @@ -265,5 +281,5 @@ ! Add N2 totmasstmp=0. - DO compt=1,4 + DO compt=1,nbedge i=edge(compt) pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep @@ -285,4 +301,5 @@ ENDDO ! ig + end subroutine spreadglacier_paleo