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