SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
topograd_2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : t o p o g r a d _ 2
4 !
5 !> @file
6 !!
7 !! Calculation of topography gradients on the staggered grid and on the grid
8 !! points (the latter by fourth-order discretization).
9 !! Length rescaling with the corresponding components of the metric tensor is
10 !! already included.
11 !!
12 !! @section Copyright
13 !!
14 !! Copyright 2009-2013 Ralf Greve
15 !!
16 !! @section License
17 !!
18 !! This file is part of SICOPOLIS.
19 !!
20 !! SICOPOLIS is free software: you can redistribute it and/or modify
21 !! it under the terms of the GNU General Public License as published by
22 !! the Free Software Foundation, either version 3 of the License, or
23 !! (at your option) any later version.
24 !!
25 !! SICOPOLIS is distributed in the hope that it will be useful,
26 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
27 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 !! GNU General Public License for more details.
29 !!
30 !! You should have received a copy of the GNU General Public License
31 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
32 !<
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 
35 !-------------------------------------------------------------------------------
36 !> Calculation of topography gradients on the staggered grid and on the grid
37 !! points (the latter by fourth-order discretization).
38 !! Length rescaling with the corresponding components of the metric tensor is
39 !! already included.
40 !<------------------------------------------------------------------------------
41 subroutine topograd_2(dxi, deta, n_switch)
42 
43 use sico_types
44 use sico_variables
45 
46 implicit none
47 
48 integer(i4b), intent(in) :: n_switch
49 real(dp), intent(in) :: dxi, deta
50 
51 integer(i4b) :: i, j
52 real(dp) :: dxi_inv, deta_inv, dxi12_inv, deta12_inv
53 real(dp) :: zs_help(0:jmax,0:imax), zm_help(0:jmax,0:imax), &
54  zb_help(0:jmax,0:imax)
55 
56 dxi_inv = 1.0_dp/dxi
57 deta_inv = 1.0_dp/deta
58 
59 dxi12_inv = 1.0_dp/(12.0_dp*dxi)
60 deta12_inv = 1.0_dp/(12.0_dp*deta)
61 
62 !-------- Distinguish between old and new topography data --------
63 
64 if (n_switch == 1) then
65 
66  do i=0, imax
67  do j=0, jmax
68  zs_help(j,i) = zs(j,i)
69  zm_help(j,i) = zm(j,i)
70  zb_help(j,i) = zb(j,i)
71  end do
72  end do
73 
74 else if (n_switch == 2) then
75 
76  do i=0, imax
77  do j=0, jmax
78  zs_help(j,i) = zs_neu(j,i)
79  zm_help(j,i) = zm_neu(j,i)
80  zb_help(j,i) = zb_neu(j,i)
81  end do
82  end do
83 
84 else
85  stop ' topograd_2: Wrong value for n_switch!'
86 end if
87 
88 !-------- Topography gradients on the staggered grid --------
89 
90 ! ------ x-derivatives
91 
92 do i=0, imax-1
93 do j=0, jmax
94  dzs_dxi(j,i) = (zs_help(j,i+1)-zs_help(j,i))*dxi_inv &
95  *insq_g11_sgx(j,i)
96  dzm_dxi(j,i) = (zm_help(j,i+1)-zm_help(j,i))*dxi_inv &
97  *insq_g11_sgx(j,i)
98  dzb_dxi(j,i) = (zb_help(j,i+1)-zb_help(j,i))*dxi_inv &
99  *insq_g11_sgx(j,i)
100  dh_c_dxi(j,i) = dzs_dxi(j,i)-dzm_dxi(j,i)
101  dh_t_dxi(j,i) = dzm_dxi(j,i)-dzb_dxi(j,i)
102 end do
103 end do
104 
105 ! ------ y-derivatives
106 
107 do i=0, imax
108 do j=0, jmax-1
109  dzs_deta(j,i) = (zs_help(j+1,i)-zs_help(j,i))*deta_inv &
110  *insq_g22_sgy(j,i)
111  dzm_deta(j,i) = (zm_help(j+1,i)-zm_help(j,i))*deta_inv &
112  *insq_g22_sgy(j,i)
113  dzb_deta(j,i) = (zb_help(j+1,i)-zb_help(j,i))*deta_inv &
114  *insq_g22_sgy(j,i)
115  dh_c_deta(j,i) = dzs_deta(j,i)-dzm_deta(j,i)
116  dh_t_deta(j,i) = dzm_deta(j,i)-dzb_deta(j,i)
117 end do
118 end do
119 
120 !-------- Topography gradients on the grid points --------
121 
122 ! ------ x-derivatives
123 
124 do i=2, imax-2
125 do j=0, jmax
126  dzs_dxi_g(j,i) &
127  = ( -zs_help(j,i+2) + 8.0_dp*zs_help(j,i+1) &
128  -8.0_dp*zs_help(j,i-1) + zs_help(j,i-2) ) &
129  *dxi12_inv &
130  *insq_g11_g(j,i)
131  dzm_dxi_g(j,i) &
132  = ( -zm_help(j,i+2) + 8.0_dp*zm_help(j,i+1) &
133  -8.0_dp*zm_help(j,i-1) + zm_help(j,i-2) ) &
134  *dxi12_inv &
135  *insq_g11_g(j,i)
136  dzb_dxi_g(j,i) &
137  = ( -zb_help(j,i+2) + 8.0_dp*zb_help(j,i+1) &
138  -8.0_dp*zb_help(j,i-1) + zb_help(j,i-2) ) &
139  *dxi12_inv &
140  *insq_g11_g(j,i)
141  dh_c_dxi_g(j,i) = dzs_dxi_g(j,i)-dzm_dxi_g(j,i)
142  dh_t_dxi_g(j,i) = dzm_dxi_g(j,i)-dzb_dxi_g(j,i)
143 end do
144 end do
145 
146 do j=0, jmax
147 
148  dzs_dxi_g(j,0) = (zs_help(j,1)-zs_help(j,0))*dxi_inv &
149  *insq_g11_g(j,0)
150  dzm_dxi_g(j,0) = (zm_help(j,1)-zm_help(j,0))*dxi_inv &
151  *insq_g11_g(j,0)
152  dzb_dxi_g(j,0) = (zb_help(j,1)-zb_help(j,0))*dxi_inv &
153  *insq_g11_g(j,0)
154  dh_c_dxi_g(j,0) = dzs_dxi_g(j,0)-dzm_dxi_g(j,0)
155  dh_t_dxi_g(j,0) = dzm_dxi_g(j,0)-dzb_dxi_g(j,0)
156 
157  dzs_dxi_g(j,1) = (zs_help(j,2)-zs_help(j,0)) &
158  *0.5_dp*dxi_inv &
159  *insq_g11_g(j,1)
160  dzm_dxi_g(j,1) = (zm_help(j,2)-zm_help(j,0)) &
161  *0.5_dp*dxi_inv &
162  *insq_g11_g(j,1)
163  dzb_dxi_g(j,1) = (zb_help(j,2)-zb_help(j,0)) &
164  *0.5_dp*dxi_inv &
165  *insq_g11_g(j,1)
166  dh_c_dxi_g(j,1) = dzs_dxi_g(j,1)-dzm_dxi_g(j,1)
167  dh_t_dxi_g(j,1) = dzm_dxi_g(j,1)-dzb_dxi_g(j,1)
168 
169  dzs_dxi_g(j,imax-1) = (zs_help(j,imax)-zs_help(j,imax-2)) &
170  *0.5_dp*dxi_inv &
171  *insq_g11_g(j,imax-1)
172  dzm_dxi_g(j,imax-1) = (zm_help(j,imax)-zm_help(j,imax-2)) &
173  *0.5_dp*dxi_inv &
174  *insq_g11_g(j,imax-1)
175  dzb_dxi_g(j,imax-1) = (zb_help(j,imax)-zb_help(j,imax-2)) &
176  *0.5_dp*dxi_inv &
177  *insq_g11_g(j,imax-1)
178  dh_c_dxi_g(j,imax-1) = dzs_dxi_g(j,imax-1) &
179  -dzm_dxi_g(j,imax-1)
180  dh_t_dxi_g(j,imax-1) = dzm_dxi_g(j,imax-1) &
181  -dzb_dxi_g(j,imax-1)
182 
183  dzs_dxi_g(j,imax) = (zs_help(j,imax)-zs_help(j,imax-1)) &
184  *dxi_inv &
185  *insq_g11_g(j,imax)
186  dzm_dxi_g(j,imax) = (zm_help(j,imax)-zm_help(j,imax-1)) &
187  *dxi_inv &
188  *insq_g11_g(j,imax)
189  dzb_dxi_g(j,imax) = (zb_help(j,imax)-zb_help(j,imax-1)) &
190  *dxi_inv &
191  *insq_g11_g(j,imax)
192  dh_c_dxi_g(j,imax) = dzs_dxi_g(j,imax)-dzm_dxi_g(j,imax)
193  dh_t_dxi_g(j,imax) = dzm_dxi_g(j,imax)-dzb_dxi_g(j,imax)
194 
195 end do
196 
197 ! ------ y-derivatives
198 
199 do i=0, imax
200 do j=2, jmax-2
201  dzs_deta_g(j,i) &
202  = ( -zs_help(j+2,i) + 8.0_dp*zs_help(j+1,i) &
203  -8.0_dp*zs_help(j-1,i) + zs_help(j-2,i) ) &
204  *deta12_inv &
205  *insq_g22_g(j,i)
206  dzm_deta_g(j,i) &
207  = ( -zm_help(j+2,i) + 8.0_dp*zm_help(j+1,i) &
208  -8.0_dp*zm_help(j-1,i) + zm_help(j-2,i) ) &
209  *deta12_inv &
210  *insq_g22_g(j,i)
211  dzb_deta_g(j,i) &
212  = ( -zb_help(j+2,i) + 8.0_dp*zb_help(j+1,i) &
213  -8.0_dp*zb_help(j-1,i) + zb_help(j-2,i) ) &
214  *deta12_inv &
215  *insq_g22_g(j,i)
216  dh_c_deta_g(j,i) = dzs_deta_g(j,i)-dzm_deta_g(j,i)
217  dh_t_deta_g(j,i) = dzm_deta_g(j,i)-dzb_deta_g(j,i)
218 end do
219 end do
220 
221 do i=0, imax
222 
223  dzs_deta_g(0,i) = (zs_help(1,i)-zs_help(0,i))*deta_inv &
224  *insq_g22_g(0,i)
225  dzm_deta_g(0,i) = (zm_help(1,i)-zm_help(0,i))*deta_inv &
226  *insq_g22_g(0,i)
227  dzb_deta_g(0,i) = (zb_help(1,i)-zb_help(0,i))*deta_inv &
228  *insq_g22_g(0,i)
229  dh_c_deta_g(0,i) = dzs_deta_g(0,i)-dzm_deta_g(0,i)
230  dh_t_deta_g(0,i) = dzm_deta_g(0,i)-dzb_deta_g(0,i)
231 
232  dzs_deta_g(1,i) = (zs_help(2,i)-zs_help(0,i)) &
233  *0.5_dp*deta_inv &
234  *insq_g22_g(1,i)
235  dzm_deta_g(1,i) = (zm_help(2,i)-zm_help(0,i)) &
236  *0.5_dp*deta_inv &
237  *insq_g22_g(1,i)
238  dzb_deta_g(1,i) = (zb_help(2,i)-zb_help(0,i)) &
239  *0.5_dp*deta_inv &
240  *insq_g22_g(1,i)
241  dh_c_deta_g(1,i) = dzs_deta_g(1,i)-dzm_deta_g(1,i)
242  dh_t_deta_g(1,i) = dzm_deta_g(1,i)-dzb_deta_g(1,i)
243 
244  dzs_deta_g(jmax-1,i) = (zs_help(jmax,i)-zs_help(jmax-2,i)) &
245  *0.5_dp*deta_inv &
246  *insq_g22_g(jmax-1,i)
247  dzm_deta_g(jmax-1,i) = (zm_help(jmax,i)-zm_help(jmax-2,i)) &
248  *0.5_dp*deta_inv &
249  *insq_g22_g(jmax-1,i)
250  dzb_deta_g(jmax-1,i) = (zb_help(jmax,i)-zb_help(jmax-2,i)) &
251  *0.5_dp*deta_inv &
252  *insq_g22_g(jmax-1,i)
253  dh_c_deta_g(jmax-1,i) = dzs_deta_g(jmax-1,i) &
254  -dzm_deta_g(jmax-1,i)
255  dh_t_deta_g(jmax-1,i) = dzm_deta_g(jmax-1,i) &
256  -dzb_deta_g(jmax-1,i)
257 
258  dzs_deta_g(jmax,i) = (zs_help(jmax,i)-zs_help(jmax-1,i)) &
259  *deta_inv &
260  *insq_g22_g(jmax,i)
261  dzm_deta_g(jmax,i) = (zm_help(jmax,i)-zm_help(jmax-1,i)) &
262  *deta_inv &
263  *insq_g22_g(jmax,i)
264  dzb_deta_g(jmax,i) = (zb_help(jmax,i)-zb_help(jmax-1,i)) &
265  *deta_inv &
266  *insq_g22_g(jmax,i)
267  dh_c_deta_g(jmax,i) = dzs_deta_g(jmax,i)-dzm_deta_g(jmax,i)
268  dh_t_deta_g(jmax,i) = dzm_deta_g(jmax,i)-dzb_deta_g(jmax,i)
269 
270 end do
271 
272 end subroutine topograd_2
273 !
Generated on Thu Sep 12 2013 09:45:20 for SICOPOLIS V3.1 by   doxygen 1.8.2