da_setup_runconstants.inc
References to this file elsewhere.
1 subroutine da_setup_runconstants(grid, xbx)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Define constants used.
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (domain), intent(inout) :: grid
10 type (xbx_type), intent(inout) :: xbx ! Header & non-gridded vars.
11
12 integer, parameter :: nrange = 50 ! Range to search for efficient FFT.
13 integer :: n ! Loop counter.
14
15 integer :: fft_size_i ! Efficient FFT 1st dimension.
16 integer :: fft_size_j ! Efficient FFT 2nd dimension.
17 logical :: found_magic ! true when efficient FFT found
18 logical :: need_pad ! True when need pad_i > 0
19
20 integer :: i, j
21 integer :: nx ! nx + 1 = ix + pad_i
22 integer :: ny ! ny + 1 = jy + pad_j
23 real :: const ! Multiplicative constants.
24 real :: coeff_nx ! Multiplicative coefficients.
25 real :: coeff_ny ! Multiplicative coefficients.
26 real :: cos_coeff_nx ! Multiplicative coefficients.
27 real :: cos_coeff_ny ! Multiplicative coefficients.
28
29 if (trace_use) call da_trace_entry("da_setup_runconstants")
30
31 !---------------------------------------------------------------------------
32 ! [1.0] Calculate padded grid-size required for efficient FFTs and factors:
33 !---------------------------------------------------------------------------
34
35 ! [1.1] In x direction:
36
37 nx = ide - ids
38
39 allocate(xbx % fft_factors_x(num_fft_factors))
40
41 do n = nx, nx+nrange
42 call da_find_fft_factors(n, found_magic, xbx % fft_factors_x)
43
44 if (found_magic) then
45 if (mod(n, 2) == 0) then
46 fft_size_i = n
47 xbx % fft_pad_i = n - nx
48 exit
49 end if
50 end if
51 end do
52
53 if (.NOT. found_magic) then
54 call da_error(__FILE__,__LINE__, &
55 (/"No factor found"/))
56 end if
57
58 allocate(xbx % trig_functs_x(1:3*fft_size_i))
59
60 call da_find_fft_trig_funcs(fft_size_i, xbx % trig_functs_x(:))
61
62 ! [1.2] In y direction:
63
64 ny = jde - jds
65
66 allocate(xbx % fft_factors_y(num_fft_factors))
67
68 do n = ny, ny+nrange
69 call da_find_fft_factors(n, found_magic, xbx % fft_factors_y)
70
71 if (found_magic) then
72 if (mod(n, 2) == 0) then
73 fft_size_j = n
74 xbx % fft_pad_j = n - ny
75 exit
76 end if
77 end if
78 end do
79
80 if (.NOT. found_magic) then
81 call da_error(__FILE__,__LINE__, &
82 (/"No factor found"/))
83 end if
84
85 allocate(xbx % trig_functs_y(1:3*fft_size_j))
86
87 call da_find_fft_trig_funcs(fft_size_j, xbx % trig_functs_y(:))
88
89 !-----Multiplicative coefficent for solution of spectral Poission eqn:
90
91 !mslee.Bgrid
92 ! const = -0.5 * grid%xb%ds * grid%xb%ds
93 const = -2.0 * grid%xb%ds * grid%xb%ds
94
95 nx = ide - ids + xbx%fft_pad_i
96 ny = jde - jds + xbx%fft_pad_j
97 ! YRG: A-grid:
98 coeff_nx = 2.0 * pi / real(nx)
99 coeff_ny = 2.0 * pi / real(ny)
100
101 ! YRG: B-grid::
102 ! coeff_nx = pi / real(nx)
103 ! coeff_ny = pi / real(ny)
104
105 xbx%fft_ix = nx + 1
106 xbx%fft_jy = ny + 1
107
108 allocate(xbx%fft_coeffs(1:xbx%fft_ix,1:xbx%fft_jy))
109
110 do j = 2, ny
111 cos_coeff_ny = COS(coeff_ny * real(j-1))
112 do i = 2, nx
113 cos_coeff_nx = COS(coeff_nx * real(i-1))
114 !mslee.Bgrid
115 ! xbx%fft_coeffs(i,j) = const / (1.0 - cos_coeff_nx * cos_coeff_ny)
116 if (cos_coeff_nx.eq.1.and.cos_coeff_ny.eq.1) then
117 xbx%fft_coeffs(i,j) = 0.0
118 else
119 xbx%fft_coeffs(i,j) = const / (2.0 - cos_coeff_nx - cos_coeff_ny)
120 end if
121 end do
122 end do
123
124 ! Set to zero all coefficients which are multiplied by sin(0.0) in FST:
125
126 !mslee xbx%fft_coeffs(1,1:xbx%fft_jy) = 0.0
127 !mslee xbx%fft_coeffs(xbx%fft_ix,1:xbx%fft_jy) = 0.0
128 !mslee xbx%fft_coeffs(1:xbx%fft_ix,1) = 0.0
129 !mslee xbx%fft_coeffs(1:xbx%fft_ix,xbx%fft_jy) = 0.0
130 !mslee. we need to check the following
131
132 xbx%fft_adjoint_factor = 4.0 / real(nx * ny)
133
134 !---------------------------------------------------------------------------
135
136 ! Calculate i increment for distributing pad region across processors.
137
138 need_pad = .true.
139 if (xbx%fft_pad_i <= 0) then
140 need_pad = .false.
141 else if (xbx%fft_pad_i > (ide-ids+1)) then
142 write(unit=message(1), fmt='(a)') &
143 'FFT xbx%fft_pad_i is too large!'
144 write(unit=message(2), fmt='(2x,a,i4)') &
145 '(ide-ids+1) = ', (ide-ids+1), &
146 'xbx%fft_pad_i = ', xbx%fft_pad_i
147 call da_error (__FILE__,__LINE__,message(1:2))
148 else
149 xbx%pad_inc = (ide-ids+1)/xbx%fft_pad_i
150 end if
151
152 ! Calculate number of pad vectors in x to be done on this processor.
153 ! Need to save xbx%pad_num, xbx%pad_inc, and xbx%pad_loc
154
155 xbx%pad_num = 0
156 if (need_pad) then
157 do n=1, xbx%fft_pad_i
158 i = (n-1)*xbx%pad_inc + 1
159 if ((i >= grid%xp%itsy) .and. (i <= grid%xp%itey)) then
160 xbx%pad_num = xbx%pad_num + 1
161 end if
162 end do
163
164 if (xbx%pad_num > 0) then
165 allocate(xbx%pad_loc(1:xbx%pad_num))
166 allocate(xbx%pad_pos(1:xbx%pad_num))
167 end if
168
169 xbx%pad_num = 0
170 do n=1, xbx%fft_pad_i
171 i = (n-1)*xbx%pad_inc + 1
172 if ((i >= grid%xp%itsy) .and. (i <= grid%xp%itey)) then
173 xbx%pad_num = xbx%pad_num + 1
174 xbx%pad_loc(xbx%pad_num) = i
175 xbx%pad_pos(xbx%pad_num) = grid%xp%ide + n
176 end if
177 end do
178 end if
179
180 !---------------------------------------------------------------------------
181
182 if (global) then
183 ! Set up Spectral transform constants:
184 xbx%inc = 1
185 xbx%ni = ide-ids+1
186 xbx%nj = jde-jds+1
187 xbx%lenr = xbx%inc * (xbx%ni - 1) + 1
188 xbx%lensav = xbx%ni + int(log(real(xbx%ni))) + 4
189 xbx%lenwrk = xbx%ni
190 xbx%max_wavenumber = xbx%ni/2 -1
191 xbx%alp_size = (xbx%nj+ 1)*(xbx%max_wavenumber+1)*(xbx%max_wavenumber+2)/4
192
193 if (print_detail_spectral) then
194 write (unit=stdout,fmt='(a)') "Spectral transform constants"
195 write (unit=stdout, fmt='(a, i8)') &
196 ' inc =', xbx%inc , &
197 ' ni =', xbx%ni , &
198 ' nj =', xbx%nj , &
199 ' lenr =', xbx%lenr , &
200 ' lensav =', xbx%lensav, &
201 ' lenwrk =', xbx%lenwrk, &
202 ' max_wavenumber =', xbx%max_wavenumber, &
203 ' alp_size =', xbx%alp_size
204 write (unit=stdout,fmt='(a)') " "
205 end if
206
207 allocate(xbx%coslat(1:xbx%nj))
208 allocate(xbx%sinlat(1:xbx%nj))
209 allocate(xbx%coslon(1:xbx%ni))
210 allocate(xbx%sinlon(1:xbx%ni))
211 allocate(xbx%int_wgts(1:xbx%nj)) ! Interpolation weights
212 allocate(xbx%alp(1:xbx%alp_size))
213 allocate(xbx%wsave(1:xbx%lensav))
214
215 call da_initialize_h(xbx%ni, xbx%nj, xbx%max_wavenumber, &
216 xbx%lensav, xbx%alp_size, &
217 xbx%wsave, grid%xb%lon, xbx%sinlon, &
218 xbx%coslon, grid%xb%lat, xbx%sinlat, &
219 xbx%coslat, xbx%int_wgts, xbx%alp)
220
221
222 end if
223
224 if (trace_use) call da_trace_exit("da_setup_runconstants")
225
226 end subroutine da_setup_runconstants
227
228