Changeset 4119 for trunk/LMDZ.PLUTO/libf/phypluto
- Timestamp:
- Mar 10, 2026, 10:58:47 PM (4 weeks ago)
- Location:
- trunk/LMDZ.PLUTO/libf/phypluto
- Files:
-
- 2 edited
-
physiq_mod.F90 (modified) (16 diffs)
-
spreadglacier_paleo.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
r4082 r4119 23 23 use aerosol_mod, only: i_haze, haze_prof 24 24 use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & 25 dryness 25 dryness, qsurfyear, phisfibed 26 26 use comdiurn_h, only: coslat, sinlat, coslon, sinlon 27 27 use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen … … 29 29 use geometry_mod, only: latitude, longitude, cell_area, & 30 30 cell_area_for_lonlat_outputs 31 use control_mod, only: iphysiq 31 32 USE comgeomfi_h, only: totarea, totarea_planet 32 33 USE tracer_h, only: noms, mmol, radius, rho_q, qext, & … … 254 255 ! ---------------------- 255 256 integer,save :: day_ini ! Initial date of the run (sol since Ls=0). 257 real,save :: time_phys ! Initial time of day of the run 256 258 integer,save :: icount ! Counter of calls to physiq during the run. 257 259 !$OMP THREADPRIVATE(day_ini,icount) 258 260 259 ! Pluto specific261 ! Pluto specific 260 262 REAL,save :: acond,bcond 261 263 REAL,save :: tcond1p4Pa … … 264 266 ! Local variables : 265 267 ! ----------------- 266 ! Tendencies for the paleoclimate mode 267 REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run 268 REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g) 269 REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi 270 REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi 268 ! Variables for the paleoclimate mode 269 REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi 270 REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi 271 271 REAL oblipal ! change of obliquity 272 272 REAL peri_daypal ! new periday … … 276 276 REAL pdaypal ! new pday = day_ini + step 277 277 REAL zdt_tot ! time range corresponding to the flux of qsurfyear 278 REAL massacc(nq) ! accumulated mass 279 REAL masslost(nq) ! accumulated mass 280 281 REAL globave ! globalaverage 2D ps 282 REAL globaveice(nq) ! globalaverage 2D ice 283 REAL globavenewice(nq) ! globalaverage 2D ice 284 INTEGER lecttsoil ! lecture of tsoil from proftsoil 285 REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2 286 REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1 287 REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1 288 REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) 289 REAL vecnull(ngrid,nq) ! null vector used to check conservation of tracer 278 REAL massacc(nq) ! accumulated mass 279 REAL masslost(nq) ! accumulated mass 280 REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2 281 REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1 282 REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1 283 REAL globaveice(nq) ! globalaverage 2D ice 284 REAL globavenewice(nq) ! globalaverage 2D ice 290 285 291 286 REAL,SAVE :: ptime0 ! store the first time 292 287 REAL dstep 293 288 REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier 289 290 ! Other variables 291 REAL globave ! globalaverage 2D ps 292 INTEGER lecttsoil ! lecture of tsoil from proftsoil 293 REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) ! pressure density levels 294 REAL vecnull(ngrid,nq) ! null vector used to check conservation of tracer 294 295 295 296 ! Aerosol (dust or ice) extinction optical depth at reference wavelength … … 343 344 real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. 344 345 real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. 345 REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine346 REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine347 !$OMP THREADPRIVATE(zdqchim,zdqschim)348 346 349 347 !! PLUTO variables … … 443 441 real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T). 444 442 real vmr(ngrid,nlayer) ! volume mixing ratio 445 real time_phys446 443 447 444 real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. … … 667 664 endif 668 665 666 ! Initialize paleo variables 667 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ 668 if (paleo) then 669 IF (.not. ALLOCATED(qsurfyear)) ALLOCATE(qsurfyear(ngrid,nq)) 670 IF (.not. ALLOCATED(phisfibed)) ALLOCATE(phisfibed(ngrid)) 671 write(*,*) 'Paleo time tpal = ',tpal 672 qsurfyear(:,:)=0. 673 DO ig=1,ngrid 674 phisfibed(ig)=phisfi(ig)-qsurf(ig,igcm_n2)*g/rho_q(igcm_n2) ! topo of bedrock below the ice 675 ENDDO 676 endif 677 669 678 ! Initialize correlated-k. 670 679 ! ~~~~~~~~~~~~~~~~~~~~~~~~ … … 783 792 glat(:) = g !AF24: removed oblateness 784 793 785 !zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) 786 !do l=1,nlayer 787 ! zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) 788 !enddo 789 !zzlev(1:ngrid,1)=0. 790 !do l=2,nlayer 791 ! do ig=1,ngrid 792 ! z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) 793 ! z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) 794 ! zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 795 ! enddo 796 !enddo 797 !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 798 !zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) 794 if (fast) then 795 zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) 796 do l=1,nlayer 797 zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) 798 enddo 799 zzlev(1:ngrid,1)=0. 800 do l=2,nlayer 801 do ig=1,ngrid 802 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) 803 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) 804 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 805 enddo 806 enddo 807 !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 808 zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) 809 810 else 799 811 800 ! Calculation zzlev & zzlay with g variable (from Mars PCM)801 do ig = 1, ngrid812 ! Calculation zzlev & zzlay with g variable (from Mars PCM) 813 do ig = 1, ngrid 802 814 ! First layer 803 815 zzlay(ig,1) = -(log(pplay(ig,1)/pplev(ig,1)))*r*pt(ig,1)/g … … 817 829 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) 818 830 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) 831 print*,ig,l,z1,z2,zzlay(ig,l-1),zzlay(ig,l) 819 832 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 820 833 enddo 821 834 zzlev(ig,nlayer+1) = 2*zzlev(ig,nlayer)-zzlev(ig,nlayer-1) 822 enddo 835 enddo 836 837 endif !fast 823 838 824 839 ! Compute potential temperature … … 1981 1996 IF (paleo) then 1982 1997 call spreadglacier_paleo(ngrid,nq,qsurf, & 1983 phisfi new,dstep,tsurf)1998 phisfibed,dstep,tsurf) 1984 1999 else 1985 2000 call spreadglacier_simple(ngrid,nq,qsurf,dstep) … … 2283 2298 2284 2299 ! update new geopotential depending on the ice reservoir 2285 phisfipal(:)=phisfi new(:)+qsurfpal(:,igcm_n2)*g/1000.2300 phisfipal(:)=phisfibed(:)+qsurfpal(:,igcm_n2)*g/rho_q(igcm_n2) 2286 2301 !phisfipal(ig)=phisfi(ig) 2287 2302 … … 2313 2328 ! create restartfi 2314 2329 if (ngrid.ne.1) then 2330 ztime_restart = ptime + ptimestep/(iphysiq*daysec) 2315 2331 call physdem0pal("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & 2316 ptimestep,pdaypal, time_phys,cell_area, &2332 ptimestep,pdaypal,ztime_restart,cell_area, & 2317 2333 albedo_bareground,zmea,zstd,zsig,zgam,zthe, & 2318 2334 oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal) 2319 2320 !call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &2321 ! ptimestep,pdaypal, &2322 ! ztime_restart,tsurf,tsoil,emis,q2,qsurfpal, &2323 ! cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &2324 ! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, &2325 ! peri_daypal)2326 2335 endif 2327 else ! 'paleo'2328 2329 2336 2330 2337 endif ! end of 'paleo' … … 2345 2352 ! in 'restart'. Between now and the writing of 'restart', 2346 2353 ! there will have been the itau=itau+1 instruction and 2347 ! a reset of 'time' (last acll = .true. when itau+1= itaufin)2354 ! a reset of 'time' (lastcall = .true. when itau+1= itaufin) 2348 2355 ! thus we store for time=time+dtvr 2349 2356 … … 2364 2371 endif ! ngrid 2365 2372 endif ! is_omp_master 2366 2373 if (lastcall.and.paleo) then 2374 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & 2375 ptimestep,ztime_restart,tsurf, & 2376 tsoil,therm_inertia,emis,albedo,q2,qsurfpal,n2frac) 2377 if (is_master) write(*,*)'PHYSIQ: writing restartfi at time =',ztime_restart 2378 endif 2367 2379 2368 2380 … … 2479 2491 ! Diagnostics of optical thickness (dtau = dtau_gas + dtau_rayaer + dtau_cont). 2480 2492 ! Warning this is exp(-dtau), I let you postproc with -log to have tau and k itself 2481 ! VI diagnostics for ALICE (/!\ for 28+3 VI bands) 2482 call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,31)) ! 0.185 um 2483 ! IR diagnostics for JWST (/!\ for 20+6 IR bands) 2484 call write_output('dtaui_25250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,3)) ! 25.250 um 2485 call write_output('dtaui_20800nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,5)) ! 20.800 um 2486 call write_output('dtaui_18000nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,7)) ! 18.000 um 2487 call write_output('dtaui_15050nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,9)) ! 15.050 um 2493 !! VI 2494 !call write_output('dtauv_4656nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,2)) ! 4.656 um (28 VIS Bands) 2495 !call write_output('dtauv_1181nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,21)) ! 1.181 um (28 VIS Bands) 2496 !call write_output('dtauv_700nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,24)) ! 0.700 um (28 VIS Bands) 2497 !call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,27)) ! 0.185 um (28 VIS Bands) 2498 !call write_output('dtauv_118nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,28)) ! 0.118 um (28 VIS Bands) 2499 !! IR 2500 !call write_output('dtaui_81250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,2)) ! 81.250 um (17 IR Bands) 2501 !call write_output('dtaui_3859nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,16)) ! 3.859 um (17 IR Bands) 2502 2503 !call write_output('dtaui_25250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,4)) ! 25.250 um (25 IR Bands) 2504 !call write_output('dtaui_20800nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,6)) ! 20.800 um (25 IR Bands) 2505 !call write_output('dtaui_18000nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,8)) ! 18.000 um (25 IR Bands) 2506 !call write_output('dtaui_15050nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,10)) ! 15.050 um (25 IR Bands) 2507 2488 2508 !if (callmufi) then 2489 ! ! Aerosol optical thickness 2490 ! call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,31,1)) 2491 ! call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,31,2)) 2492 ! ! Aerosols single scattering albedo 2493 ! call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,31,1)) 2494 ! call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,31,2)) 2509 ! Aerosol optical thickness 2510 !call write_output('dtauv_aers_4656nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,1)) 2511 !call write_output('dtauv_aerf_4656nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,2)) 2512 !call write_output('dtauv_aers_1181nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,1)) 2513 !call write_output('dtauv_aerf_1181nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,2)) 2514 !call write_output('dtauv_aers_700nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,1)) 2515 !call write_output('dtauv_aerf_700nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,2)) 2516 !call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,1)) 2517 !call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,2)) 2518 !call write_output('dtauv_aers_118nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,1)) 2519 !call write_output('dtauv_aerf_118nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,2)) 2520 !! Aerosols single scattering albedo 2521 !call write_output('wbarv_aers_4656nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,1)) 2522 !call write_output('wbarv_aerf_4656nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,2)) 2523 !call write_output('wbarv_aers_1181nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,1)) 2524 !call write_output('wbarv_aerf_1181nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,2)) 2525 !call write_output('wbarv_aers_700nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,1)) 2526 !call write_output('wbarv_aerf_700nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,2)) 2527 !call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,1)) 2528 !call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,2)) 2529 !call write_output('wbarv_aers_118nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,1)) 2530 !call write_output('wbarv_aerf_118nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,2)) 2495 2531 !endif ! end callmufi 2496 2532 -
trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90
r4040 r4119 48 48 REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf 49 49 REAL :: zdqsurf(ngrid) ! tendency of qsurf 50 REAL :: diff(ngrid) ! height difference (topography+ice layer) between two adjacent points 50 51 REAL massflow ! function 51 REAL hlim,qlim,difflim, diff,stock,H0,totmasstmp52 INTEGER compt !compteur52 REAL hlim,qlim,difflim,stock,H0,totmasstmp 53 INTEGER compt,nbedge 53 54 INTEGER slid ! option slid 0 or 1 54 55 INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W … … 63 64 difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier 64 65 zdqsurf(:)=0. 66 67 ! We set the number of edges to make it applicable for 1D case 68 IF (iim.eq.1) THEN 69 nbedge=2 ! 1D Case: Only check North/South 70 ELSE 71 nbedge=4 ! 3D Cse: Check N,S,E,W 72 ENDIF 73 65 74 !--------------- Dimensions ------------------------------------- 66 75 … … 73 82 qlim=hlim*1000. ! kg 74 83 !! Security for not depleted all ice too fast in one timestep 75 secu= 476 84 secu=2 85 77 86 !************************************* 78 87 ! Loop over all points 79 88 !************************************* 80 89 DO ig=1,ngrid 81 !if (ig.eq.ngrid) then82 ! print*, 'qpole=',pqsurf(ig,igcm_n2),qlim83 !endif84 90 85 91 !************************************* 86 92 ! if significant amount of ice, where qsurf > qlim 87 93 !************************************* 94 88 95 IF (pqsurf(ig,igcm_n2).gt.qlim) THEN 89 96 … … 94 101 tgla(ig)=ptsurf(ig) 95 102 ! temperature at the base (we neglect convection) 96 tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. 97 if (tbase(ig).gt.63.) then 103 ! tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. 104 105 ! DEBUG test: paleo simulation without basal melting 106 tbase(ig) = tgla(ig) 107 108 if (tbase(ig).gt.63.) then 98 109 slid=1 ! activate slide 99 110 else … … 150 161 151 162 masstmp(:)=0. ! mass to be transferred 163 diff(:)=0. ! height difference between adjacent points (m) 152 164 totmasstmp=0. ! total mass to be transferred 153 165 H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m) … … 161 173 !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point 162 174 DO i=inddeb,indfin 163 diff =H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m)164 if (diff .gt.difflim) then165 masstmp(i)=massflow(ig,i,tgla(ig),diff ,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s)175 diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) 176 if (diff(i).gt.difflim) then 177 masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) 166 178 else 167 179 masstmp(i)=0. … … 171 183 if (totmasstmp.gt.0.) then 172 184 !! 2) Available mass to be transferred 173 stock=m axval(masstmp(:))/secu! kg/m2/s185 stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s 174 186 !! 3) Real amounts of ice to be transferred : 175 187 masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb) !kg/m2/s all area around the pole are the same … … 217 229 218 230 !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point 219 DO compt=1, 4231 DO compt=1,nbedge 220 232 i=edge(compt) 221 diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) 222 if (diff.gt.difflim) then 223 masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) 233 diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) 234 235 236 if (diff(i).gt.difflim) then 237 masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid) 238 239 224 240 else 225 241 masstmp(i)=0. … … 229 245 if (totmasstmp.gt.0) then 230 246 !! 2) Available mass to be transferred 231 stock=m axval(masstmp(:))/secu! kg/m2/s247 stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s 232 248 !! 3) Real amounts of ice to be transferred : 233 DO compt=1, 4249 DO compt=1,nbedge 234 250 i=edge(compt) 235 251 masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i) !kg/m2/s … … 246 262 pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 247 263 !! add CH4 248 DO compt=1, 4264 DO compt=1,nbedge 249 265 i=edge(compt) 250 266 pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 … … 257 273 pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco 258 274 !! add CO 259 DO compt=1, 4275 DO compt=1,nbedge 260 276 i=edge(compt) 261 277 pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco … … 265 281 ! Add N2 266 282 totmasstmp=0. 267 DO compt=1, 4283 DO compt=1,nbedge 268 284 i=edge(compt) 269 285 pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep … … 285 301 286 302 ENDDO ! ig 303 287 304 288 305 end subroutine spreadglacier_paleo
Note: See TracChangeset
for help on using the changeset viewer.
