SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
sico_main_loop.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : s i c o _ m a i n _ l o o p
4 !
5 !> @file
6 !!
7 !! Main loop of SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Main loop of SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 subroutine sico_main_loop(delta_ts, glac_index, &
36  sum_accum, mean_accum, mean_accum_inv, &
37  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38  time, time_init, time_end, time_output, &
39  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40  z_sl, dzsl_dtau, z_mar, &
41  ii, jj, nn, &
42  ndat2d, ndat3d, n_output, &
43  runname)
44 
45 use sico_types
46 use sico_variables
47 
48 implicit none
49 integer(i4b) :: i, j, kc, kt, kr, n
50 integer(i4b) :: itercount, itercount_max
51 integer(i4b) :: iter_temp, iter_wss, iter_ser, iter_out, iter_output(100)
52 integer(i4b) :: ndat2d, ndat3d
53 integer(i4b) :: n_output
54 integer(i4b) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1)), nn(0:jmax,0:imax)
55 real(dp) :: delta_ts, glac_index
56 real(dp) :: sum_accum, mean_accum, mean_accum_inv
57 real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
58 real(dp) :: time, time_init, time_end, time_output(100)
59 real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
60 real(dp) :: dtime_temp_inv
61 real(dp) :: z_sl, dzsl_dtau, z_mar
62 character (len=100) :: runname
63 logical :: flag_3d_output
64 
65 !-------- Begin of main loop --------
66 
67 write(6,*) ' '
68 
69 itercount_max = nint((time_end-time_init)/dtime)
70 
71 iter_temp = nint(dtime_temp/dtime)
72 #if REBOUND==2
73 iter_wss = nint(dtime_wss/dtime)
74 #endif
75 iter_ser = nint(dtime_ser/dtime)
76 #if ( OUTPUT==1 || OUTPUT==3 )
77 iter_out = nint(dtime_out/dtime)
78 #endif
79 #if ( OUTPUT==2 || OUTPUT==3 )
80 do n=1, n_output
81  iter_output(n) = nint((time_output(n)-time_init)/dtime)
82 end do
83 #endif
84 
85 main_loop : do itercount=1, itercount_max
86 
87 write(6,'(i7)') itercount
88 
89 !-------- Update of time --------
90 
91 time = time_init + real(itercount,dp)*dtime
92 
93 !-------- Save old mask --------
94 
95 maske_help = maske
96 
97 !-------- Boundary conditions --------
98 
99 call boundary(time, dtime, dxi, deta, delta_ts, glac_index, &
100  z_sl, dzsl_dtau, z_mar)
101 
102 !-------- Temperature, water content, age, flow enhancement factor --------
103 
104 if ( mod(itercount, iter_temp) == 0 ) then
105 
106  write(6,*) ' -------- Computation of T'
107 
108 ! ------ Temperature, water content, age
109 
110 #if CALCMOD==1
111  call calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
112  dtime_temp, mean_accum_inv) ! polythermal mode
113 
114 #elif CALCMOD==0
115  call calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
116  dtime_temp, mean_accum_inv) ! cold-ice mode
117 #endif
118 
119 ! ------ Time derivative of H_t (further time derivatives are
120 ! computed in subroutine calc_top)
121 
122  dtime_temp_inv = 1.0_dp/dtime_temp
123 
124  dh_t_dtau = (h_t_neu - h_t)*dtime_temp_inv
125 
126 ! ------ New values --> old values
127 
128  n_cts = n_cts_neu
129  zm = zm_neu
130  h_c = h_c_neu
131  h_t = h_t_neu
132  temp_c = temp_c_neu
133  age_c = age_c_neu
134  omega_t= omega_t_neu
135  age_t = age_t_neu
136  temp_r = temp_r_neu
137 
138 ! ------ Flow enhancement factor
139 
140  call calc_enhance(time)
141 
142 end if
143 ! End of computation of temperature, water content, age and
144 ! enhancement factor
145 
146 !-------- Velocity --------
147 
148 call calc_vxy_b_sia(time, z_sl)
149 call calc_vxy_sia(dzeta_c, dzeta_t)
150 
151 #if MARGIN==3
152 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, ii, jj, nn)
153 #endif
154 
155 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
156 
157 #if MARGIN==3
158 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
159 #endif
160 
161 !-------- Topography --------
162 
163 call calc_top(time, time_init, &
164  z_sl, dzsl_dtau, z_mar, &
165  dtime, dxi, deta, dzeta_t, &
166  mean_accum, &
167  ii, jj, nn, itercount, iter_wss)
168 
169 ! ------ New values --> old values
170 
171 zs = zs_neu
172 zm = zm_neu
173 zb = zb_neu
174 zl = zl_neu
175 h_c= h_c_neu
176 h_t= h_t_neu
177 
178 !-------- Melting temperature --------
179 
180 call calc_temp_melt
181 
182 !-------- Basal temperature --------
183 
184 call calc_temp_bas
185 
186 !-------- Basal melting rate --------
187 
188 call calc_qbm(time, z_sl, dzeta_c, dzeta_r)
189 
190 !-------- Effective thickness of subglacial water --------
191 
192 call calc_water_bas
193 
194 !-------- Data output --------
195 
196 #if OUTPUT==1
197 
198 if ( mod(itercount, iter_out) == 0 ) then
199 
200 #if ERGDAT==0
201  flag_3d_output = .false.
202 #elif ERGDAT==1
203  flag_3d_output = .true.
204 #endif
205 
206 #if NETCDF==1
207  call output1(runname, time, delta_ts, glac_index, z_sl, &
208  flag_3d_output, ndat2d, ndat3d)
209 #elif NETCDF==2
210  call output_nc(runname, time, delta_ts, glac_index, z_sl, &
211  flag_3d_output, ndat2d, ndat3d)
212 #endif
213 
214 end if
215 
216 #elif OUTPUT==2
217 
218 do n=1, n_output
219 
220  if (itercount == iter_output(n)) then
221 
222 #if ERGDAT==0
223  flag_3d_output = .false.
224 #elif ERGDAT==1
225  flag_3d_output = .true.
226 #endif
227 
228 #if NETCDF==1
229  call output1(runname, time, delta_ts, glac_index, z_sl, &
230  flag_3d_output, ndat2d, ndat3d)
231 #elif NETCDF==2
232  call output_nc(runname, time, delta_ts, glac_index, z_sl, &
233  flag_3d_output, ndat2d, ndat3d)
234 #endif
235 
236  end if
237 
238 end do
239 
240 #elif OUTPUT==3
241 
242 if ( mod(itercount, iter_out) == 0 ) then
243 
244  flag_3d_output = .false.
245 
246 #if NETCDF==1
247  call output1(runname, time, delta_ts, glac_index, z_sl, &
248  flag_3d_output, ndat2d, ndat3d)
249 #elif NETCDF==2
250  call output_nc(runname, time, delta_ts, glac_index, z_sl, &
251  flag_3d_output, ndat2d, ndat3d)
252 #endif
253 
254 end if
255 
256 do n=1, n_output
257 
258  if (itercount == iter_output(n)) then
259 
260  flag_3d_output = .true.
261 
262 #if NETCDF==1
263  call output1(runname, time, delta_ts, glac_index, z_sl, &
264  flag_3d_output, ndat2d, ndat3d)
265 #elif NETCDF==2
266  call output_nc(runname, time, delta_ts, glac_index, z_sl, &
267  flag_3d_output, ndat2d, ndat3d)
268 #endif
269 
270  end if
271 
272 end do
273 
274 #endif
275 
276 if ( mod(itercount, iter_ser) == 0 ) then
277  call output2(time, dxi, deta, delta_ts, glac_index, z_sl)
278 #if OUTSER==2
279  call output3(time, delta_ts, glac_index, z_sl)
280 #endif
281 #if OUTSER==3
282  call output4(time, dxi, deta, delta_ts, glac_index, z_sl)
283 #endif
284 #if ( defined(ASF) && OUTSER==4 )
285  call output5(time, dxi, deta, delta_ts, glac_index, z_sl)
286 #endif
287 
288 end if
289 
290 end do main_loop ! End of main loop
291 
292 end subroutine sico_main_loop
293 !
Generated on Mon May 20 2013 10:18:43 for SICOPOLIS V3.0 by   doxygen 1.8.1