1
1
PROGRAM allofmoon
2
2
! ***********************************************************************
3
- ! surface temperatures of airless body with H2O exosphere
3
+ ! H2O exosphere and surface temperatures for the Moon
4
4
! ***********************************************************************
5
5
use grid
6
6
use body, only: solarDay, dtsec, Rmoon = > Rbody
@@ -12,19 +12,19 @@ PROGRAM allofmoon
12
12
integer cc(6 ), cc_prod, cc_prod_total, ccc(4 )
13
13
integer :: idum=- 92309 ! random number seed
14
14
real (8 ) tmax, time, tequil
15
- real (8 ), dimension (veclen) :: Tsurf, Qn ! , sigma
15
+ real (8 ), dimension (veclen) :: Tsurf, Qn ! , sigma
16
16
real (8 ) residencetime, HAi
17
17
character (10 ) ext
18
18
19
- real (8 ), dimension (np ,2 ) :: p_r ! longitude(1) and latitude(2)
20
- integer , dimension (np ) :: p_s ! status 0=on surface, 1=inflight, <0= destroyed or trapped
21
- real (8 ), dimension (np ) :: p_t ! time
22
- integer , dimension (np ) :: p_n ! # of hops (diagnostic only)
19
+ real (8 ), dimension (Np ,2 ) :: p_r ! longitude(1) and latitude(2)
20
+ integer , dimension (Np ) :: p_s ! status 0=on surface, 1=inflight, <0= destroyed or trapped
21
+ real (8 ), dimension (Np ) :: p_t ! time
22
+ integer , dimension (Np ) :: p_n ! # of hops (diagnostic only)
23
23
24
24
integer , external :: inbox
25
25
real (8 ), external :: residence_time
26
26
real (8 ), external :: flux_noatm, ran2
27
-
27
+
28
28
! if used with Diviner surface temperature maps as input
29
29
! integer, parameter :: NT=708 ! tied to timestep of 708/(29.53*86400.) = 3596 sec
30
30
! real(8) lon(nlon),lat(nlat),Tbig(nlon,nlat,24)
@@ -46,7 +46,7 @@ PROGRAM allofmoon
46
46
print * ,' Equilibration time' ,tequil,' lunar days'
47
47
print * ,' Maximum time' ,tmax,' lunar days' ,tmax* solarDay/ (86400 .* 365.242 ),' years'
48
48
print * ,' grid size' ,nlon,' x' ,nlat,' veclen=' ,veclen
49
- print * ,' number of molecules' ,np
49
+ print * ,' number of molecules' ,Np
50
50
51
51
HAi = 0 . ! noon
52
52
@@ -57,10 +57,10 @@ PROGRAM allofmoon
57
57
! initial configuration
58
58
! p_t(:)=0.; p_s(:)=0 ! all particles on surface
59
59
p_t(:)= 1d99 ; p_s(:)=- 9 ! no particles
60
- ! do i=1,np
61
- ! p_r(i,1)=360.*ran2(idum);
60
+ ! do i=1,Np
61
+ ! p_r(i,1)=360.*ran2(idum)
62
62
! p_r(i,1)=0.
63
- ! p_r(i,2)=0.;
63
+ ! p_r(i,2)=0.
64
64
! enddo
65
65
66
66
open (unit= 20 ,file= ' Tsurface' ,status= ' unknown' ,action= ' write' )
@@ -75,9 +75,9 @@ PROGRAM allofmoon
75
75
! if Diviner surface temperatures are read in
76
76
! call readTmaps(nlon,nlat,lon,lat,Tbig)
77
77
78
- !- ------loop over time steps
78
+ !- ------loop over time steps
79
79
do n=- nequil,1000000
80
- time = (n+1 )* dtsec ! time at n+1
80
+ time = (n+1 )* dtsec ! time at n+1
81
81
if (time> tmax* solarDay) exit
82
82
83
83
!- - Temperature
@@ -86,7 +86,7 @@ PROGRAM allofmoon
86
86
if (n< 0 ) cycle ! skip remainder
87
87
88
88
! some output
89
- call totalnrs(np ,p_s,cc)
89
+ call totalnrs(Np ,p_s,cc)
90
90
print * ,time/ 3600 .,' One hour call' ,sum (cc(1 :2 ))
91
91
write (30 ,* ) time/ 3600 .,cc(1 :2 ),ccc(1 :4 ),cc_prod_total
92
92
@@ -96,7 +96,7 @@ PROGRAM allofmoon
96
96
cc_prod_total = cc_prod_total + cc_prod
97
97
98
98
! update residence times with new temperature
99
- do i= 1 ,np
99
+ do i= 1 ,Np
100
100
if (p_s(i)==0 ) then ! on surface
101
101
k = inbox(p_r(i,:))
102
102
residencetime = residence_time(Tsurf(k))
@@ -105,14 +105,14 @@ PROGRAM allofmoon
105
105
enddo
106
106
107
107
if (n== 0 ) then ! write out initial distribution
108
- ! call writeparticles(50,np ,p_r,p_s,p_t,p_n)
108
+ ! call writeparticles(50,Np ,p_r,p_s,p_t,p_n)
109
109
call writeglobe(20 ,Tsurf)
110
110
endif
111
111
112
112
write (ext,' (i0.4)' ) n
113
113
call deblank(ext)
114
114
! open(unit=27,file='particles.'//ext,status='unknown',action='write')
115
- ! call writeparticles(27,np ,p_r,p_s,p_t,p_n)
115
+ ! call writeparticles(27,Np ,p_r,p_s,p_t,p_n)
116
116
! close(27)
117
117
! open(unit=28,file='tsurf.'//ext,status='unknown',action='write')
118
118
! call writeglobe(28,Tsurf)
@@ -122,7 +122,7 @@ PROGRAM allofmoon
122
122
! endif
123
123
124
124
! 1 hour of hopping
125
- call montecarlo(np ,idum,p_r,p_s,p_t,p_n,Tsurf,dtsec,ccc,Qn)
125
+ call montecarlo(Np ,idum,p_r,p_s,p_t,p_n,Tsurf,dtsec,ccc,Qn)
126
126
127
127
call totalnrs(Np,p_s,cc)
128
128
@@ -205,8 +205,8 @@ subroutine SurfaceTemperature(dtsec,HAi,time,Tsurf,Qn)
205
205
! toy orbit
206
206
decl = 0 .
207
207
sunR = semia
208
- ! better orbit
209
- ! eps = 1.54*d2r ! lunar obliquity to ecliptic
208
+ ! better orbit
209
+ ! eps = 1.54*d2r ! lunar obliquity to ecliptic
210
210
! ecc=0.; omega=0.
211
211
! call generalorbit(time/86400.,semia,ecc,omega,eps,Ls,decl,sunR)
212
212
@@ -272,11 +272,14 @@ subroutine particles2sigma(Np, p_r, p_s, sigma)
272
272
integer , intent (IN ) :: Np, p_s(Np)
273
273
real (8 ), intent (IN ) :: p_r(Np,2 )
274
274
real (8 ), intent (OUT ) :: sigma(veclen)
275
- integer i, k, nr0(veclen), nr1(veclen), totalnr0, totalnr1
276
- real (8 ) dA(veclen)
275
+ integer i, k, totalnr0, totalnr1
276
+ integer , allocatable :: nr0(:), nr1(:)
277
+ real (8 ), allocatable :: dA(:) ! allocation avoids max-stack-var-size warning
277
278
logical , save :: FirstCall = .TRUE.
278
279
integer , external :: inbox
279
280
281
+ allocate ( nr0(veclen), nr1(veclen), dA(veclen) )
282
+
280
283
nr0(:)= 0 ; nr1(:)= 0
281
284
do i= 1 ,Np
282
285
if (p_s(i)==0 ) then
@@ -293,15 +296,15 @@ subroutine particles2sigma(Np, p_r, p_s, sigma)
293
296
dA = dA* Rbody** 2
294
297
295
298
sigma(:) = nr0(:)/ dA(:)
296
- ! sigma(:) = nr1(:)/dA(:)
299
+ ! sigma(:) = nr1(:)/dA(:)
297
300
! print *,'total area',sum(dA)/(4*pi*Rbody**2) ! test
298
301
totalnr0 = sum (nr0(:))
299
302
totalnr1 = sum (nr1(:))
300
303
! print *,'total # particles in particles2sigma:',Np,totalnr0 ! for checking purposes
301
304
302
305
! where(sigma>maxsigma) maxsigma=sigma
303
306
304
- if (FirstCall) then
307
+ if (FirstCall) then
305
308
open (unit= 40 ,file= ' sigma.dat' ,action= ' write' )
306
309
FirstCall = .FALSE.
307
310
! maxsigma=0.
@@ -310,6 +313,7 @@ subroutine particles2sigma(Np, p_r, p_s, sigma)
310
313
endif
311
314
write (40 ,' (*(1x,g11.5))' ) sigma
312
315
close (40 )
316
+ deallocate (nr0,nr1,dA)
313
317
end subroutine particles2sigma
314
318
315
319
@@ -321,9 +325,11 @@ subroutine fallmap(unit, fall)
321
325
implicit none
322
326
integer , intent (IN ) :: unit, fall(veclen)
323
327
integer i,j,k
324
- real (8 ) sigma(veclen), dA(veclen)
328
+ real (8 ), dimension (:), allocatable :: sigma(:), dA(:)
329
+ ! allocation avoids max-stack-var-size warning
325
330
real (8 ) longitude(nlon), latitude(nlat)
326
331
332
+ allocate ( sigma(veclen), dA(veclen) )
327
333
call areas(dA)
328
334
dA = dA* Rmoon** 2
329
335
@@ -334,13 +340,14 @@ subroutine fallmap(unit, fall)
334
340
write (unit,110 ) 0 .,90 .,sigma(1 )
335
341
do j= 1 ,nlat
336
342
do i= 1 ,nlon
337
- k = 1 + i + (j-1 )* nlon
343
+ k = 1 + i + (j-1 )* nlon
338
344
write (unit,110 ) longitude(i),latitude(j),sigma(k)
339
345
enddo
340
346
enddo
341
347
write (unit,110 ) 0 .,- 90 .,sigma(veclen)
342
348
110 format (f5.1 ,1x ,f6.2 ,1x ,g10.4 )
343
349
344
350
print * ,' Total area' ,sum (dA)
351
+ deallocate ( sigma, dA)
345
352
end subroutine fallmap
346
353
0 commit comments