SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
topography3_nc.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : t o p o g r a p h y 3 _ n c
4 !
5 !> @file
6 !!
7 !! Definition of the initial surface and bedrock topography
8 !! (including gradients) and of the horizontal grid spacings dxi, deta.
9 !! For initial topography from previous simulation (read in NetCDF format).
10 !!
11 !! @section Copyright
12 !!
13 !! Copyright 2009-2013 Ralf Greve
14 !!
15 !! @section License
16 !!
17 !! This file is part of SICOPOLIS.
18 !!
19 !! SICOPOLIS is free software: you can redistribute it and/or modify
20 !! it under the terms of the GNU General Public License as published by
21 !! the Free Software Foundation, either version 3 of the License, or
22 !! (at your option) any later version.
23 !!
24 !! SICOPOLIS is distributed in the hope that it will be useful,
25 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
26 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 !! GNU General Public License for more details.
28 !!
29 !! You should have received a copy of the GNU General Public License
30 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !<
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 
34 !-------------------------------------------------------------------------------
35 !> Definition of the initial surface and bedrock topography
36 !! (including gradients) and of the horizontal grid spacings dxi, deta.
37 !! For initial topography from previous simulation (read in NetCDF format).
38 !<------------------------------------------------------------------------------
39 subroutine topography3_nc(dxi, deta, z_sl, anfdatname)
40 
41 use sico_types
42 use sico_variables
43 use netcdf
44 
45 implicit none
46 
47 character (len=100), intent (in) :: anfdatname
48 
49 real(dp), intent (out) :: dxi, deta, z_sl
50 
51 integer(i4b) :: i, j, kc, kt, kr
52 integer(i4b) :: ios
53 real(sp) :: time_conv, dummy_conv, z_sl_conv, &
54  h_r_conv, xi_conv(0:imax), eta_conv(0:jmax), &
55  sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
56  sigma_level_r_conv(0:krmax), &
57  lambda_conv(0:jmax,0:imax), phi_conv(0:jmax,0:imax), &
58  temp_s_conv(0:jmax,0:imax), as_perp_conv(0:jmax,0:imax), &
59  zs_conv(0:jmax,0:imax), zm_conv(0:jmax,0:imax), &
60  zb_conv(0:jmax,0:imax), zl_conv(0:jmax,0:imax), &
61  h_c_conv(0:jmax,0:imax), h_t_conv(0:jmax,0:imax), &
62  h_conv(0:jmax,0:imax), &
63  q_bm_conv(0:jmax,0:imax), q_tld_conv(0:jmax,0:imax), &
64  am_perp_conv(0:jmax,0:imax), &
65  qx_conv(0:jmax,0:imax), qy_conv(0:jmax,0:imax), &
66  temp_r_conv(0:krmax,0:jmax,0:imax), &
67  vx_c_conv(0:kcmax,0:jmax,0:imax), &
68  vy_c_conv(0:kcmax,0:jmax,0:imax), &
69  vz_c_conv(0:kcmax,0:jmax,0:imax), &
70  vx_t_conv(0:ktmax,0:jmax,0:imax), &
71  vy_t_conv(0:ktmax,0:jmax,0:imax), &
72  vz_t_conv(0:ktmax,0:jmax,0:imax), &
73  temp_c_conv(0:kcmax,0:jmax,0:imax), &
74  age_c_conv(0:kcmax,0:jmax,0:imax), &
75  omega_t_conv(0:ktmax,0:jmax,0:imax), &
76  age_t_conv(0:ktmax,0:jmax,0:imax), &
77  dzs_dtau_conv(0:jmax,0:imax), dzm_dtau_conv(0:jmax,0:imax), &
78  dzb_dtau_conv(0:jmax,0:imax), dzl_dtau_conv(0:jmax,0:imax), &
79  dh_c_dtau_conv(0:jmax,0:imax), dh_t_dtau_conv(0:jmax,0:imax), &
80  dh_dtau_conv(0:jmax,0:imax), &
81  vx_b_g_conv(0:jmax,0:imax), vy_b_g_conv(0:jmax,0:imax), &
82  vz_b_conv(0:jmax,0:imax), &
83  vx_s_g_conv(0:jmax,0:imax), vy_s_g_conv(0:jmax,0:imax), &
84  vz_s_conv(0:jmax,0:imax), &
85  temp_b_conv(0:jmax,0:imax), temph_b_conv(0:jmax,0:imax), &
86  p_b_w_conv(0:jmax,0:imax), h_w_conv(0:jmax,0:imax), &
87  q_gl_g_conv(0:jmax,0:imax)
88 
89 integer(i4b) :: ncid, ncv
90 ! ncid: ID of the output file
91 ! ncv: Variable ID
92 
93 !-------- Read data from time-slice file of previous simulation --------
94 
95 call check( nf90_open(anfdatpath//'/'//trim(anfdatname), nf90_nowrite, ncid) )
96 
97 call check( nf90_inq_varid(ncid, 'time', ncv) )
98 call check( nf90_get_var(ncid, ncv, time_conv) )
99 
100 if ( nf90_inq_varid(ncid, 'delta_ts', ncv) == nf90_noerr ) then
101  call check( nf90_get_var(ncid, ncv, dummy_conv) )
102 else if ( nf90_inq_varid(ncid, 'glac_index', ncv) == nf90_noerr ) then
103  call check( nf90_get_var(ncid, ncv, dummy_conv) )
104 else
105  continue
106 end if
107 
108 call check( nf90_inq_varid(ncid, 'z_sl', ncv) )
109 call check( nf90_get_var(ncid, ncv, z_sl_conv) )
110 
111 call check( nf90_inq_varid(ncid, 'xi', ncv) )
112 call check( nf90_get_var(ncid, ncv, xi_conv) )
113 
114 call check( nf90_inq_varid(ncid, 'eta', ncv) )
115 call check( nf90_get_var(ncid, ncv, eta_conv) )
116 
117 call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv) )
118 call check( nf90_get_var(ncid, ncv, sigma_level_c_conv) )
119 
120 call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv) )
121 call check( nf90_get_var(ncid, ncv, sigma_level_t_conv) )
122 
123 call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv) )
124 call check( nf90_get_var(ncid, ncv, sigma_level_r_conv) )
125 
126 call check( nf90_inq_varid(ncid, 'lambda', ncv) )
127 call check( nf90_get_var(ncid, ncv, lambda_conv) )
128 
129 call check( nf90_inq_varid(ncid, 'phi', ncv) )
130 call check( nf90_get_var(ncid, ncv, phi_conv) )
131 
132 call check( nf90_inq_varid(ncid, 'temp_s', ncv) )
133 call check( nf90_get_var(ncid, ncv, temp_s_conv) )
134 
135 call check( nf90_inq_varid(ncid, 'as_perp', ncv) )
136 call check( nf90_get_var(ncid, ncv, as_perp_conv) )
137 
138 call check( nf90_inq_varid(ncid, 'maske', ncv) )
139 call check( nf90_get_var(ncid, ncv, maske) )
140 
141 call check( nf90_inq_varid(ncid, 'n_cts', ncv) )
142 call check( nf90_get_var(ncid, ncv, n_cts) )
143 
144 call check( nf90_inq_varid(ncid, 'zs', ncv) )
145 call check( nf90_get_var(ncid, ncv, zs_conv) )
146 
147 call check( nf90_inq_varid(ncid, 'zm', ncv) )
148 call check( nf90_get_var(ncid, ncv, zm_conv) )
149 
150 call check( nf90_inq_varid(ncid, 'zb', ncv) )
151 call check( nf90_get_var(ncid, ncv, zb_conv) )
152 
153 call check( nf90_inq_varid(ncid, 'zl', ncv) )
154 call check( nf90_get_var(ncid, ncv, zl_conv) )
155 
156 call check( nf90_inq_varid(ncid, 'H_c', ncv) )
157 call check( nf90_get_var(ncid, ncv, h_c_conv) )
158 
159 call check( nf90_inq_varid(ncid, 'H_t', ncv) )
160 call check( nf90_get_var(ncid, ncv, h_t_conv) )
161 
162 call check( nf90_inq_varid(ncid, 'H', ncv) )
163 call check( nf90_get_var(ncid, ncv, h_conv) )
164 
165 call check( nf90_inq_varid(ncid, 'H_R', ncv) )
166 call check( nf90_get_var(ncid, ncv, h_r_conv) )
167 
168 call check( nf90_inq_varid(ncid, 'vx_c', ncv) )
169 call check( nf90_get_var(ncid, ncv, vx_c_conv) )
170 
171 call check( nf90_inq_varid(ncid, 'vy_c', ncv) )
172 call check( nf90_get_var(ncid, ncv, vy_c_conv) )
173 
174 call check( nf90_inq_varid(ncid, 'vz_c', ncv) )
175 call check( nf90_get_var(ncid, ncv, vz_c_conv) )
176 
177 call check( nf90_inq_varid(ncid, 'vx_t', ncv) )
178 call check( nf90_get_var(ncid, ncv, vx_t_conv) )
179 
180 call check( nf90_inq_varid(ncid, 'vy_t', ncv) )
181 call check( nf90_get_var(ncid, ncv, vy_t_conv) )
182 
183 call check( nf90_inq_varid(ncid, 'vz_t', ncv) )
184 call check( nf90_get_var(ncid, ncv, vz_t_conv) )
185 
186 call check( nf90_inq_varid(ncid, 'temp_c', ncv) )
187 call check( nf90_get_var(ncid, ncv, temp_c_conv) )
188 
189 call check( nf90_inq_varid(ncid, 'omega_t', ncv) )
190 call check( nf90_get_var(ncid, ncv, omega_t_conv) )
191 
192 call check( nf90_inq_varid(ncid, 'temp_r', ncv) )
193 call check( nf90_get_var(ncid, ncv, temp_r_conv) )
194 
195 call check( nf90_inq_varid(ncid, 'Q_bm', ncv) )
196 call check( nf90_get_var(ncid, ncv, q_bm_conv) )
197 
198 call check( nf90_inq_varid(ncid, 'Q_tld', ncv) )
199 call check( nf90_get_var(ncid, ncv, q_tld_conv) )
200 
201 call check( nf90_inq_varid(ncid, 'am_perp', ncv) )
202 call check( nf90_get_var(ncid, ncv, am_perp_conv) )
203 
204 call check( nf90_inq_varid(ncid, 'qx', ncv) )
205 call check( nf90_get_var(ncid, ncv, qx_conv) )
206 
207 call check( nf90_inq_varid(ncid, 'qy', ncv) )
208 call check( nf90_get_var(ncid, ncv, qy_conv) )
209 
210 call check( nf90_inq_varid(ncid, 'age_c', ncv) )
211 call check( nf90_get_var(ncid, ncv, age_c_conv) )
212 
213 call check( nf90_inq_varid(ncid, 'age_t', ncv) )
214 call check( nf90_get_var(ncid, ncv, age_t_conv) )
215 
216 call check( nf90_inq_varid(ncid, 'dzs_dtau', ncv) )
217 call check( nf90_get_var(ncid, ncv, dzs_dtau_conv) )
218 
219 call check( nf90_inq_varid(ncid, 'dzm_dtau', ncv) )
220 call check( nf90_get_var(ncid, ncv, dzm_dtau_conv) )
221 
222 call check( nf90_inq_varid(ncid, 'dzb_dtau', ncv) )
223 call check( nf90_get_var(ncid, ncv, dzb_dtau_conv) )
224 
225 call check( nf90_inq_varid(ncid, 'dzl_dtau', ncv) )
226 call check( nf90_get_var(ncid, ncv, dzl_dtau_conv) )
227 
228 call check( nf90_inq_varid(ncid, 'dH_c_dtau', ncv) )
229 call check( nf90_get_var(ncid, ncv, dh_c_dtau_conv) )
230 
231 call check( nf90_inq_varid(ncid, 'dH_t_dtau', ncv) )
232 call check( nf90_get_var(ncid, ncv, dh_t_dtau_conv) )
233 
234 call check( nf90_inq_varid(ncid, 'dH_dtau', ncv) )
235 call check( nf90_get_var(ncid, ncv, dh_dtau_conv) )
236 
237 if ( nf90_inq_varid(ncid, 'vx_b_g', ncv) == nf90_noerr ) then
238  call check( nf90_get_var(ncid, ncv, vx_b_g_conv) )
239 else
240  go to 110
241 end if
242 
243 if ( nf90_inq_varid(ncid, 'vy_b_g', ncv) == nf90_noerr ) then
244  call check( nf90_get_var(ncid, ncv, vy_b_g_conv) )
245 else
246  go to 110
247 end if
248 
249 if ( nf90_inq_varid(ncid, 'vz_b', ncv) == nf90_noerr ) then
250  call check( nf90_get_var(ncid, ncv, vz_b_conv) )
251 else
252  go to 110
253 end if
254 
255 if ( nf90_inq_varid(ncid, 'vx_s_g', ncv) == nf90_noerr ) then
256  call check( nf90_get_var(ncid, ncv, vx_s_g_conv) )
257 else
258  go to 110
259 end if
260 
261 if ( nf90_inq_varid(ncid, 'vy_s_g', ncv) == nf90_noerr ) then
262  call check( nf90_get_var(ncid, ncv, vy_s_g_conv) )
263 else
264  go to 110
265 end if
266 
267 if ( nf90_inq_varid(ncid, 'vz_s', ncv) == nf90_noerr ) then
268  call check( nf90_get_var(ncid, ncv, vz_s_conv) )
269 else
270  go to 110
271 end if
272 
273 if ( nf90_inq_varid(ncid, 'temp_b', ncv) == nf90_noerr ) then
274  call check( nf90_get_var(ncid, ncv, temp_b_conv) )
275 else
276  go to 110
277 end if
278 
279 if ( nf90_inq_varid(ncid, 'temph_b', ncv) == nf90_noerr ) then
280  call check( nf90_get_var(ncid, ncv, temph_b_conv) )
281 else
282  go to 110
283 end if
284 
285 if ( nf90_inq_varid(ncid, 'p_b_w', ncv) == nf90_noerr ) then
286  call check( nf90_get_var(ncid, ncv, p_b_w_conv) )
287 else
288  go to 110
289 end if
290 
291 if ( nf90_inq_varid(ncid, 'H_w', ncv) == nf90_noerr ) then
292  call check( nf90_get_var(ncid, ncv, h_w_conv) )
293 else
294  go to 110
295 end if
296 
297 if ( nf90_inq_varid(ncid, 'q_gl_g', ncv) == nf90_noerr ) then
298  call check( nf90_get_var(ncid, ncv, q_gl_g_conv) )
299 else
300  go to 110
301 end if
302 
303 go to 120
304 
305 110 continue
306 write(6,'(/1x,a)') 'topography3_nc: End-of-file condition while reading *.nc.'
307 write(6, '(1x,a)') ' Some variables will be undefined.'
308 
309 120 continue
310 
311 call check( nf90_close(ncid) )
312 
313 !-------- Convert data to double precision --------
314 
315 z_sl = real(z_sl_conv,dp)
316 
317 h_r = real(h_r_conv,dp)
318 
319 do i=0, imax
320 do j=0, jmax
321 
322  xi(i) = real(xi_conv(i),dp)
323  eta(j) = real(eta_conv(j),dp)
324  zs(j,i) = real(zs_conv(j,i),dp)
325  zm(j,i) = real(zm_conv(j,i),dp)
326  zb(j,i) = real(zb_conv(j,i),dp)
327  zl(j,i) = real(zl_conv(j,i),dp)
328  h_c(j,i) = real(H_c_conv(j,i),dp)
329  h_t(j,i) = real(H_t_conv(j,i),dp)
330  q_bm(j,i) = real(Q_bm_conv(j,i),dp)
331  q_tld(j,i) = real(Q_tld_conv(j,i),dp)
332  am_perp(j,i) = real(am_perp_conv(j,i),dp)
333  qx(j,i) = real(qx_conv(j,i),dp)
334  qy(j,i) = real(qy_conv(j,i),dp)
335  dzs_dtau(j,i) = real(dzs_dtau_conv(j,i),dp)
336  dzm_dtau(j,i) = real(dzm_dtau_conv(j,i),dp)
337  dzb_dtau(j,i) = real(dzb_dtau_conv(j,i),dp)
338  dzl_dtau(j,i) = real(dzl_dtau_conv(j,i),dp)
339  dh_c_dtau(j,i) = real(dH_c_dtau_conv(j,i),dp)
340  dh_t_dtau(j,i) = real(dH_t_dtau_conv(j,i),dp)
341  vx_b_g(j,i) = real(vx_b_g_conv(j,i),dp)
342  vy_b_g(j,i) = real(vy_b_g_conv(j,i),dp)
343  vz_b(j,i) = real(vz_b_conv(j,i),dp)
344  vx_s_g(j,i) = real(vx_s_g_conv(j,i),dp)
345  vy_s_g(j,i) = real(vy_s_g_conv(j,i),dp)
346  vz_s(j,i) = real(vz_s_conv(j,i),dp)
347  temp_b(j,i) = real(temp_b_conv(j,i),dp)
348  temph_b(j,i) = real(temph_b_conv(j,i),dp)
349  p_b_w(j,i) = real(p_b_w_conv(j,i),dp)
350  h_w(j,i) = real(H_w_conv(j,i),dp)
351  q_gl_g(j,i) = real(q_gl_g_conv(j,i),dp)
352 
353  do kr=0, krmax
354  temp_r(kr,j,i) = real(temp_r_conv(kr,j,i),dp)
355  end do
356 
357  do kt=0, ktmax
358  vx_t(kt,j,i) = real(vx_t_conv(kt,j,i),dp)
359  vy_t(kt,j,i) = real(vy_t_conv(kt,j,i),dp)
360  vz_t(kt,j,i) = real(vz_t_conv(kt,j,i),dp)
361  omega_t(kt,j,i) = real(omega_t_conv(kt,j,i),dp)
362  age_t(kt,j,i) = real(age_t_conv(kt,j,i),dp)
363  end do
364 
365  do kc=0, kcmax
366  vx_c(kc,j,i) = real(vx_c_conv(kc,j,i),dp)
367  vy_c(kc,j,i) = real(vy_c_conv(kc,j,i),dp)
368  vz_c(kc,j,i) = real(vz_c_conv(kc,j,i),dp)
369  temp_c(kc,j,i) = real(temp_c_conv(kc,j,i),dp)
370  age_c(kc,j,i) = real(age_c_conv(kc,j,i),dp)
371  end do
372 
373 end do
374 end do
375 
376 !-------- Set topography of the relaxed bedrock --------
377 
378 zl0 = 0.0_dp
379 
380 !-------- Further stuff --------
381 
382 dxi = dx *1000.0_dp ! km -> m
383 deta = dx *1000.0_dp ! km -> m
384 
385 do i=0, imax
386 do j=0, jmax
387  call geo_coord(phi(j,i), lambda(j,i), xi(i), eta(j))
388 end do
389 end do
390 
391 !-------- Metric tensor, gradients of the topography --------
392 
393 call metric
394 
395 #if TOPOGRAD==0
396 call topograd_1(dxi, deta, 1)
397 #elif TOPOGRAD==1
398 call topograd_2(dxi, deta, 1)
399 #endif
400 
401 !-------- Corresponding area of grid points --------
402 
403 do i=0, imax
404 do j=0, jmax
405  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
406 end do
407 end do
408 
409 contains
410 
411 !-------------------------------------------------------------------------------
412 !> NetCDF error capturing.
413 !<------------------------------------------------------------------------------
414  subroutine check(status)
415  integer(i4b), intent (in) :: status
416  if(status /= nf90_noerr) then
417  write(6,'(a)') trim(nf90_strerror(status))
418  stop ' topography3_nc: Stopped due to netcdf error!'
419  end if
420  end subroutine check
421 
422 end subroutine topography3_nc
423 !
Generated on Mon May 20 2013 10:18:39 for SICOPOLIS V3.0 by   doxygen 1.8.1