SICOPOLIS V5-dev  Revision 1264
ctrl_m.F90
Go to the documentation of this file.
1 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 
3 ! Module : c t r l
4 
5 ! Purpose : Declarations of control variables for adjointing
6 
7 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 
9 module ctrl_m
10 
11  use sico_types_m
13  use sico_vars_m
14 
15  implicit none
16 
17  public :: ctrl_init
19  public :: read_ad_data
20  public :: cost_final
21 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
22  public :: myceiling, myfloor
23 #endif
24 
25 contains
26 
27 !-------------------------------------------------------------------------------
28 !> Initialiation of control variable.
29 !! Recognized by OpenAD with the prefix xxVar,
30 !! where Var is the variable of choice (normally
31 !! something in sico_variables_m
32 !<------------------------------------------------------------------------------
33  subroutine ctrl_init()
34 
35  implicit none
36 
37  integer(i4b) :: i, j
38 
39  ! CHOICE OF CONTROL VARIABLE IS MADE HERE:
40  !double precision, dimension(0:JMAX,0:IMAX) :: xxH_t
41  !double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxtemp_c
42  !double precision, dimension(0:JMAX,0:IMAX,12) :: xxprecip_present
43  !double precision, dimension(0:JMAX,0:IMAX) :: xxc_slide
44  double precision, dimension(0:JMAX,0:IMAX) :: xxH_c
45 
46  ! 2D controls:
47 
48  do i=0, imax
49  do j=0, jmax
50 
51  ! BC: sliding parameter
52  !xxc_slide(j,i) = 0.0d0
53 
54  ! BC: mean annual precip_presentulation at the ice surface
55  !xxtemp_c(0,j,i) = 0.0d0
56  !xxprecip_present(j,i,7) = 0.0d0
57 
58  ! IC: ice geometry
59  xxh_c(j,i) = 0.0d0
60  !xxH_t(j,i) = 0.0d0
61 
62  end do
63  end do
64 
65  end subroutine ctrl_init
66 
67 !-------------------------------------------------------------------------------
68 !> The independent variable is perturbed by the xxVar
69 !! at all places during the adjoint mode.
70 !<------------------------------------------------------------------------------
71  subroutine cost_independent_init()
72 
73  implicit none
74 
75  integer(i4b) :: i, j
76  !double precision, dimension(0:JMAX,0:IMAX) :: xxH_t
77  !double precision, dimension(0:JMAX,0:IMAX) :: xxc_slide
78  !double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxtemp_c
79  !double precision, dimension(0:JMAX,0:IMAX,12) :: xxprecip_present
80  double precision, dimension(0:JMAX,0:IMAX) :: xxH_c
81 
82  !-------- Initialize independent ctrl variables --------
83  do j=0,jmax
84  do i=0,imax
85 
86  !c_slide(j,i) = c_slide(j,i) + xxc_slide(j,i)
87 
88  ! BC: mean annual precip_presentulation at the ice surface
89  !temp_c(0,j,i) = temp_c(0,j,i) + xxtemp_c(0,j,i)
90  !precip_present(j,i,7) = precip_present(j,i,7) + xxprecip_present(j,i,7)
91 
92  ! IC: ice geometry
93  h_c(j,i) = h_c(j,i) + xxh_c(j,i)
94  !H_t(j,i) = H_t(j,i) + xxH_t(j,i)
95 
96  end do
97  end do
98 
99  end subroutine cost_independent_init
100 
101 !-------------------------------------------------------------------------------
102 !> The dependent variable for the cost routine.
103 !! This must be a scalar variable, although option for
104 !! a summation or weighted average of multiple cost
105 !! targets is possible.
106 !<------------------------------------------------------------------------------
107  subroutine cost_dependent_init()
108 
109  implicit none
110 
111  !------- Final cost value:
112  fc = 0.0_dp
113 
114  !------- Can be sum of multiple cost test values (later):
115  objf_test = 0.0_dp
116  mult_test = 1.0_dp
117 
118  end subroutine cost_dependent_init
119 
120 !-------------------------------------------------------------------------------
121 !> Reading in of data.
122 !<------------------------------------------------------------------------------
123  subroutine read_ad_data()
124 
125  use sico_types_m
126  use sico_variables_m
127  use sico_vars_m
128 
129  implicit none
130 #ifdef AGE_COST
131  integer(i4b) :: ios, i, j, k, kc, kt
132  !real(dp), dimension(0:JMAX,0:IMAX) :: H_data
133 
134  ! If using ice thickness data as a target:
135  !open(unit=90, iostat=ios, file='ice_thick.dat', status='old')
136  !do j=0,JMAX
137  ! read(unit=90, fmt=*) (H_obs(j,i), i=0,IMAX)
138  !end do
139  !close(unit=90)
140 
141  ! Else if using the gridded ages as a target:
142  ! [Note: this file starts with the KT ages first, then the KC]
143  open(unit=90, iostat=ios, file='subroutines/openad/gridded_age.dat', status='old')
144  do k=0,kcmax+ktmax+1
145  do j=0,jmax
146  read(unit=90, fmt=*) (age_data(k,j,i), i=0,imax)
147  end do
148  end do
149  close(unit=90)
150 
151  ! Now we have an option where we can initialize the ages
152  ! (for testing purposes only):
153 
154  do k=0,kcmax+ktmax+1
155  do j=0,jmax
156  do i=0,imax
157  if (k.le.ktmax) then
158  kt=k
159  age_t (kt,j,i) = age_data(k,j,i)*365*24*3600
160  age_t_neu(kt,j,i) = age_data(k,j,i)*365*24*3600
161  else
162  kc=k - (ktmax+1)
163  age_c (kc,j,i) = age_data(k,j,i)*365*24*3600
164  age_c_neu(kc,j,i) = age_data(k,j,i)*365*24*3600
165  end if
166  end do
167  end do
168  end do
169 #endif
170  end subroutine read_ad_data
171 
172 !-------------------------------------------------------------------------------
173 !> This is the final cost calculation.
174 !! The cost function structure is defined here. Currently
175 !! is a "observed age" - modeled age summed over the entire
176 !! domain. The "observed age" is a fake, generated age
177 !! field performed by the 125 ka run in headers.
178 !! Other cost functions (e.g., total cold ice volume,
179 !! commented out below) are certainly possible, and
180 !! recommended!
181 !<------------------------------------------------------------------------------
182  subroutine cost_final(runname)
183 
184  implicit none
185 
186  integer(i4b) :: i, j, k, kc, kt, ios
187  character(len=100), intent(out) :: runname
188 
189  !-------- Calculate the difference between the modeled and 'observed' ages:
190 #ifdef AGE_COST
191  do k=0, ktmax+kcmax+1
192  do i=0, imax
193  do j=0, jmax
194 
195  if (k.lt.ktmax) then
196  kt = k
197  objf_test = objf_test + (age_data(k,j,i)*365*24*3600 - age_t(kt,j,i))
198  else
199  kc = k - (ktmax+1)
200  objf_test = objf_test + (age_data(k,j,i)*365*24*3600 - age_c(kc,j,i))
201  endif
202 
203  end do
204  end do
205  end do
206 #else
207  do i=0, imax
208  do j=0, jmax
209  !--- Other cost functions:
210  objf_test = objf_test + (h_c(j,i) + h_t(j,i))*area(j,i)
211  end do
212  end do
213 #endif
214 
215  fc = mult_test * (objf_test)
216 
217  !-------- Print to screen just in case something gets
218  ! crazy with the file outputting:
219  print *, 'Final cost, fc = ', fc
220  print *, trim(outpath)
221 
222  !-------- Write final cost to a file:
223  open(unit=97, iostat=ios, &
224 #ifndef ALLOW_OPENAD
225  file=trim(outpath)//'/'//trim(runname)//'_COST.dat', &
226 #else
227  file='AD_COST', &
228 #endif
229  status='new')
230  write(unit=97, fmt='(f26.6)') fc
231  write(unit=97, fmt='(2a)') 'Final cost, fc = ', runname
232  close(unit=97)
233 
234  end subroutine cost_final
235 
236 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
237  subroutine myceiling(num, output)
238 
239  use sico_types_m
240  use sico_variables_m
241 
242  implicit none
243 
244  integer(i4b) :: inum
245  real(dp), intent(in) :: num
246  integer(i4b), intent(out) :: output
247 
248  inum = int(num);
249 
250  inum = int(num);
251  if (abs(num-real(inum,dp)) <= &
252  (abs(num+real(inum,dp))*myepsilon_dp) ) then
253  output = inum;
254  else if (num>0) then
255  output = inum + 1;
256  else if (num<0) then
257  output = inum;
258  end if
259  end subroutine myceiling
260 
261  subroutine myfloor(num, output)
262 
263  use sico_types_m
264  use sico_variables_m
265 
266  implicit none
267 
268  integer(i4b) :: inum
269  real(dp), intent(in) :: num
270  integer(i4b), intent(out) :: output
271 
272  inum = int(num);
273 
274  if (abs(num-real(inum,dp)) <= &
275  (abs(num+real(inum,dp))*myepsilon_dp) ) then
276  output = inum;
277  else if (num>0) then
278  output = inum;
279  else if (num<0) then
280  output = inum - 1;
281  end if
282 
283  end subroutine myfloor
284 #endif
285 
286 end module ctrl_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
Definition: ctrl_m.F90:9
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
subroutine, public cost_independent_init()
The independent variable is perturbed by the xxVar at all places during the adjoint mode...
Definition: ctrl_m.F90:72
subroutine, public cost_dependent_init()
The dependent variable for the cost routine. This must be a scalar variable, although option for a su...
Definition: ctrl_m.F90:108
subroutine read_ad_data()
Reading in of data.
subroutine, public ctrl_init()
Initialiation of control variable. Recognized by OpenAD with the prefix xxVar, where Var is the varia...
Definition: ctrl_m.F90:34
subroutine, public cost_final(runname)
This is the final cost calculation. The cost function structure is defined here. Currently is a "obse...
Definition: ctrl_m.F90:183
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.