-
Notifications
You must be signed in to change notification settings - Fork 0
/
init_param_bio.F90
180 lines (144 loc) · 5.71 KB
/
init_param_bio.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
program main
!
! Function: generate ensemble parameters files based on prescribed list of default values.
!
! Usage: ./init_param_bio nens perr styp
! ex) ./init_param_bio 10 10 NORM
!
! nens: number of ensemble members [integer]
! perr: size of perturbation [0-100: integer]
! styp: type of distribution ["NORM": string]
!
! Note: mem000 : default value
! mem001-: perturbed values
!
! Author: [email protected]
!
implicit none
type gen_params
character(len=9) :: pnam
real :: pval, pmin, pmax
end type gen_params
integer(kind=4), parameter :: npar= 8
type(gen_params) :: bio(npar)
integer(kind=4) :: narg ! number of command line arguments
integer(kind=4) :: nens, perr
character(len=10) :: styp
integer(kind=4) :: ip,ie
real :: pstd
real :: rng(1)
real, allocatable :: prm(:,:)
character(len=255) :: filename
character(len=10) :: tmpchar,tmp1,tmp2,tmp3,tmp4
!--- define parameter set
bio(1)%pnam = "grPl" ; bio(1)%pmin = 0.8 ; bio(1)%pmax = 3.0 ! growth rate of Pl [1/day]
bio(2)%pnam = "grPs" ; bio(2)%pmin = 0.5 ; bio(2)%pmax = 2.0 ! growth rate of Ps [1/day]
bio(3)%pnam = "mrPl" ; bio(3)%pmin = 0.02 ; bio(3)%pmax = 0.25 ! mortality rate of Pl [1/day]
bio(4)%pnam = "mrPs" ; bio(4)%pmin = 0.02 ; bio(4)%pmax = 0.25 ! mortality rate of Ps [1/day]
bio(5)%pnam = "mrZl" ; bio(5)%pmin = 0.05 ; bio(5)%pmax = 0.4 ! mortality rate of Zl [1/day]
bio(6)%pnam = "mrZs" ; bio(6)%pmin = 0.05 ; bio(6)%pmax = 0.4 ! mortality rate of Zs [1/day]
bio(7)%pnam = "srDO" ; bio(7)%pmin = 3.5 ; bio(7)%pmax = 20.0 ! sinking rate of Detritus & Opal [m/day]
bio(8)%pnam = "CrSi" ; bio(8)%pmin = 3.0 ; bio(8)%pmax = 13.0 ! Carbon to Silicate ratio [mol C/mol Si]
!--- set default parameter values
bio(1)%pval = 1.30 ! grPl
bio(2)%pval = 1.10 ! grPs
bio(3)%pval = 0.04 ! mrPl
bio(4)%pval = 0.08 ! mrPs
bio(5)%pval = 0.10 ! mrZl
bio(6)%pval = 0.20 ! mrZs
bio(7)%pval = 5.00 ! srDO
bio(8)%pval = 6.625! CrSi
!--- read the command line arguments
narg = command_argument_count() ! number of the command line arguments
if ( narg.lt.3 ) then
write ( *,'(a)' ) "Usage: Ensemble size, Error size [%], Sampling method [NORM/LOGN]"
print *,"Error: reading environment variables"
stop
endif
call get_command_argument(1,tmpchar) ; read(tmpchar,'(i3)') nens
call get_command_argument(2,tmpchar) ; read(tmpchar,'(i3)') perr
call get_command_argument(3,styp)
write ( *,'(a,i3)' ) "Ensemble size : ", nens
write ( *,'(a,i3)' ) "Initial parameter variance [%]: ", perr
write ( *,'(a,a4)' ) "Sampling method [NORM or LOGN]: ", TRIM(styp)
write ( *,'(a)' ) "PARAMETERS:"
do ip = 1,npar
write ( *,'(x,a5)' ) "-"//TRIM(bio(ip)%pnam)
enddo
!--- generate sets of perturbed parameters
allocate( prm(0:nens,npar) )
prm(0,:) = bio(:)%pval
if (nens == 1) then
prm(1,:) = bio(:)%pval
else
call set_random_seed ! set a seed for random number generator
do ie = 1, nens
do ip = 1,npar
call random(rng(1),1)
if ( styp.eq."LOGN" ) then
print *,"Error: sampling method:"//TRIM(styp)//" not supported yet"
stop
!prm(ie,ip) = bio(ip)%pval*exp(0.01*perr*rng(1))
elseif ( styp.eq."NORM" ) then
pstd = 0.01*perr*bio(ip)%pval
prm(ie,ip) = max( 0.0, bio(ip)%pval + pstd*rng(1) )
else
print *,"Error: sampling method:"//TRIM(styp)//" not supported"
stop
endif
if (prm(ie,ip).gt.bio(ip)%pmax) then
prm(ie,ip) = bio(ip)%pmax
elseif (prm(ie,ip).lt.bio(ip)%pmin) then
prm(ie,ip) = bio(ip)%pmin
endif
!print *, rng(1), bio(ip)%pval, bio(ip)%pmin, bio(ip)%pmax, prm(ie,ip)
enddo
enddo
endif
do ie = 0,nens
write( tmp1,'(i0.3)' ) ie
filename="Parameter_bio_mem"//trim(adjustl(tmp1))//".txt"
open( unit=10,file=filename,status='replace',form='formatted' )
do ip = 1,npar
write(10,'(a4,1x,a1,1x,f7.3)') bio(ip)%pnam,"=",prm(ie,ip)
enddo
close( unit=10 )
enddo
! store parameter ensemble in column wise
write(tmp2,'(i3)') npar
tmp3 = '('//trim(adjustl(tmp2))//'f10.5)'
tmp4 = '('//trim(adjustl(tmp2))//'a10)'
open( unit=10,file='Parameter_bio_mem_all.txt',status='replace' )
write( 10,tmp4 ) ( TRIM(bio(ip)%pnam), ip = 1,npar )
do ie = 0,nens
write( 10,tmp3 ) ( prm(ie,ip), ip = 1,npar )
enddo
close( unit=10 )
end program main
subroutine set_random_seed
implicit none
integer :: i, n, clock
integer, dimension(:), allocatable :: seed
call random_seed(size = n)
allocate(seed(n))
call system_clock(count=clock)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
call random_seed(put = seed)
deallocate(seed)
end subroutine set_random_seed
subroutine random(work1,n)
!
! Returns a vector of normally distributed (variance=1,mean=0) random value
! generated with the Box-Muller transformation
!
implicit none
integer, intent(in) :: n
real, intent(out) :: work1(n)
real, allocatable :: work2(:)
real, parameter :: pi=3.141592653589
allocate (work2(n))
call random_number(work1)
call random_number(work2)
work1= sqrt(-2.0*log(work1))*cos(2.0*pi*work2)
deallocate(work2)
end subroutine random