SICOPOLIS V5-dev  Revision 1368
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 :: cost_final
20 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
21  public :: myceiling, myfloor
22 #endif
23 
24 contains
25 
26 !-------------------------------------------------------------------------------
27 !> Initialiation of control variable.
28 !! Recognized by OpenAD with the prefix xxVar,
29 !! where Var is the variable of choice (normally
30 !! something in sico_variables_m
31 !<------------------------------------------------------------------------------
32  subroutine ctrl_init()
33 
34  implicit none
35 
36  integer(i4b) :: i, j, k
37 
38  ! CHOICE OF CONTROL VARIABLE IS MADE HERE:
39  double precision, dimension(0:JMAX,0:IMAX) :: xxH_c
40  double precision, dimension(0:JMAX,0:IMAX) :: xxc_drag
41  double precision, dimension(0:JMAX,0:IMAX) :: xxc_slide
42  double precision, dimension(0:JMAX,0:IMAX) :: xxtemp_ma_present
43  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxtemp_c
44  double precision, dimension(0:JMAX,0:IMAX) :: xxq_geo
45  double precision, dimension(0:JMAX,0:IMAX) :: xxvis_int_g
46  double precision, dimension(0:JMAX,0:IMAX) :: xxdzs_dxi_g
47  double precision, dimension(0:JMAX,0:IMAX) :: xxdzs_deta_g
48  double precision, dimension(0:JMAX,0:IMAX) :: xxcalv_grounded
49  double precision, dimension(0:JMAX,0:IMAX) :: xxQ_bm
50  double precision, dimension(0:JMAX,0:IMAX) :: xxQ_tld
51  double precision, dimension(0:JMAX,0:IMAX) :: xxacc_fact
52  double precision, dimension(0:JMAX,0:IMAX,12) :: xxprecip_present
53  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxsigma_c
54  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxvx_c
55 
56  ! 3D controls:
57  do k=0, kcmax
58  do j=0, jmax
59  do i=0, imax
60  xxtemp_c(k,j,i) = 0.0d0
61  xxvx_c(k,j,i) = 0.0d0
62  xxsigma_c(k,j,i) = 0.0d0
63  end do
64  end do
65  end do
66 
67  ! 2D controls:
68  do i=0, imax
69  do j=0, jmax
70 
71  ! IC: ice geometry
72  xxh_c(j,i) = 0.0d0
73 
74  ! BC: sliding parameter
75  xxc_drag(j,i) = 0.0d0
76  xxc_slide(j,i) = 0.0d0
77 
78  ! BC: temperature-related quantities:
79  xxq_geo(j,i) = 0.0d0
80  xxtemp_ma_present(j,i) = 0.0d0
81  xxvis_int_g(j,i) = 0.0d0
82  xxdzs_dxi_g(j,i) = 0.0d0
83  xxdzs_deta_g(j,i) = 0.0d0
84  xxcalv_grounded(j,i) = 0.0d0
85  xxq_bm(j,i) = 0.0d0
86  xxq_tld(j,i) = 0.0d0
87 
88  ! BC: SMB-related quantities
89  xxprecip_present(j,i,1) = 0.0d0
90  xxacc_fact(j,i) = 0.0d0
91 
92  end do
93  end do
94 
95  end subroutine ctrl_init
96 
97 !-------------------------------------------------------------------------------
98 !> The independent variable is perturbed by the xxVar
99 !! at all places during the adjoint mode.
100 !<------------------------------------------------------------------------------
101  subroutine cost_independent_init()
102 
103  implicit none
104 
105  integer(i4b) :: i, j, k
106  double precision, dimension(0:JMAX,0:IMAX) :: xxH_c
107  double precision, dimension(0:JMAX,0:IMAX) :: xxtemp_ma_present
108  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxtemp_c
109  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxsigma_c
110  double precision, dimension(0:KCMAX,0:JMAX,0:IMAX) :: xxvx_c
111  double precision, dimension(0:JMAX,0:IMAX,12) :: xxprecip_present
112  double precision, dimension(0:JMAX,0:IMAX) :: xxacc_fact
113  double precision, dimension(0:JMAX,0:IMAX) :: xxq_geo
114  double precision, dimension(0:JMAX,0:IMAX) :: xxc_drag
115  double precision, dimension(0:JMAX,0:IMAX) :: xxc_slide
116  double precision, dimension(0:JMAX,0:IMAX) :: xxvis_int_g
117  double precision, dimension(0:JMAX,0:IMAX) :: xxdzs_dxi_g
118  double precision, dimension(0:JMAX,0:IMAX) :: xxdzs_deta_g
119  double precision, dimension(0:JMAX,0:IMAX) :: xxcalv_grounded
120  double precision, dimension(0:JMAX,0:IMAX) :: xxQ_bm
121  double precision, dimension(0:JMAX,0:IMAX) :: xxQ_tld
122 
123  !-------- Initialize independent ctrl variables --------
124 
125  ! 3D controls:
126  do k=0, kcmax
127  do j=0, jmax
128  do i=0, imax
129  temp_c(k,j,i) = temp_c(k,j,i) + xxtemp_c(k,j,i)
130  sigma_c(k,j,i) = sigma_c(k,j,i) + xxsigma_c(k,j,i)
131  vx_c(k,j,i) = vx_c(k,j,i) + xxvx_c(k,j,i)
132  end do
133  end do
134  end do
135 
136  ! 2D controls:
137  do j=0,jmax
138  do i=0,imax
139 
140  ! IC: ice geometry
141  h_c(j,i) = h_c(j,i) + xxh_c(j,i)
142 
143  ! Basal processes:
144  c_drag(j,i) = c_drag(j,i) + xxc_drag(j,i)
145  c_slide(j,i) = c_slide(j,i) + xxc_slide(j,i)
146 
147  ! BC: temperature-related quantities
148  q_geo(j,i) = q_geo(j,i) + xxq_geo(j,i)
149  temp_ma_present(j,i) = temp_ma_present(j,i) + xxtemp_ma_present(j,i)
150  vis_int_g(j,i) = vis_int_g(j,i) + xxvis_int_g(j,i)
151  dzs_dxi_g(j,i) = dzs_dxi_g(j,i) + xxdzs_dxi_g(j,i)
152  dzs_deta_g(j,i) = dzs_deta_g(j,i) + xxdzs_deta_g(j,i)
153  calv_grounded(j,i) = calv_grounded(j,i) + xxcalv_grounded(j,i)
154  q_bm(j,i) = q_bm(j,i) + xxq_bm(j,i)
155  q_tld(j,i) = q_tld(j,i) + xxq_tld(j,i)
156 
157  ! BC: SMB-related quantities
158  precip_present(j,i,1) = precip_present(j,i,1) + xxprecip_present(j,i,1)
159  acc_fact(j,i) = acc_fact(j,i) + xxacc_fact(j,i)
160 
161  end do
162  end do
163 
164  end subroutine cost_independent_init
165 
166 !-------------------------------------------------------------------------------
167 !> The dependent variable for the cost routine.
168 !! This must be a scalar variable, although option for
169 !! a summation or weighted average of multiple cost
170 !! targets is possible.
171 !<------------------------------------------------------------------------------
172  subroutine cost_dependent_init()
173 
174  implicit none
175 
176  !------- Final cost value:
177  fc = 0.0_dp
178 
179  !------- Can be sum of multiple cost test values (later):
180  objf_test = 0.0_dp
181  mult_test = 1.0_dp
182 
183  end subroutine cost_dependent_init
184 
185 !-------------------------------------------------------------------------------
186 !> This is the final cost calculation.
187 !! The cost function structure is defined here. Currently
188 !! is a "observed age" - modeled age summed over the entire
189 !! domain. The "observed age" is a fake, generated age
190 !! field performed by the 125 ka run in headers.
191 !! Other cost functions (e.g., total cold ice volume,
192 !! commented out below) are certainly possible, and
193 !! recommended!
194 !<------------------------------------------------------------------------------
195  subroutine cost_final(runname)
196 
197  implicit none
198 
199  integer(i4b) :: i, j, k, kc, kt, ios, KDATA
200  character(len=100), intent(out) :: runname
201 
202  !-------- Calculate the difference between the modeled and 'observed' ages:
203 #ifdef AGE_COST
204 
205 #if (CALCMOD!=1)
206  kdata = kcmax
207 #else
208  kdata = kcmax + ktmax
209 print *, '>>> error: CALCMOD == 1 but final cost not properly working for '
210 print *, ' AGE_COST simulations'
211 #endif
212 
213  do k=0, kdata
214  do i=0, imax
215  do j=0, jmax
216 
217  ! only counting points that are real in the data:
218  if ( age_data(k,j,i).lt. 0.0 &
219  .and.age_unc(k,j,i) .lt. 0.0 ) then
220 
221  objf_test = objf_test + &
222  sqrt(((age_data(k,j,i) - age_c(k,j,i)))**2 &
223  / ((age_unc(k,j,i))**2))
224 
225  end if
226  end do
227  end do
228  end do
229 
230 #else
231 
232  do i=0, imax
233  do j=0, jmax
234  !--- Other cost functions:
235  objf_test = objf_test + (h_c(j,i) + h_t(j,i))*area(j,i)
236  end do
237  end do
238 
239 #endif
240 
241  fc = mult_test * (objf_test)
242 
243  !-------- Print to screen just in case something gets
244  ! crazy with the file outputting:
245  print *, 'Final cost, fc = ', fc
246  print *, trim(outpath)
247 
248  !-------- Write final cost to a file:
249  open(unit=97, iostat=ios, &
250 #ifndef ALLOW_OPENAD
251  file=trim(outpath)//'/'//trim(runname)//'_COST.dat', &
252 #else
253  file='AD_COST', &
254 #endif
255  status='new')
256  write(unit=97, fmt='(f26.6)') fc
257  write(unit=97, fmt='(2a)') 'Final cost, fc = ', runname
258  close(unit=97)
259 
260  end subroutine cost_final
261 
262 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
263  subroutine myceiling(num, output)
264 
265  implicit none
266 
267  integer(i4b) :: inum
268  real(dp), intent(in) :: num
269  integer(i4b), intent(out) :: output
270 
271  inum = int(num);
272 
273  inum = int(num);
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 + 1;
279  else if (num<0) then
280  output = inum;
281  end if
282  end subroutine myceiling
283 
284  subroutine myfloor(num, output)
285 
286  implicit none
287 
288  integer(i4b) :: inum
289  real(dp), intent(in) :: num
290  integer(i4b), intent(out) :: output
291 
292  inum = int(num);
293 
294  if (abs(num-real(inum,dp)) <= &
295  (abs(num+real(inum,dp))*myepsilon_dp) ) then
296  output = inum;
297  else if (num>0) then
298  output = inum;
299  else if (num<0) then
300  output = inum - 1;
301  end if
302 
303  end subroutine myfloor
304 #endif
305 
306 end module ctrl_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
Definition: ctrl_m.F90:9
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SSA, at (i,j)
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:102
real(dp), dimension(0:jmax, 0:imax) temp_ma_present
temp_ma_present(j,i): Present-day mean annual surface temperature
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:173
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine, public ctrl_init()
Initialiation of control variable. Recognized by OpenAD with the prefix xxVar, where Var is the varia...
Definition: ctrl_m.F90:33
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
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:196
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
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.