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