SICOPOLIS V3.0
 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 integer(i4b) :: i, j
48 integer(i4b) :: n_switch
49 real(dp) :: dxi, deta, dxi_inv, deta_inv, dxi12_inv, deta12_inv
50 real(dp) :: zs_help(0:jmax,0:imax), zm_help(0:jmax,0:imax), &
51  zb_help(0:jmax,0:imax)
52 
53 dxi_inv = 1.0_dp/dxi
54 deta_inv = 1.0_dp/deta
55 
56 dxi12_inv = 1.0_dp/(12.0_dp*dxi)
57 deta12_inv = 1.0_dp/(12.0_dp*deta)
58 
59 !-------- Distinguish between old and new topography data --------
60 
61 if (n_switch.eq.1) then
62 
63  do i=0, imax
64  do j=0, jmax
65  zs_help(j,i) = zs(j,i)
66  zm_help(j,i) = zm(j,i)
67  zb_help(j,i) = zb(j,i)
68  end do
69  end do
70 
71 else if (n_switch.eq.2) then
72 
73  do i=0, imax
74  do j=0, jmax
75  zs_help(j,i) = zs_neu(j,i)
76  zm_help(j,i) = zm_neu(j,i)
77  zb_help(j,i) = zb_neu(j,i)
78  end do
79  end do
80 
81 else
82  stop ' topograd_2: Wrong value for n_switch!'
83 end if
84 
85 !-------- Topography gradients on the staggered grid --------
86 
87 ! ------ x-derivatives
88 
89 do i=0, imax-1
90 do j=0, jmax
91  dzs_dxi(j,i) = (zs_help(j,i+1)-zs_help(j,i))*dxi_inv &
92  *insq_g11_sgx(j,i)
93  dzm_dxi(j,i) = (zm_help(j,i+1)-zm_help(j,i))*dxi_inv &
94  *insq_g11_sgx(j,i)
95  dzb_dxi(j,i) = (zb_help(j,i+1)-zb_help(j,i))*dxi_inv &
96  *insq_g11_sgx(j,i)
97  dh_c_dxi(j,i) = dzs_dxi(j,i)-dzm_dxi(j,i)
98  dh_t_dxi(j,i) = dzm_dxi(j,i)-dzb_dxi(j,i)
99 end do
100 end do
101 
102 ! ------ y-derivatives
103 
104 do i=0, imax
105 do j=0, jmax-1
106  dzs_deta(j,i) = (zs_help(j+1,i)-zs_help(j,i))*deta_inv &
107  *insq_g22_sgy(j,i)
108  dzm_deta(j,i) = (zm_help(j+1,i)-zm_help(j,i))*deta_inv &
109  *insq_g22_sgy(j,i)
110  dzb_deta(j,i) = (zb_help(j+1,i)-zb_help(j,i))*deta_inv &
111  *insq_g22_sgy(j,i)
112  dh_c_deta(j,i) = dzs_deta(j,i)-dzm_deta(j,i)
113  dh_t_deta(j,i) = dzm_deta(j,i)-dzb_deta(j,i)
114 end do
115 end do
116 
117 !-------- Topography gradients on the grid points --------
118 
119 ! ------ x-derivatives
120 
121 do i=2, imax-2
122 do j=0, jmax
123  dzs_dxi_g(j,i) &
124  = ( -zs_help(j,i+2) + 8.0_dp*zs_help(j,i+1) &
125  -8.0_dp*zs_help(j,i-1) + zs_help(j,i-2) ) &
126  *dxi12_inv &
127  *insq_g11_g(j,i)
128  dzm_dxi_g(j,i) &
129  = ( -zm_help(j,i+2) + 8.0_dp*zm_help(j,i+1) &
130  -8.0_dp*zm_help(j,i-1) + zm_help(j,i-2) ) &
131  *dxi12_inv &
132  *insq_g11_g(j,i)
133  dzb_dxi_g(j,i) &
134  = ( -zb_help(j,i+2) + 8.0_dp*zb_help(j,i+1) &
135  -8.0_dp*zb_help(j,i-1) + zb_help(j,i-2) ) &
136  *dxi12_inv &
137  *insq_g11_g(j,i)
138  dh_c_dxi_g(j,i) = dzs_dxi_g(j,i)-dzm_dxi_g(j,i)
139  dh_t_dxi_g(j,i) = dzm_dxi_g(j,i)-dzb_dxi_g(j,i)
140 end do
141 end do
142 
143 do j=0, jmax
144 
145  dzs_dxi_g(j,0) = (zs_help(j,1)-zs_help(j,0))*dxi_inv &
146  *insq_g11_g(j,0)
147  dzm_dxi_g(j,0) = (zm_help(j,1)-zm_help(j,0))*dxi_inv &
148  *insq_g11_g(j,0)
149  dzb_dxi_g(j,0) = (zb_help(j,1)-zb_help(j,0))*dxi_inv &
150  *insq_g11_g(j,0)
151  dh_c_dxi_g(j,0) = dzs_dxi_g(j,0)-dzm_dxi_g(j,0)
152  dh_t_dxi_g(j,0) = dzm_dxi_g(j,0)-dzb_dxi_g(j,0)
153 
154  dzs_dxi_g(j,1) = (zs_help(j,2)-zs_help(j,0)) &
155  *0.5_dp*dxi_inv &
156  *insq_g11_g(j,1)
157  dzm_dxi_g(j,1) = (zm_help(j,2)-zm_help(j,0)) &
158  *0.5_dp*dxi_inv &
159  *insq_g11_g(j,1)
160  dzb_dxi_g(j,1) = (zb_help(j,2)-zb_help(j,0)) &
161  *0.5_dp*dxi_inv &
162  *insq_g11_g(j,1)
163  dh_c_dxi_g(j,1) = dzs_dxi_g(j,1)-dzm_dxi_g(j,1)
164  dh_t_dxi_g(j,1) = dzm_dxi_g(j,1)-dzb_dxi_g(j,1)
165 
166  dzs_dxi_g(j,imax-1) = (zs_help(j,imax)-zs_help(j,imax-2)) &
167  *0.5_dp*dxi_inv &
168  *insq_g11_g(j,imax-1)
169  dzm_dxi_g(j,imax-1) = (zm_help(j,imax)-zm_help(j,imax-2)) &
170  *0.5_dp*dxi_inv &
171  *insq_g11_g(j,imax-1)
172  dzb_dxi_g(j,imax-1) = (zb_help(j,imax)-zb_help(j,imax-2)) &
173  *0.5_dp*dxi_inv &
174  *insq_g11_g(j,imax-1)
175  dh_c_dxi_g(j,imax-1) = dzs_dxi_g(j,imax-1) &
176  -dzm_dxi_g(j,imax-1)
177  dh_t_dxi_g(j,imax-1) = dzm_dxi_g(j,imax-1) &
178  -dzb_dxi_g(j,imax-1)
179 
180  dzs_dxi_g(j,imax) = (zs_help(j,imax)-zs_help(j,imax-1)) &
181  *dxi_inv &
182  *insq_g11_g(j,imax)
183  dzm_dxi_g(j,imax) = (zm_help(j,imax)-zm_help(j,imax-1)) &
184  *dxi_inv &
185  *insq_g11_g(j,imax)
186  dzb_dxi_g(j,imax) = (zb_help(j,imax)-zb_help(j,imax-1)) &
187  *dxi_inv &
188  *insq_g11_g(j,imax)
189  dh_c_dxi_g(j,imax) = dzs_dxi_g(j,imax)-dzm_dxi_g(j,imax)
190  dh_t_dxi_g(j,imax) = dzm_dxi_g(j,imax)-dzb_dxi_g(j,imax)
191 
192 end do
193 
194 ! ------ y-derivatives
195 
196 do i=0, imax
197 do j=2, jmax-2
198  dzs_deta_g(j,i) &
199  = ( -zs_help(j+2,i) + 8.0_dp*zs_help(j+1,i) &
200  -8.0_dp*zs_help(j-1,i) + zs_help(j-2,i) ) &
201  *deta12_inv &
202  *insq_g22_g(j,i)
203  dzm_deta_g(j,i) &
204  = ( -zm_help(j+2,i) + 8.0_dp*zm_help(j+1,i) &
205  -8.0_dp*zm_help(j-1,i) + zm_help(j-2,i) ) &
206  *deta12_inv &
207  *insq_g22_g(j,i)
208  dzb_deta_g(j,i) &
209  = ( -zb_help(j+2,i) + 8.0_dp*zb_help(j+1,i) &
210  -8.0_dp*zb_help(j-1,i) + zb_help(j-2,i) ) &
211  *deta12_inv &
212  *insq_g22_g(j,i)
213  dh_c_deta_g(j,i) = dzs_deta_g(j,i)-dzm_deta_g(j,i)
214  dh_t_deta_g(j,i) = dzm_deta_g(j,i)-dzb_deta_g(j,i)
215 end do
216 end do
217 
218 do i=0, imax
219 
220  dzs_deta_g(0,i) = (zs_help(1,i)-zs_help(0,i))*deta_inv &
221  *insq_g22_g(0,i)
222  dzm_deta_g(0,i) = (zm_help(1,i)-zm_help(0,i))*deta_inv &
223  *insq_g22_g(0,i)
224  dzb_deta_g(0,i) = (zb_help(1,i)-zb_help(0,i))*deta_inv &
225  *insq_g22_g(0,i)
226  dh_c_deta_g(0,i) = dzs_deta_g(0,i)-dzm_deta_g(0,i)
227  dh_t_deta_g(0,i) = dzm_deta_g(0,i)-dzb_deta_g(0,i)
228 
229  dzs_deta_g(1,i) = (zs_help(2,i)-zs_help(0,i)) &
230  *0.5_dp*deta_inv &
231  *insq_g22_g(1,i)
232  dzm_deta_g(1,i) = (zm_help(2,i)-zm_help(0,i)) &
233  *0.5_dp*deta_inv &
234  *insq_g22_g(1,i)
235  dzb_deta_g(1,i) = (zb_help(2,i)-zb_help(0,i)) &
236  *0.5_dp*deta_inv &
237  *insq_g22_g(1,i)
238  dh_c_deta_g(1,i) = dzs_deta_g(1,i)-dzm_deta_g(1,i)
239  dh_t_deta_g(1,i) = dzm_deta_g(1,i)-dzb_deta_g(1,i)
240 
241  dzs_deta_g(jmax-1,i) = (zs_help(jmax,i)-zs_help(jmax-2,i)) &
242  *0.5_dp*deta_inv &
243  *insq_g22_g(jmax-1,i)
244  dzm_deta_g(jmax-1,i) = (zm_help(jmax,i)-zm_help(jmax-2,i)) &
245  *0.5_dp*deta_inv &
246  *insq_g22_g(jmax-1,i)
247  dzb_deta_g(jmax-1,i) = (zb_help(jmax,i)-zb_help(jmax-2,i)) &
248  *0.5_dp*deta_inv &
249  *insq_g22_g(jmax-1,i)
250  dh_c_deta_g(jmax-1,i) = dzs_deta_g(jmax-1,i) &
251  -dzm_deta_g(jmax-1,i)
252  dh_t_deta_g(jmax-1,i) = dzm_deta_g(jmax-1,i) &
253  -dzb_deta_g(jmax-1,i)
254 
255  dzs_deta_g(jmax,i) = (zs_help(jmax,i)-zs_help(jmax-1,i)) &
256  *deta_inv &
257  *insq_g22_g(jmax,i)
258  dzm_deta_g(jmax,i) = (zm_help(jmax,i)-zm_help(jmax-1,i)) &
259  *deta_inv &
260  *insq_g22_g(jmax,i)
261  dzb_deta_g(jmax,i) = (zb_help(jmax,i)-zb_help(jmax-1,i)) &
262  *deta_inv &
263  *insq_g22_g(jmax,i)
264  dh_c_deta_g(jmax,i) = dzs_deta_g(jmax,i)-dzm_deta_g(jmax,i)
265  dh_t_deta_g(jmax,i) = dzm_deta_g(jmax,i)-dzb_deta_g(jmax,i)
266 
267 end do
268 
269 end subroutine topograd_2
270 !
Generated on Mon May 20 2013 10:18:43 for SICOPOLIS V3.0 by   doxygen 1.8.1