Line data Source code
1 : !> \file mo_neutrons.f90
2 : !> \brief \copybrief mo_neutrons
3 : !> \details \copydetails mo_neutrons
4 :
5 : !> \brief Models to predict neutron intensities above soils
6 : !> \details The number of neutrons above the ground is directly related to
7 : !! the number soil water content in the ground, air, vegetation and/or snow.
8 : !! This module forward-models neutron abundance as a state variable for each cell.
9 : !> \authors Martin Schroen
10 : !> \date Mar 2015
11 : !> \warning THIS MODULE IS WORK IN PROGRESS, DO NOT USE FOR RESEARCH.
12 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
13 : !! mHM is released under the LGPLv3+ license \license_note
14 : !> \ingroup f_mhm
15 : MODULE mo_neutrons
16 :
17 : ! TODO make it faster with pre-calculated horizons and variables
18 : ! TODO use global parameters as linear model
19 :
20 : USE mo_kind, ONLY: i4, dp
21 : IMPLICIT NONE
22 :
23 : ! Neutron forward model using particle transport physics, see Shuttleworth et al. 2013
24 : PUBLIC :: COSMIC
25 :
26 : ! inverse \theta(N) relation based on Desilets et al. 2010
27 : PUBLIC :: DesiletsN0
28 :
29 : ! integration tabular for approximating the neutron flux integral
30 : PUBLIC :: TabularIntegralAFast
31 :
32 : PRIVATE
33 :
34 : CONTAINS
35 :
36 : ! -----------------------------------------------------------------------------------
37 : ! NAME
38 : ! DesiletsN0
39 : !
40 : ! PURPOSE
41 : !> \brief Calculate neutrons from soil moisture for effective soil layer
42 : !> \details Using the N0-relation derived by Desilets, neutron
43 : !> counts above the ground (one value per cell in mHM) can be
44 : !> derived by a semi-empirical, semi-physical relation.
45 : !> The result depends on N0, the neutron counts for 0% soil mositure.
46 : !> This variable is site-specific and is a global parameter in mHM.
47 : ! ------------------------------------------------------------------
48 : ! N0 formula based on Desilets et al. 2010
49 : ! ------------------------------------------------------------------
50 : !
51 : ! CALLING SEQUENCE
52 : ! call DesiletsN0( Moisture(cells,layers), Depths(layers), &
53 : ! N0-parameter, output(cells) )
54 : !
55 : ! INTENT(IN)
56 : !> \param[in] "real(dp), dimension(:) :: SoilMoisture" Soil Moisture
57 : !> \param[in] "real(dp), dimension(:) :: Horizon_depth" Horizon depths
58 : !> \param[in] "real(dp), dimension(:) :: Bd" Bulk density
59 : !> \param[in] "real(dp), dimension(:) :: LatWater" Lattice water
60 : !> \param[in] "real(dp) :: N0" dry neutron counts
61 : !
62 : ! INTENT(INOUT)
63 : ! None
64 : !
65 : ! INTENT(OUT)
66 : !> \param[out] "real(dp), dimension(size(SoilMoisture,1)) :: neutrons" Neutron counts
67 : !
68 : ! INTENT(IN), OPTIONAL
69 : ! None
70 : !
71 : ! INTENT(INOUT), OPTIONAL
72 : ! None
73 : !
74 : ! INTENT(OUT), OPTIONAL
75 : ! None
76 : !
77 : ! RETURN
78 : ! None
79 : !
80 : ! RESTRICTIONS
81 : ! Horizons(1) must not be zero.
82 : !
83 : ! EXAMPLE
84 : ! N0=1500cph, SoilMoisture(1,1)=70mm, Horizons(1)=200mm
85 : ! 1500*(0.372+0.0808/ (70mm/200mm + 0.115))
86 : ! DesiletsN0 = 819cph
87 : !
88 : ! LITERATURE
89 : ! Desilets, D., M. Zreda, and T. P. A. Ferre (2010),
90 : ! Nature's neutron probe: Land surface hydrology at an elusive scale
91 : ! with cosmic rays, WRR, 46, W11505, doi:10.1029/2009WR008726.
92 : !
93 : ! HISTORY
94 : !> \author Martin Schroen
95 : !> \date Mar 2015
96 : !> Modified Rohini Kumar Oct 2021 - Vertical weighting approach
97 : !> for the Neutron count module to mHM - develop branch
98 :
99 17140032 : subroutine DesiletsN0(SoilMoisture, Horizon_depth, Bd, latWater, N0, neutrons)
100 :
101 : use mo_mhm_constants, only: Desilets_a0, Desilets_a1, Desilets_a2
102 : use mo_moment, only: average
103 : implicit none
104 :
105 : real(dp), dimension(:), intent(in) :: SoilMoisture
106 : real(dp), dimension(:), intent(in) :: Horizon_depth
107 : real(dp), dimension(:), intent(in) :: Bd
108 : real(dp), dimension(:), intent(in) :: latWater
109 : real(dp), intent(in) :: N0 ! from global parameters
110 : real(dp), intent(inout) :: neutrons
111 : ! local variables
112 : integer(i4) :: nLayers, LL, nn, nIntervals
113 8570016 : real(dp), dimension(:), allocatable :: Layer_min, Layer_max, Layer_depth
114 8570016 : real(dp) :: average_swc, average_bd
115 8570016 : real(dp) :: D_86_in_cm, D_86_in_mm
116 8570016 : real(dp), dimension(:), allocatable :: cummulative_Layer_weight
117 8570016 : real(dp) :: depth, weight_10mm_spacing, grav_swc
118 :
119 : ! get the # of soil hoizons
120 8570016 : nLayers = size(SoilMoisture)
121 25710048 : allocate( Layer_min(nLayers) )
122 17140032 : allocate( Layer_max(nLayers) )
123 17140032 : allocate( Layer_depth(nLayers) )
124 :
125 : ! assign layer-1
126 8570016 : Layer_min(1) = 0.0_dp
127 8570016 : Layer_max(1) = Horizon_depth(1)
128 : ! rest layers
129 8570016 : do LL = 2, nLayers
130 0 : Layer_min(LL) = Layer_max(LL-1)
131 8570016 : Layer_max(LL) = Horizon_depth(LL)
132 : end do
133 :
134 : ! estimate layer depth [mm]
135 17140032 : Layer_depth(:) = Layer_max(:) - Layer_min(:)
136 :
137 : ! average soil water content (volumetric ones) and Bulk density
138 17140032 : average_swc = average( SoilMoisture(:)/Layer_depth(:) )
139 8570016 : average_bd = average( Bd(:) )
140 :
141 : ! estimate D86 [in cm and mm ]
142 : D_86_in_cm = ( 1.0_dp/average_bd ) * &
143 8570016 : ( 8.321_dp + 0.14249_dp * (0.96655_dp + exp(-0.01_dp)) * (20.0_dp + average_swc)/(0.0429_dp + average_swc) )
144 8570016 : D_86_in_mm = D_86_in_cm * 10.0_dp !# convert cm to mm
145 :
146 :
147 : ! initalise cummulative layer specific weight
148 17140032 : allocate( cummulative_Layer_weight(nLayers) )
149 17140032 : cummulative_Layer_weight(:) = 0.0_dp
150 :
151 : !! only once to calculate weight at equal space at every 10 mm spacing
152 25710048 : nIntervals = nint( maxval(Horizon_depth(:))/10.0_dp )
153 :
154 : ! calculate 10 mm spacing and on fly the cummulative weights
155 8570016 : depth = 0.0_dp
156 179970336 : do nn = 1, nIntervals
157 171400320 : weight_10mm_spacing = exp(-2.0_dp * depth / D_86_in_mm )
158 :
159 : ! estimate the layer specific cummulative weights
160 342800640 : do LL = 1, nLayers
161 342800640 : if( (depth .GE. Layer_min(LL)) .AND. (depth .LT. Layer_max(LL)) ) then
162 171400320 : cummulative_Layer_weight(LL) = cummulative_Layer_weight(LL) + weight_10mm_spacing
163 : end if
164 : end do
165 : ! update depth
166 179970336 : depth = depth + 10.0_dp
167 : end do
168 :
169 :
170 : ! estimate weightage SWC
171 : ! here average_swc is gravemtric swc = vol.swc/Bulk density
172 : average_swc = 0.0_dp
173 17140032 : do LL = 1, nLayers
174 : !! >> add here if organic water is coming !! >> SoilMoisture(LL) + latWater(LL) + organic_water(LL)...
175 8570016 : grav_swc = ( SoilMoisture(LL) + latWater(LL) )/Layer_depth(LL)/Bd(LL)
176 17140032 : average_swc = average_swc + ( grav_swc * cummulative_Layer_weight(LL) )
177 : end do
178 17140032 : average_swc = average_swc / sum( cummulative_Layer_weighT(:) )
179 :
180 : ! calculate neutron count based on depth weighted SM of *D86*
181 8570016 : neutrons = N0 * ( Desilets_a1 + Desilets_a0 / (average_swc + Desilets_a2) )
182 :
183 : !! deallocate variables
184 8570016 : deallocate(Layer_min, Layer_max, Layer_depth, cummulative_Layer_weight)
185 :
186 8570016 : end subroutine DesiletsN0
187 : ! -----------------------------------------------------------------------------------
188 : ! NAME
189 : ! COSMIC
190 : !
191 : ! PURPOSE
192 : !> \brief Calculate neutrons from soil moisture in all layers.
193 : !> \details Neutron counts above the ground (one value per cell in mHM)
194 : !> can be derived by a simplified physical neutron transport simulation.
195 : !> Fast cosmic-Ray neutrons are generated in the soil and attenuated
196 : !> differently in water and soil. The remaining neutrons that reached
197 : !> the surface relate to the profile of soil water content below.
198 : !> Variables like N, alpha and L3 are site-specific and need to be calibrated.
199 : ! ------------------------------------------------------------------
200 : ! COSMIC model based on Shuttleworth et al. 2013
201 : ! ------------------------------------------------------------------
202 : !
203 : ! CALLING SEQUENCE
204 : ! call COSMIC( Moisture(cells,layers), Depths(layers), &
205 : ! COSMIC-parameterset, neutron_integral_AFast, output(cells) )
206 : !
207 : ! INTENT(IN)
208 : !> \param[in] "real(dp), dimension(:) :: SoilMoisture" Soil Moisture
209 : !> \param[in] "real(dp), dimension(:) :: Horizons" Horizon depths
210 : !> \param[in] "real(dp), dimension(:) :: neutron_integral_AFast" Tabular for Int Approx
211 : !> \param[in] "real(dp) :: interc" interception
212 : !> \param[in] "real(dp) :: snowpack" snowpack
213 : !> \param[in] "real(dp) :: L1_No_count" L1_No_count
214 : !> \param[in] "real(dp), dimension(:) :: L1_bulkDens" Bulk Density
215 : !> \param[in] "real(dp), dimension(:) :: L1_latticeWater" Lattice Water
216 : !> \param[in] "real(dp), dimension(:) :: L1_COSMICL3" L3 from the COSMIC module
217 :
218 : !
219 : ! INTENT(INOUT)
220 : !> \param[inout] "real(dp) :: neutrons" Neutron counts
221 : !
222 : ! INTENT(OUT)
223 : ! None
224 : !
225 : ! INTENT(IN), OPTIONAL
226 : ! None
227 : !
228 : ! INTENT(INOUT), OPTIONAL
229 : ! None
230 : !
231 : ! INTENT(OUT), OPTIONAL
232 : ! None
233 : !
234 : ! RETURN
235 : ! None
236 : !
237 : ! RESTRICTIONS
238 : ! Horizons(:) must not be zero.
239 : !
240 : ! EXAMPLE
241 : ! see supplementaries in literature
242 : !
243 : ! LITERATURE
244 : ! J. Shuttleworth, R. Rosolem, M. Zreda, and T. Franz,
245 : ! The COsmic-ray Soil Moisture Interaction Code (COSMIC) for use in data assimilation,
246 : ! HESS, 17, 3205-3217, 2013, doi:10.5194/hess-17-3205-2013
247 : ! Support and Code: http://cosmos.hwr.arizona.edu/Software/cosmic.html
248 : !
249 : ! HISTORY
250 : !> \author Martin Schroen, originally written by Rafael Rosolem
251 : !> \date Mar 2015
252 : !> Rohini Kumar Oct 2021 - Neutron count module to mHM (5.11.2) - develop branch
253 0 : subroutine COSMIC(SoilMoisture, Horizons, neutron_integral_AFast, &
254 0 : interc, snowpack, L1_N0, L1_bulkDens, L1_latticeWater, L1_COSMICL3, &
255 : neutrons )
256 :
257 8570016 : use mo_mhm_constants, only: H2Odens, COSMIC_alpha, COSMIC_L1, COSMIC_L2, COSMIC_L4 !, COSMIC_N
258 : use mo_constants, only: PI_dp
259 : implicit none
260 :
261 : real(dp), dimension(:), intent(in) :: SoilMoisture
262 : real(dp), dimension(:), intent(in) :: Horizons
263 : real(dp), dimension(:), intent(in) :: neutron_integral_AFast
264 : real(dp), intent(in) :: interc
265 : real(dp), intent(in) :: snowpack
266 : real(dp), intent(in) :: L1_N0
267 : real(dp), dimension(:), intent(in) :: L1_bulkDens
268 : real(dp), dimension(:), intent(in) :: L1_latticeWater
269 : real(dp), dimension(:), intent(in) :: L1_COSMICL3
270 : real(dp), intent(inout) :: neutrons
271 :
272 0 : real(dp) :: lambdaHigh
273 0 : real(dp) :: lambdaFast
274 0 : real(dp) :: totflux
275 0 : real(dp) :: sm ! SoilMoisture
276 0 : real(dp) :: lw ! lattice water
277 0 : real(dp) :: bd ! bulk density
278 0 : real(dp) :: L3
279 : integer(i4):: iFlag_snowlayer_intecept ! 1 if snowlayer is active, 0 else
280 :
281 0 : real(dp), dimension(size(Horizons)+1) :: zthick ! Soil layer thickness (cm)
282 0 : real(dp), dimension(:), allocatable :: isoimass ! Integrated dry soil mass above layer (g)
283 0 : real(dp), dimension(:), allocatable :: iwatmass ! Integrated water mass above layer (g)
284 0 : real(dp), dimension(:), allocatable :: hiflux ! High energy neutron flux
285 0 : real(dp), dimension(:), allocatable :: xeff ! Fast neutron source strength of layer
286 0 : real(dp), dimension(:), allocatable :: h2oeffheight! "Effective" height of water in layer (g/cm3)
287 0 : real(dp), dimension(:), allocatable :: h2oeffdens ! "Effective" density of water in layer (g/cm3)
288 0 : real(dp), dimension(:), allocatable :: fastflux ! Contribution to above-ground neutron flux
289 :
290 : integer(i4) :: layers ! Total number of soil layers
291 : integer(i4) :: ll
292 :
293 :
294 : ! # of layers
295 0 : layers = size(SoilMoisture) + 1 ! soil horizons + one additional snowpack and interception layer
296 :
297 : ! allocate and initalize
298 : allocate(hiflux(layers),xeff(layers),&
299 0 : h2oeffdens(layers),h2oeffheight(layers),fastflux(layers),&
300 0 : isoimass(layers),iwatmass(layers))
301 :
302 0 : zthick(:) = 0.0_dp
303 0 : isoimass(:) = 0.0_dp
304 0 : iwatmass(:) = 0.0_dp
305 0 : hiflux(:) = 0.0_dp
306 0 : xeff(:) = 0.0_dp
307 0 : h2oeffdens(:) = 0.0_dp
308 0 : h2oeffheight(:)= 0.0_dp
309 0 : fastflux(:) = 0.0_dp
310 0 : totflux = 0.0_dp
311 0 : lambdaHigh = 0.0_dp
312 0 : lambdaFast = 0.0_dp
313 0 : sm = 0.0_dp
314 0 : lw = 0.0_dp
315 0 : bd = 0.0_dp
316 0 : L3 = 1.0_dp
317 :
318 : ! switch of snowlayer+Intercept ON and off
319 : ! ON = 1; OFF = 0
320 : ! its hardocded here but can be later part of *.nml file
321 0 : iFlag_snowlayer_intecept = 0
322 :
323 : ! layer 1 is the surface layer. layer 2 up to layers are the usual layers
324 0 : do ll = 1,layers
325 :
326 : ! High energy neutron downward flux
327 : ! The integration is now performed at the node of each layer (i.e., center of the layer)
328 :
329 : !ToDo: maybe put zthick into global constants, so it is an input paramter
330 : ! Soil Layers and Thicknesses are constant in mHM, they could be defined outside of this function
331 : ! except the top layer thickness, which is dependend on the snow for example
332 : ! zthick will be in cm, as all heigths are in cm in this module
333 0 : call layerThickness(ll,Horizons,interc,snowpack,zthick)
334 :
335 0 : if( ( zthick(ll) .GT. 0.0_dp ) .AND. ( (iFlag_snowlayer_intecept .GT. 0 ) .OR. (ll.NE.1) ) ) then
336 :
337 : call loopConstants(ll, SoilMoisture(:), L1_bulkDens(:), &
338 0 : L1_latticeWater(:), L1_COSMICL3(:), sm, bd, lw, L3)
339 :
340 0 : if (ll .EQ. 1) then
341 0 : h2oeffdens(ll) = H2Odens/1000.0_dp
342 : else
343 : ! calculate the effective height of water in each layer in cm
344 : ! because neutron standard measurements are in cm
345 0 : call layerWaterHeight(ll,sm,h2oeffheight)
346 : ! divided by the thickness of the layers,we get the effective density
347 0 : h2oeffdens(ll) = (h2oeffheight(ll) + lw/10.0_dp)/zthick(ll)*H2Odens/1000.0_dp
348 : endif
349 :
350 : ! Assuming an area of 1 cm2
351 : ! we integrate the bulkdensity/h2oeffdens down to the middle of the layer ll:
352 0 : isoimass(ll) = bd*(0.5_dp*zthick(ll))*1.0_dp
353 0 : iwatmass(ll) = h2oeffdens(ll)*(0.5_dp*zthick(ll))*1.0_dp
354 0 : if ( ll .gt. 1 ) then
355 0 : isoimass(ll) = isoimass(ll)+isoimass(ll-1)+bd*(0.5_dp*zthick(ll-1))*1.0_dp
356 0 : iwatmass(ll) = iwatmass(ll)+iwatmass(ll-1)+h2oeffdens(ll-1)*(0.5_dp*zthick(ll-1))*1.0_dp
357 : endif
358 :
359 0 : lambdaHigh = isoimass(ll)/COSMIC_L1 + iwatmass(ll)/COSMIC_L2
360 0 : lambdaFast = isoimass(ll)/L3 + iwatmass(ll)/COSMIC_L4
361 :
362 0 : hiflux(ll) = exp(-lambdaHigh)
363 0 : xeff(ll) = zthick(ll)*(COSMIC_alpha*bd + h2oeffdens(ll))
364 :
365 0 : call lookUpIntegral(fastflux(ll),neutron_integral_AFast,lambdaFast)
366 :
367 : ! After contribution from all directions are taken into account,
368 : ! need to multiply fastflux by 2/pi
369 0 : fastflux(ll) = (2.0_dp/PI_dp) * fastflux(ll)
370 :
371 : ! Low energy (fast) neutron upward flux
372 0 : totflux = totflux + hiflux(ll) * xeff(ll) * fastflux(ll)
373 : endif
374 :
375 : enddo
376 :
377 : ! neutrons=COSMIC_N*totflux
378 : !! >> now based on global parameter given in mhm_paramater.nml
379 0 : neutrons = L1_N0*totflux
380 :
381 : !! free space
382 0 : deallocate( hiflux, xeff, h2oeffheight, h2oeffdens, fastflux, isoimass, iwatmass)
383 :
384 0 : end subroutine COSMIC
385 :
386 :
387 :
388 :
389 : ! >>> Loop constants
390 0 : subroutine loopConstants(ll, SoilMoisture,L1_bulkDens,L1_latticeWater,&
391 0 : L1_COSMICL3,sm,bd,lw,L3 )
392 : implicit none
393 : integer(i4), intent(in) :: ll
394 : real(dp), dimension(:), intent(in) :: SoilMoisture
395 : real(dp), dimension(:), intent(in) :: L1_bulkDens
396 : real(dp), dimension(:), intent(in) :: L1_latticeWater
397 : real(dp), dimension(:), intent(in) :: L1_COSMICL3
398 : real(dp), intent(out) :: sm ! SoilMoisture
399 : real(dp), intent(out) :: bd ! Bulk density
400 : real(dp), intent(out) :: lw ! Lattice water
401 : real(dp), intent(out) :: L3 ! L3
402 :
403 0 : if( ll .EQ. 1 ) then
404 : !ToDo
405 0 : sm = 0.0_dp
406 0 : bd = 0.0_dp
407 0 : lw = 0.0_dp
408 0 : L3 = 1.0_dp
409 : else
410 0 : sm = SoilMoisture(ll-1)
411 0 : bd = L1_bulkDens(ll-1)
412 0 : lw = L1_latticeWater(ll-1)
413 0 : L3 = L1_COSMICL3(ll-1)
414 : endif
415 0 : end subroutine loopConstants
416 :
417 :
418 :
419 : ! >> layer thickness
420 0 : subroutine layerThickness(ll,Horizons,interc,snowpack,zthick)
421 : implicit none
422 : integer(i4), intent(in) :: ll
423 : real(dp),dimension(:), intent(in) :: Horizons
424 : real(dp), intent(in) :: interc
425 : real(dp), intent(in) :: snowpack
426 : real(dp),dimension(:), intent(out) :: zthick
427 :
428 0 : if (ll.eq.1) then
429 0 : zthick(ll)=(snowpack+interc)/10.0_dp
430 0 : else if (ll.eq.2) then
431 0 : zthick(ll)=Horizons(ll-1)/10.0_dp
432 : else
433 0 : zthick(ll)=(Horizons(ll-1) - Horizons(ll-2))/10.0_dp
434 : endif
435 0 : end subroutine
436 :
437 :
438 : ! >> layer specific water height
439 0 : subroutine layerWaterHeight(ll,sm,h2oeffheight)
440 : implicit none
441 : integer(i4), intent(in) :: ll
442 : real(dp), intent(in) :: sm
443 : real(dp),dimension(:), intent(out) :: h2oeffheight
444 : ! The effective water height in each layer in each profile:
445 : ! ToDo:This should include in future: roots, soil organic matter
446 0 : h2oeffheight(ll) = sm/10.0_dp
447 0 : end subroutine
448 :
449 : ! integrade a monotonuous function f, dependend on two parameters c and phi
450 : ! xmin and xmax are the borders for the integration
451 : ! if the values for f(xmin) or f(xmax) are undefinde (like exp(-1/0)), they
452 : ! can be set with fxmin, fxmax.
453 : ! eps is for the accuracy of the result. If the function f is monotonuous, the
454 : ! error is at most eps.
455 : ! steps is the maximum number of interpolation points. It is overriding the
456 : ! error and is the maximum number of steps. A specification of the error
457 : ! though still has an impact. If the function is interpolated well enough
458 : ! in a specific flat region regarding the error it can be interpolated better
459 : ! in a less flat region.
460 : !
461 : !For the specific given integral it is very precise with steps=1024
462 0 : subroutine approx_mon_int(res,f,c,xmin,xmax,eps,steps,fxmin,fxmax)
463 : implicit none
464 : real(dp) :: res
465 : real(dp), external :: f
466 : real(dp), intent(in) :: c
467 : real(dp), intent(in) :: xmax
468 : real(dp), intent(in) :: xmin
469 : real(dp), optional :: eps
470 : integer(i4), optional :: steps
471 : real(dp), optional :: fxmin
472 : real(dp), optional :: fxmax
473 :
474 : !locale variables
475 0 : real(dp) :: epstemp
476 : integer(i4) :: stepstemp
477 0 : real(dp) :: fxmintemp
478 0 : real(dp) :: fxmaxtemp
479 :
480 : ! init
481 0 : if (.not. present(eps)) then
482 0 : epstemp=0.001_dp
483 : else
484 0 : epstemp=eps
485 : endif
486 :
487 0 : if (.not. present(steps)) then
488 0 : stepstemp=0
489 : else
490 0 : stepstemp=steps
491 : endif
492 :
493 0 : if (.not. present(fxmin)) then
494 0 : fxmintemp=f(c,xmin)
495 : else
496 0 : fxmintemp=fxmin
497 : endif
498 :
499 0 : if (.not. present(fxmax)) then
500 0 : fxmaxtemp=f(c,xmax)
501 : else
502 0 : fxmaxtemp=fxmax
503 : endif
504 :
505 0 : res=0.0_dp
506 :
507 0 : if (stepstemp .gt. 0) then
508 0 : call approx_mon_int_steps(res,f,c,xmin,xmax,epstemp,stepstemp,fxmintemp,fxmaxtemp)
509 : else
510 0 : call approx_mon_int_eps(res,f,c,xmin,xmax,epstemp,fxmintemp,fxmaxtemp)
511 : endif
512 :
513 0 : end subroutine
514 :
515 0 : recursive subroutine approx_mon_int_steps(res,f,c,xmin,xmax,eps,steps,fxmin,fxmax)
516 : implicit none
517 : real(dp) :: res
518 : real(dp), external :: f
519 : real(dp), intent(in) :: c
520 : real(dp), intent(in) :: xmax
521 : real(dp), intent(in) :: xmin
522 : real(dp), intent(in) :: eps
523 : integer(i4), intent(in) :: steps
524 : real(dp), intent(in) :: fxmin
525 : real(dp), intent(in) :: fxmax
526 :
527 : !locale variables
528 0 : real(dp) :: xm
529 0 : real(dp) :: fxm
530 0 : real(dp) :: err
531 :
532 0 : xm = (xmax+xmin)/2.0_dp
533 0 : fxm= f(c,xm)
534 :
535 0 : err=abs((fxmax-fxm)*(xmax-xm))
536 0 : if ((err .gt. eps).and.(steps .gt. 1)) then
537 0 : call approx_mon_int_steps(res,f,c,xm,xmax,eps/2.0,steps-steps/2,fxm,fxmax)
538 : else
539 0 : res=res+(xmax-xm)*(fxmax+fxm)/2.0_dp
540 : endif
541 :
542 0 : err=abs((fxm-fxmin)*(xm-xmin))
543 0 : if ((err .gt. eps).and.(steps .gt. 1)) then
544 0 : call approx_mon_int_steps(res,f,c,xmin,xm,eps/2.0,steps/2,fxmin,fxm)
545 : else
546 0 : res=res+(xm-xmin)*(fxm+fxmin)/2.0_dp
547 : endif
548 0 : end subroutine
549 :
550 0 : recursive subroutine approx_mon_int_eps(res,f,c,xmin,xmax,eps,fxmin,fxmax)
551 : implicit none
552 : real(dp) :: res
553 : real(dp), external :: f
554 : real(dp), intent(in) :: c
555 : real(dp), intent(in) :: xmax
556 : real(dp), intent(in) :: xmin
557 : real(dp), intent(in) :: eps
558 : real(dp), intent(in) :: fxmin
559 : real(dp), intent(in) :: fxmax
560 :
561 : !locale variables
562 0 : real(dp) :: xm
563 0 : real(dp) :: fxm
564 0 : real(dp) :: err
565 :
566 0 : xm = (xmax+xmin)/2.0_dp
567 0 : fxm= f(c,xm)
568 :
569 0 : err=abs((fxmax-fxm)*(xmax-xm))
570 0 : if (err .gt. eps) then
571 0 : call approx_mon_int_eps(res,f,c,xm,xmax,eps/2.0,fxm,fxmax)
572 : else
573 0 : res=res+(xmax-xm)*(fxmax+fxm)/2.0_dp
574 : endif
575 :
576 0 : err=abs((fxm-fxmin)*(xm-xmin))
577 0 : if (err .gt. eps) then
578 0 : call approx_mon_int_eps(res,f,c,xmin,xm,eps/2.0,fxmin,fxm)
579 : else
580 0 : res=res+(xm-xmin)*(fxm+fxmin)/2.0_dp
581 : endif
582 0 : end subroutine
583 :
584 :
585 : ! -----------------------------------------------------------------------------------
586 : ! NAME
587 : ! TabularIntegralAFast
588 : !
589 : ! PURPOSE
590 : !> \brief Save approximation data for A_fast
591 : !> \details The COSMIC subroutine needs A_fast to be calculated.
592 : !> A_fast=int_{0}^{pi/2} exp(-Lambda_fast(z)/cos(phi) dphi)
593 : !> This subroutine stores data for intsize values for
594 : !> c:=Lambda_fast(z) between 0 and maxC, and will be written
595 : !> into the global array variable neutron_integral_AFast.
596 : !> The calculation of the values is done with a very precise
597 : !> recursive approximation subroutine. That recursive subroutine
598 : !> should not be used inside the time, cells and layer loops, because
599 : !> it is slow.
600 : !> Inside the loops in the module COSMIC the tabular is used to
601 : !> estimate A_fast, if 0<c<maxC, otherwise the recursive
602 : !> approximation is used.
603 : ! ------------------------------------------------------------------
604 : ! TabularIntegralAFast: a tabular for calculations with splines
605 : ! ------------------------------------------------------------------
606 : !
607 : ! CALLING SEQUENCE
608 : ! call TabularIntegralAFast(neutron_integral_AFast,intsize,maxC)
609 : !
610 : ! INTENT(IN)
611 : !> \param[in] "real(dp), dimension(:,:) :: SoilMoisture" Soil Moisture
612 : !> \param[in] "real(dp), dimension(:) :: Horizons" Horizon depths
613 : !> \param[in] "real(dp), dimension(:) :: params" ! N0, N1, N2, alpha0, alpha1, L30, L31
614 : !> \param[in] "integer(i4) :: intsize" ! number of values for the approximation
615 : !> \param[in] "real(dp) :: maxC" ! maximum value for A_fast
616 : !
617 : ! INTENT(INOUT)
618 : !> \param[out] "real(dp), dimension(intsize) :: neutron_integral_AFast" approximation values
619 : !
620 : ! INTENT(OUT)
621 : ! None
622 : !
623 : ! INTENT(IN), OPTIONAL
624 : ! None
625 : !
626 : ! INTENT(INOUT), OPTIONAL
627 : ! None
628 : !
629 : ! INTENT(OUT), OPTIONAL
630 : ! None
631 : !
632 : ! RETURN
633 : ! None
634 : !
635 : ! RESTRICTIONS
636 : ! intsize and maxC must be positive
637 : !
638 : ! EXAMPLE
639 : ! intsize=8000, maxC=20.0_dp
640 : !
641 : ! LITERATURE
642 : ! see splines for example
643 : !
644 : ! HISTORY
645 : !> \author Maren Kaluza
646 : !> \date Nov 2017
647 :
648 0 : subroutine TabularIntegralAFast(integral,maxC)
649 : use mo_constants, only: PI_dp
650 : implicit none
651 : real(dp), dimension(:) :: integral
652 : real(dp), intent(in) :: maxC
653 :
654 : !local variables
655 : integer(i4) :: i
656 0 : real(dp) :: c
657 : integer(i4) :: intsize
658 :
659 0 : intsize=size(integral)-2
660 :
661 0 : do i=1,intsize+1
662 0 : c =real(i-1,dp)*maxC/real(intsize,dp)
663 0 : call approx_mon_int(integral(i),&
664 0 : intgrandFast,c,0.0_dp,PI_dp/2.0_dp,steps=1024,fxmax=0.0_dp)
665 : enddo
666 0 : integral(intsize+2)=maxC
667 0 : end subroutine
668 :
669 : ! if c>1.0, the function can be fitted very nice with gnuplot
670 : ! pi/2*exp(a*x**b)
671 : subroutine TabularIntegralHermAFast(integral,intsize,maxC)
672 0 : use mo_constants, only: PI_dp
673 : implicit none
674 : real(dp), dimension(:) :: integral
675 : integer(i4), intent(in) :: intsize
676 : real(dp), intent(in) :: maxC
677 :
678 : !local variables
679 : integer(i4) :: i
680 : real(dp) :: c
681 :
682 : do i=1,intsize/2+1
683 : c =real(i-1,dp)*maxC/real(intsize/2,dp)
684 : call approx_mon_int(integral(2*i-1),&
685 : intgrandFast,c,0.0_dp,PI_dp/2.0_dp,steps=1024,fxmax=0.0_dp)
686 : call approx_mon_int(integral(2*i),&
687 : intgrandDerivFast,c,0.0_dp,PI_dp/2.0_dp,steps=1024,fxmax=0.0_dp)
688 : integral(2*i)=integral(2*i)*maxC/real(intsize/2,dp)
689 : enddo
690 : end subroutine
691 :
692 : ! if c>1.0, the function can be fitted very nice with gnuplot
693 : ! pi/2*exp(a*x**b)
694 0 : subroutine lookUpIntegral(res,integral,c)
695 : use mo_constants, only: PI_dp
696 : implicit none
697 : real(dp), intent(out):: res
698 : real(dp), dimension(:),intent(in) :: integral
699 : real(dp), intent(in) :: c
700 :
701 : !local variables
702 : integer(i4) :: place
703 0 : real(dp) :: mu
704 : integer(i4) :: intsize
705 0 : real(dp) :: maxC
706 :
707 0 : intsize=size(integral)-2
708 0 : maxC=integral(intsize+2)
709 0 : mu=c*real(intsize,dp)/maxC
710 0 : place=int(mu,i4)+1
711 0 : if (place .gt. intsize) then
712 : !call approx_mon_int(res,intgrandFast,c,0.0_dp,PI_dp/2.0_dp,steps=1024,fxmax=0.0_dp)
713 : !write(*,*) 'Warning: Lambda_Fast is huge. Slow integration used.'
714 0 : res=(PI_dp/2.0_dp)*exp(-1.57406_dp*c**0.815488_dp)
715 : else
716 0 : mu=mu-real(place-1,dp)
717 0 : res=(1.0_dp-mu)*integral(place)+mu*integral(place+1)
718 : res=res
719 : end if
720 0 : end subroutine
721 :
722 : subroutine lookUpHermiteIntegral(res,integral,intsize,c,maxC)
723 0 : use mo_constants, only: PI_dp
724 : implicit none
725 : real(dp) :: res
726 : real(dp), dimension(:),intent(in):: integral
727 : integer(i4), intent(in) :: intsize
728 : real(dp), intent(in) :: c
729 : real(dp), intent(in) :: maxC
730 :
731 : !local variables
732 : integer(i4) :: place
733 : real(dp) :: mu
734 :
735 : mu=c*real(intsize/2,dp)/maxC
736 : place=int(mu,i4)+1
737 : if (place .gt. intsize) then
738 : !call approx_mon_int(res,intgrandFast,c,0.0_dp,PI_dp/2.0_dp,steps=1024,fxmax=0.0_dp)
739 : !write(*,*) 'Warning: Lambda_Fast is huge. Slow integration used.'
740 : res=(PI_dp/2.0_dp)*exp(-1.57406_dp*c**0.815488_dp)
741 : else
742 : mu=mu-real(place-1,dp)
743 : res=h00(mu)*integral(2*place-1)+h01(mu)*integral(2*place+1)+&
744 : h10(mu)*integral(2*place )+h11(mu)*integral(2*place+2)
745 : end if
746 : end subroutine
747 :
748 : !very bad approximation for the integral from the COSMIC paper. Maybe there is
749 : !some copying error? They claim, this function would have an error of less
750 : !then 1/1000 (which can not be true anyway, because the integral goes to zero
751 : !for c->\infty, and the last if case is a polynome with some coefficients
752 : !unequal to zero and therefore tends to \pm \infty), but even in the first 5
753 : !cases this approximation has sometimes an error of about 1/3 in case 4.
754 : subroutine COSMICeffIntegration(res,x)
755 : implicit none
756 : real(dp) :: res
757 : real(dp), intent(in) :: x
758 : !local variables
759 : real(dp) :: a,b,c,d
760 :
761 : ! write(*,*) x
762 : if (x .le. 0.05_dp) then
763 : a=-347.86105_dp
764 : b=41.64233_dp
765 : c=-4.018_dp
766 : d=-0.00018_dp
767 : res=expPolynomDeg3(x,a,b,c,d)
768 : ! write(*,*) 'c1'
769 : else if (x .le. 0.1_dp) then
770 : a=-16.24066_dp
771 : b=6.64468_dp
772 : c=-2.82003_dp
773 : d=-0.01389_dp
774 : res=expPolynomDeg3(x,a,b,c,d)
775 : ! write(*,*) 'c2'
776 : else if (x .le. 0.5_dp) then
777 : a=-0.95245_dp
778 : b=1.44751_dp
779 : c=-2.18933_dp
780 : d=-0.04034_dp
781 : res=expPolynomDeg3(x,a,b,c,d)
782 : ! write(*,*) 'c3'
783 : else if (x .le. 1.0_dp) then
784 : a=-0.09781_dp
785 : b=0.36907_dp
786 : c=-1.72912_dp
787 : d=-0.10761_dp
788 : res=expPolynomDeg3(x,a,b,c,d)
789 : ! write(*,*) 'c4'
790 : else if (x .le. 5.0_dp) then
791 : a=-0.00416_dp
792 : b=0.05808_dp
793 : c=-1.361482_dp
794 : d=-0.25822_dp
795 : res=expPolynomDeg3(x,a,b,c,d)
796 : ! write(*,*) 'c5'
797 : else
798 : a=0.0_dp
799 : b=0.00061_dp
800 : c=-1.04847_dp
801 : d=-0.96617_dp
802 : res=expPolynomDeg3(x,a,b,c,d)
803 : ! write(*,*) 'c6'
804 : endif
805 :
806 : end subroutine
807 :
808 : subroutine oldIntegration(res,c)
809 : use mo_constants, only: PI_dp
810 : implicit none
811 : real(dp) :: res
812 : real(dp), intent(in) :: c
813 :
814 : ! local variables
815 : real(dp) :: zdeg
816 : real(dp) :: zrad
817 : real(dp) :: costheta
818 : real(dp) :: dtheta
819 :
820 : integer(i4) :: angle ! loop indices for an integration interval
821 : ! Angle distribution parameters (HARDWIRED)
822 : ! rr: Using 0.5 deg angle intervals appears to be sufficient
823 : ! rr: (smaller angles increase the computing time for COSMIC)
824 :
825 : dtheta = 0.5_dp*(PI_dp/180.0_dp)
826 :
827 : ! This second loop needs to be done for the distribution of angles for fast neutron release
828 : ! the intent is to loop from 0 to 89.5 by 0.5 degrees - or similar.
829 : ! Because Fortran loop indices are integers, we have to divide the indices by 10 - you get the idea.
830 : res = 0.0_dp
831 : do angle=0,179
832 : zdeg = real(angle,dp)*0.5_dp
833 : zrad = (zdeg*PI_dp)/180.0_dp
834 : costheta = cos(zrad)
835 : ! Angle-dependent low energy (fast) neutron upward flux
836 : res = res + exp(-c/costheta)*dtheta
837 : enddo
838 : end subroutine
839 :
840 0 : function intgrandFast(c,phi)
841 : implicit none
842 : real(dp) :: intgrandFast
843 : real(dp), intent(in) :: c
844 : real(dp), intent(in) :: phi
845 0 : intgrandFast=exp(-c/cos(phi))
846 : return
847 0 : end function
848 :
849 : function intgrandDerivFast(c,phi)
850 : implicit none
851 : real(dp) :: intgrandDerivFast
852 : real(dp), intent(in) :: c
853 : real(dp), intent(in) :: phi
854 : intgrandDerivFast=(-1.0_dp/cos(phi))*exp(-c/cos(phi))
855 : return
856 0 : end function
857 :
858 : function expPolynomDeg3(x,a,b,c,d)
859 : implicit none
860 : real(dp) :: expPolynomDeg3
861 : real(dp), intent(in) :: x
862 : real(dp), intent(in) :: a,b,c,d
863 :
864 : expPolynomDeg3=exp(a*x**3+b*x**2+c*x+d)
865 : return
866 : end function
867 :
868 : ! hermite polynoms
869 : function h00(t)
870 : implicit none
871 : real(dp) :: h00
872 : real(dp), intent(in) :: t
873 : h00=2.0_dp*t**3-3.0_dp*t**2+1.0_dp
874 : return
875 : end function
876 :
877 : function h01(t)
878 : implicit none
879 : real(dp) :: h01
880 : real(dp), intent(in) :: t
881 : h01=-2.0_dp*t**3+3.0_dp*t**2
882 : return
883 : end function
884 :
885 : function h10(t)
886 : implicit none
887 : real(dp) :: h10
888 : real(dp), intent(in) :: t
889 : h10=t**3-2.0_dp*t**2+t
890 : return
891 : end function
892 :
893 : function h11(t)
894 : implicit none
895 : real(dp) :: h11
896 : real(dp), intent(in) :: t
897 : h11=t**3-t**2
898 : return
899 : end function
900 :
901 : END MODULE mo_neutrons
|
