gen_be_stage0_wrf.f90
References to this file elsewhere.
1 program gen_be_stage0_wrf
2 !
3 !----------------------------------------------------------------------
4 ! Purpose : To convert WRF output to "standard perturbation fields"
5 ! required by the WRF-Var BE covariance generation utility "gen_be".
6 !
7 ! Input : WRF forecasts for a specified date (NETCDF format).
8 !
9 ! Output : Binary files for use in gen_be_stage1.
10 !
11 ! Author: Dale Barker (NCAR/MMM)
12 ! Please acknowledge author/institute in work that uses this code.
13 !
14 !----------------------------------------------------------------------
15
16 #ifdef crayx1
17 #define iargc ipxfargc
18 #endif
19
20 use da_control, only : num_fft_factors, pi, stdout, stderr, &
21 filename_len, base_pres, base_temp, base_lapse
22 use da_gen_be, only : da_get_field, da_get_height, da_get_trh, &
23 da_stage0_initialize
24 use da_tools_serial, only : da_get_unit, da_free_unit,da_find_fft_factors, &
25 da_find_fft_trig_funcs
26 use module_ffts, only : fft551, fft661
27
28 implicit none
29
30 integer :: gen_be_iunit, gen_be_ounit
31
32 integer, parameter :: nrange = 50 ! Range to search for efficient FFT.
33
34 character (len=filename_len) :: filestub ! General filename stub.
35 character (len=filename_len) :: input_file ! Input file.
36 character (len=filename_len) :: output_file ! Output file.
37 character (len=10) :: date ! Character date.
38 character (len=10) :: var ! Variable to search for.
39 character (len=3) :: be_method ! "NMC" or "ENS".
40 character (len=3) :: cne ! Ensemble size.
41 character (len=3) :: ce ! Member index -> character.
42
43 integer, external :: iargc
44 integer :: numarg
45 integer :: ne ! Ensemble size.
46 integer :: i, j, k, member ! Loop counters.
47 integer :: dim1 ! Dimensions of grid (T points).
48 integer :: dim1s ! Dimensions of grid (vor/psi pts).
49 integer :: dim2 ! Dimensions of grid (T points).
50 integer :: dim2s ! Dimensions of grid (vor/psi pts).
51 integer :: dim3 ! Dimensions of grid (T points).
52 integer :: n1, n2 ! Padded dimensions (n=dim-1+pad).
53 integer :: n1s, n2s ! Padded dimensions (n=dim-1+pad).
54 integer :: poisson_method ! 1=Spectral, 2=SOR.
55 integer :: fft_method ! For poisson_method=1: 1=FCT, 2=FST.
56 integer :: ktest ! Test level.
57 real :: member_inv ! 1 / member.
58 real :: ds ! Grid resolution.
59 real :: dim12_inv_u ! 1 / (dim1*dim2).
60 real :: dim12_inv_v ! 1 / (dim1*dim2).
61 real :: dim12_inv ! 1 / (dim1*dim2).
62 logical :: test_inverse ! Test conversion by performing inverse.
63
64 integer :: ifax1(1:num_fft_factors) ! FFT factors.
65 integer :: ifax2(1:num_fft_factors) ! FFT factors.
66 integer :: ifax1s(1:num_fft_factors) ! FFT factors.
67 integer :: ifax2s(1:num_fft_factors) ! FFT factors.
68 real, allocatable :: xlat(:,:) ! Latitude of mass points.
69 real, allocatable :: mapfac_m(:,:) ! Map factor - mass pts.
70 real, allocatable :: mapfac_u(:,:) ! Map factor - u points.
71 real, allocatable :: mapfac_v(:,:) ! Map factor - v points.
72
73 real, allocatable :: u(:,:) ! u-wind.
74 real, allocatable :: v(:,:) ! v-wind.
75 real, allocatable :: div(:,:) ! Divergence.
76 real, allocatable :: vor(:,:) ! Vorticity.
77 real, allocatable :: psi2d(:,:) ! Streamfunction copy.
78 real, allocatable :: chi2d(:,:) ! Velocity potential copy.
79 real, allocatable :: temp2d(:,:) ! Temperature.
80 real, allocatable :: rh2d(:,:) ! Relative humidity.
81
82 real, allocatable :: trigs1(:) ! FFT trig functions.
83 real, allocatable :: trigs2(:) ! FFT trig functions.
84 real, allocatable :: fft_coeffs(:,:) ! FFT coefficients.
85 real, allocatable :: trigs1s(:) ! FFT trig functions.
86 real, allocatable :: trigs2s(:) ! FFT trig functions.
87 real, allocatable :: fft_coeffss(:,:) ! FFT coefficients.
88
89 ! Standard fields:
90 real, allocatable :: psi(:,:,:) ! Streamfunction.
91 real, allocatable :: chi(:,:,:) ! Velocity Potential.
92 real, allocatable :: temp(:,:,:) ! Temperature.
93 real, allocatable :: rh(:,:,:) ! Relative humidity.
94 real, allocatable :: psfc(:,:) ! Surface pressure.
95 real, allocatable :: height(:,:,:) ! Height.
96 real, allocatable :: psi_mean(:,:,:) ! Streamfunction.
97 real, allocatable :: chi_mean(:,:,:) ! Velocity Potential.
98 real, allocatable :: temp_mean(:,:,:) ! Temperature.
99 real, allocatable :: rh_mean(:,:,:) ! Relative humidity.
100 real, allocatable :: psfc_mean(:,:) ! Surface pressure.
101
102
103 stderr = 0
104 stdout = 6
105
106 base_pres=100000.0
107 base_temp=290.0
108 base_lapse=50.0
109
110 write(6,'(/a)')' [1] Initialize information.'
111
112 call da_get_unit(gen_be_iunit)
113 call da_get_unit(gen_be_ounit)
114
115 test_inverse = .true.
116 ktest = 1
117 poisson_method = 1
118 fft_method = 2
119
120 numarg = iargc()
121 if ( numarg /= 4 )then
122 write(UNIT=6,FMT='(a)') &
123 "Usage: gen_be_stage0_wrf be_method date ne <wrf_file_stub> Stop"
124 stop
125 end if
126
127 ! Initialise to stop false Cray compiler warnings
128 be_method=""
129 date=""
130 cne=""
131 filestub=""
132
133 call getarg( 1, be_method )
134 call getarg( 2, date )
135 call getarg( 3, cne )
136 read(cne,'(i3)')ne
137 call getarg( 4, filestub )
138
139 if ( be_method == "NMC" ) then
140 write(6,'(a,a)')' Computing gen_be NMC-method forecast difference files for date ', date
141 ne = 2 ! NMC-method uses differences between 2 forecasts.
142 else if ( be_method == "ENS" ) then
143 write(6,'(a,a)')' Computing gen_be ensemble perturbation files for date ', date
144 write(6,'(a,i4)')' Ensemble Size = ', ne
145 else
146 write(6,'(a,a)')trim(be_method), ' is an invalid value of be_method. Stop'
147 stop
148 end if
149 write(6,'(a,a)')' Input filestub = ', trim(filestub)
150
151 !---------------------------------------------------------------------------------------------
152 write(6,'(/a)')' [2] Setup up ancillary fields using 1st member values.'
153 !---------------------------------------------------------------------------------------------
154
155 var = "T"
156 input_file = trim(filestub)//'.e001'
157 call da_stage0_initialize( input_file, var, dim1, dim2, dim3, ds )
158 dim1s = dim1+1 ! Vorticity/streamfunction array 1 larger.
159 dim2s = dim2+1 ! Vorticity/streamfunction array 1 larger.
160 dim12_inv_u = 1.0 / real((dim1+1) * dim2)
161 dim12_inv_v = 1.0 / real(dim1 * (dim2+1))
162 dim12_inv = 1.0 / real(dim1 * dim2)
163
164 allocate( xlat(1:dim1,1:dim2) )
165 allocate( mapfac_m(1:dim1,1:dim2) )
166 allocate( mapfac_u(1:dim1+1,1:dim2) )
167 allocate( mapfac_v(1:dim1,1:dim2+1) )
168
169 var = "XLAT"
170 call da_get_field( input_file, var, 2, dim1, dim2, 1, 1, xlat )
171 var = "MAPFAC_M"
172 call da_get_field( input_file, var, 2, dim1, dim2, 1, 1, mapfac_m )
173 var = "MAPFAC_U"
174 call da_get_field( input_file, var, 2, dim1+1, dim2, 1, 1, mapfac_u )
175 var = "MAPFAC_V"
176 call da_get_field( input_file, var, 2, dim1, dim2+1, 1, 1, mapfac_v )
177
178 ! Initialize FFT coefficients:
179
180 if ( poisson_method == 1 ) then
181 call da_fft_initialize1( dim1, dim2, n1, n2, ifax1, ifax2 )
182 call da_fft_initialize1( dim1s, dim2s, n1s, n2s, ifax1s, ifax2s )
183
184 allocate( trigs1(1:3*n1) )
185 allocate( trigs2(1:3*n2) )
186 allocate( fft_coeffs(1:n1+1,1:n2+1) )
187 call da_fft_initialize2( n1, n2, ds, trigs1, trigs2, fft_coeffs )
188 allocate( trigs1s(1:3*n1s) )
189 allocate( trigs2s(1:3*n2s) )
190 allocate( fft_coeffss(1:n1s+1,1:n2s+1) )
191 call da_fft_initialize2( n1s, n2s, ds, trigs1s, trigs2s, fft_coeffss )
192 end if
193
194 ! Allocate arrays in output fields:
195 allocate( psi(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
196 allocate( chi(1:dim1,1:dim2,1:dim3) )
197 allocate( temp(1:dim1,1:dim2,1:dim3) )
198 allocate( rh(1:dim1,1:dim2,1:dim3) )
199 allocate( psfc(1:dim1,1:dim2) )
200 allocate( height(1:dim1,1:dim2,1:dim3) )
201
202 allocate( psi_mean(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
203 allocate( chi_mean(1:dim1,1:dim2,1:dim3) )
204 allocate( temp_mean(1:dim1,1:dim2,1:dim3) )
205 allocate( rh_mean(1:dim1,1:dim2,1:dim3) )
206 allocate( psfc_mean(1:dim1,1:dim2) )
207 psi_mean = 0.0
208 chi_mean = 0.0
209 temp_mean = 0.0
210 rh_mean = 0.0
211 psfc_mean = 0.0
212
213 !---------------------------------------------------------------------------------------------
214 write(6,'(/a)')' [3] Convert WRF forecast fields to standard fields and output'
215 !---------------------------------------------------------------------------------------------
216
217 ! Allocate temporary arrays:
218 allocate( u(1:dim1s,1:dim2) )
219 allocate( v(1:dim1,1:dim2s) )
220 allocate( vor(1:dim1s,1:dim2s) )
221 allocate( div(1:dim1,1:dim2) )
222 allocate( psi2d(1:dim1s,1:dim2s) )
223 allocate( chi2d(1:dim1,1:dim2) )
224 allocate( temp2d(1:dim1,1:dim2) )
225 allocate( rh2d(1:dim1,1:dim2) )
226
227 do member = 1, ne
228
229 write(UNIT=ce,FMT='(i3.3)')member
230 input_file = trim(filestub)//'.e'//ce
231
232 do k = 1, dim3
233
234 ! Read u, v:
235 var = "U"
236 call da_get_field( input_file, var, 3, dim1s, dim2, dim3, k, u )
237 var = "V"
238 call da_get_field( input_file, var, 3, dim1, dim2s, dim3, k, v )
239
240 ! Calculate vorticity (in center of mass grid on WRF's Arakawa C-grid):
241 call da_uv_to_vor_c( dim1, dim2, ds, &
242 mapfac_m, mapfac_u, mapfac_v, u, v, vor )
243
244 ! Calculate divergence (at mass pts. on WRF's Arakawa C-grid):
245
246 call da_uv_to_div_c( dim1, dim2, ds, &
247 mapfac_m, mapfac_u, mapfac_v, u, v, div )
248
249 ! Calculate streamfunction and potential
250 ! Assumes vor/div converted to Del**2 psi/chi):
251
252 if ( poisson_method == 1 ) then
253 call da_del2a_to_a( dim1s, dim2s, n1s, n2s, fft_method, ifax1s, ifax2s, &
254 trigs1s, trigs2s, fft_coeffss, vor, psi2d )
255
256 call da_del2a_to_a( dim1, dim2, n1, n2, fft_method, ifax1, ifax2, &
257 trigs1, trigs2, fft_coeffs, div, chi2d )
258 else if ( poisson_method == 2 ) then
259 call da_sor( dim1s, dim2s, ds, vor, psi2d )
260 call da_sor( dim1, dim2, ds, div, chi2d )
261 end if
262
263 if ( test_inverse .and. k == ktest .and. member == 1 ) then
264 call da_test_inverse( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
265 n1, n2, fft_method, ifax1, ifax2, trigs1, trigs2, fft_coeffs, &
266 n1s, n2s, ifax1s, ifax2s, trigs1s, trigs2s, fft_coeffss, &
267 u, v, psi2d, chi2d )
268 end if
269
270 ! Interpolate psi to mass pts ready for output:
271 do j = 1, dim2
272 do i = 1, dim1
273 psi(i,j,k) = 0.25 * ( psi2d(i,j) + psi2d(i+1,j) + &
274 psi2d(i,j+1) + psi2d(i+1,j+1) )
275 end do
276 end do
277 chi(:,:,k) = chi2d(:,:)
278
279 ! Read mass fields, and convert to T and rh:
280
281 call da_get_trh( input_file, dim1, dim2, dim3, k, temp2d, rh2d )
282 temp(:,:,k) = temp2d(:,:)
283 rh(:,:,k) = 0.01 * rh2d(:,:) ! *0.01 to conform with WRF-Var units.
284
285 end do
286
287 call da_get_height( input_file, dim1, dim2, dim3, height )
288
289 ! Finally, extract surface pressure:
290 var = "PSFC"
291 call da_get_field( input_file, var, 2, dim1, dim2, dim3, 1, psfc )
292
293 ! Write out ensemble forecasts for this member:
294 output_file = 'tmp.e'//ce
295 open (gen_be_ounit, file = output_file, form='unformatted')
296 write(gen_be_ounit)date, dim1, dim2, dim3
297 write(gen_be_ounit)psi
298 write(gen_be_ounit)chi
299 write(gen_be_ounit)temp
300 write(gen_be_ounit)rh
301 write(gen_be_ounit)psfc
302 write(gen_be_ounit)height
303 write(gen_be_ounit)xlat
304 close(gen_be_ounit)
305
306 ! Calculate accumulating mean:
307 member_inv = 1.0 / real(member)
308
309 psi_mean = ( real( member-1 ) * psi_mean + psi ) * member_inv
310 chi_mean = ( real( member-1 ) * chi_mean + chi ) * member_inv
311 temp_mean = ( real( member-1 ) * temp_mean + temp ) * member_inv
312 rh_mean = ( real( member-1 ) * rh_mean + rh ) * member_inv
313 psfc_mean = ( real( member-1 ) * psfc_mean + psfc ) * member_inv
314
315 end do
316
317 deallocate( u )
318 deallocate( v )
319 deallocate( vor )
320 deallocate( div )
321 deallocate( psi2d )
322 deallocate( chi2d )
323 deallocate( temp2d )
324 deallocate( rh2d )
325
326 !---------------------------------------------------------------------------------------------
327 write(6,'(/a)')' [4] Compute perturbations and output'
328 !---------------------------------------------------------------------------------------------
329
330 if ( be_method == "NMC" ) then
331 write(6,'(/a)')' Compute perturbation as a difference between two forecasts'
332
333 ! Re-read input forecast standard fields (ne=2 hard-wired above for NMC-method):
334 input_file = 'tmp.e001'
335 open (gen_be_iunit, file = input_file, form='unformatted')
336 read(gen_be_iunit)date, dim1, dim2, dim3
337 read(gen_be_iunit)psi
338 read(gen_be_iunit)chi
339 read(gen_be_iunit)temp
340 read(gen_be_iunit)rh
341 read(gen_be_iunit)psfc
342 read(gen_be_iunit)height
343 read(gen_be_iunit)xlat
344 close(gen_be_iunit)
345 call da_free_unit(gen_be_iunit)
346
347 ! Note overwriting mean diagnostic to save memory!
348 input_file = 'tmp.e002'
349 open (gen_be_iunit, file = input_file, form='unformatted')
350 read(gen_be_iunit)date, dim1, dim2, dim3
351 read(gen_be_iunit)psi_mean
352 read(gen_be_iunit)chi_mean
353 read(gen_be_iunit)temp_mean
354 read(gen_be_iunit)rh_mean
355 read(gen_be_iunit)psfc_mean
356 read(gen_be_iunit)height
357 read(gen_be_iunit)xlat
358 close(gen_be_iunit)
359
360 ! Take forecast difference:
361 psi = psi - psi_mean
362 chi = chi - chi_mean
363 temp = temp - temp_mean
364 rh = rh - rh_mean
365 psfc = psfc - psfc_mean
366
367 ! Write out NMC-method standard perturbations:
368 output_file = 'pert.'//date(1:10)//'.e001'
369 open (gen_be_ounit, file = output_file, form='unformatted')
370 write(gen_be_ounit)date, dim1, dim2, dim3
371 write(gen_be_ounit)psi
372 write(gen_be_ounit)chi
373 write(gen_be_ounit)temp
374 write(gen_be_ounit)rh
375 write(gen_be_ounit)psfc
376 write(gen_be_ounit)height
377 write(gen_be_ounit)xlat
378 close(gen_be_ounit)
379
380 ! Restore mean (1 member only for NMC-method):
381 psi_mean = psi
382 chi_mean = chi
383 temp_mean = temp
384 rh_mean = rh
385 psfc_mean = psfc
386
387 else ! be_method = "ENS"
388
389 write(6,'(/a)') " [4.1] Convert ensemble of standard fields to perturbations"
390
391 do member = 1, ne
392 write(UNIT=ce,FMT='(i3.3)')member
393
394 ! Re-read ensemble member standard fields:
395 input_file = 'tmp.e'//ce
396 open (gen_be_iunit, file = input_file, form='unformatted')
397 read(gen_be_iunit)date, dim1, dim2, dim3
398 read(gen_be_iunit)psi
399 read(gen_be_iunit)chi
400 read(gen_be_iunit)temp
401 read(gen_be_iunit)rh
402 read(gen_be_iunit)psfc
403 read(gen_be_iunit)height
404 read(gen_be_iunit)xlat
405 close(gen_be_iunit)
406
407 psi = psi - psi_mean
408 chi = chi - chi_mean
409 temp = temp - temp_mean
410 rh = rh - rh_mean
411 psfc = psfc - psfc_mean
412
413 ! Write out standard perturbations for this member:
414 output_file = 'pert.'//date(1:10)//'.e'//ce
415 open (gen_be_ounit, file = output_file, form='unformatted')
416 write(gen_be_ounit)date, dim1, dim2, dim3
417 write(gen_be_ounit)psi
418 write(gen_be_ounit)chi
419 write(gen_be_ounit)temp
420 write(gen_be_ounit)rh
421 write(gen_be_ounit)psfc
422 write(gen_be_ounit)height ! Full height
423 write(gen_be_ounit)xlat ! Full latitude
424 close(gen_be_ounit)
425 end do
426
427 end if
428
429 ! Write out mean/mean square fields:
430
431 output_file = 'psi/'//date(1:10)//'.psi.mean'
432 open (gen_be_ounit, file = output_file, form='unformatted')
433 write(gen_be_ounit)dim1, dim2, dim3
434 write(gen_be_ounit)psi_mean
435 close(gen_be_ounit)
436
437 output_file = 'chi/'//date(1:10)//'.chi.mean'
438 open (gen_be_ounit, file = output_file, form='unformatted')
439 write(gen_be_ounit)dim1, dim2, dim3
440 write(gen_be_ounit)chi_mean
441 close(gen_be_ounit)
442
443 output_file = 't/'//date(1:10)//'.t.mean'
444 open (gen_be_ounit, file = output_file, form='unformatted')
445 write(gen_be_ounit)dim1, dim2, dim3
446 write(gen_be_ounit)temp_mean
447 close(gen_be_ounit)
448
449 output_file = 'rh/'//date(1:10)//'.rh.mean'
450 open (gen_be_ounit, file = output_file, form='unformatted')
451 write(gen_be_ounit)dim1, dim2, dim3
452 write(gen_be_ounit)rh_mean
453 close(gen_be_ounit)
454
455 output_file = 'ps/'//date(1:10)//'.ps.mean'
456 open (gen_be_ounit, file = output_file, form='unformatted')
457 write(gen_be_ounit)dim1, dim2, dim3
458 write(gen_be_ounit)psfc_mean
459 close(gen_be_ounit)
460
461 deallocate( psi )
462 deallocate( chi )
463 deallocate( temp )
464 deallocate( rh )
465 deallocate( psfc )
466 deallocate( height )
467 deallocate( xlat )
468 deallocate( psi_mean )
469 deallocate( chi_mean )
470 deallocate( temp_mean )
471 deallocate( rh_mean )
472 deallocate( psfc_mean )
473
474 CONTAINS
475
476 !------------------------------------------------------------------------------
477
478 subroutine da_fft_initialize1( dim1, dim2, n1, n2, ifax1, ifax2 )
479
480 implicit none
481
482 integer, intent(in):: dim1, dim2 ! Dimensions.
483 integer, intent(out):: n1, n2 ! Padded dimensions (n=dim-1+pad).
484 integer, intent(out):: ifax1(1:num_fft_factors) ! FFT factors.
485 integer, intent(out):: ifax2(1:num_fft_factors) ! FFT factors.
486
487 integer :: n ! n+1 is the length of the data.
488 integer :: fft_pad1, fft_pad2 ! Range to search for efficient FFT.
489 logical :: found_magic ! True if 2**p 3**p 5**r dimension found..
490
491 integer :: fft_factors(1:num_fft_factors)! FFT factors.
492
493 ! Ensure efficient FFT dimensions by padding if necessary:
494 n1 = dim1 - 1
495 do n = n1, n1 + nrange
496 call da_find_fft_factors( n, found_magic, fft_factors )
497 if ( found_magic .and. mod(n,2) == 0 ) then ! Even magic number found.
498 fft_pad1 = n - n1
499 ifax1 = fft_factors
500 exit
501 end if
502 end do
503 n1 = n1 + fft_pad1
504
505 n2 = dim2 - 1
506 do n = n2, n2 + nrange
507 call da_find_fft_factors( n, found_magic, fft_factors )
508 if ( found_magic .and. mod(n,2) == 0 ) then ! Even magic number found.
509 fft_pad2 = n - n2
510 ifax2 = fft_factors
511 exit
512 end if
513 end do
514 n2 = n2 + fft_pad2
515
516 end subroutine da_fft_initialize1
517
518 !------------------------------------------------------------------------------
519
520 subroutine da_fft_initialize2( n1, n2, ds, trigs1, trigs2, fft_coeffs )
521
522 ! Need to split fft_initialize as array dimensions need to be calculated first.
523
524 implicit none
525
526 integer, intent(in):: n1, n2 ! Padded dimensions (n=dim-1+pad).
527 real, intent(in) :: ds ! Grid resolution.
528 real, intent(out) :: trigs1(1:3*n1) ! FFT trig functions.
529 real, intent(out) :: trigs2(1:3*n2) ! FFT trig functions.
530 real, intent(out) :: fft_coeffs(1:n1+1,1:n2+1) ! FFT coefficients.
531
532 integer :: i, j ! Loop counters.
533 real :: const ! Multiplicative constant.
534 real :: coeff_nx ! Multiplicative constant.
535 real :: coeff_ny ! Multiplicative constant.
536 real :: cos_coeff_nx ! Multiplicative constant.
537 real :: cos_coeff_ny ! Multiplicative constant.
538
539 const = -0.5 * ds * ds
540 coeff_nx = pi / real(n1)
541 coeff_ny = pi / real(n2)
542
543 ! Calculate spectral Del**2 coefficients for C-grid (all pts. except i=j=1):
544 fft_coeffs(1,1) = 0.0 ! Not used?
545 do j = 2, n2+1
546 cos_coeff_ny = cos(coeff_ny * real(j - 1))
547 do i = 1, n1+1
548 cos_coeff_nx = cos(coeff_nx * real(i - 1))
549 fft_coeffs(i,j) = const / ( 2.0 - cos_coeff_nx - cos_coeff_ny)
550 end do
551 end do
552 j = 1
553 cos_coeff_ny = cos(coeff_ny * real(j - 1))
554 do i = 2, n1+1
555 cos_coeff_nx = cos(coeff_nx * real(i - 1))
556 fft_coeffs(i,j) = const / ( 2.0 - cos_coeff_nx - cos_coeff_ny)
557 end do
558
559 call da_find_fft_trig_funcs( n1, trigs1 )
560 call da_find_fft_trig_funcs( n2, trigs2 )
561
562 end subroutine da_fft_initialize2
563
564 !------------------------------------------------------------------------------
565
566 subroutine da_uv_to_div_c( dim1, dim2, ds, &
567 mapfac_m, mapfac_u, mapfac_v, &
568 u, v, div )
569
570 !------------------------------------------------------------------------------
571 ! PURPOSE: Calculate divergence on a co-ordinate surface, given an input
572 ! wind field on an Arakawa C-grid.
573 !
574 ! NOTE: No boundary conditions required on the WRF Arakawa C-grid as
575 ! divergence (mass) points are all within the outer u/v pts.
576 !
577 !
578 ! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker
579 ! d U d V
580 ! Div = m^2 *[---(---) + ---(---) ]
581 ! dx m dy M
582 !------------------------------------------------------------------------------
583
584 implicit none
585
586 integer, intent(in):: dim1, dim2 ! Dimensions.
587 real, intent(in) :: ds ! Resolution.
588 real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts.
589 real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points.
590 real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points.
591 real, intent(in) :: u(1:dim1+1,1:dim2) ! v wind.
592 real, intent(in) :: v(1:dim1,1:dim2+1) ! v wind.
593 real, intent(out) :: div(1:dim1,1:dim2) ! Divergence.
594
595 integer :: i, j ! Loop counters.
596 real :: ds_inv ! 1/ds.
597 real :: coeff(1:dim1,1:dim2) ! Coefficient.
598 real :: um(1:dim1+1,1:dim2) ! u-wind copy.
599 real :: vm(1:dim1,1:dim2+1) ! v-wind copy.
600
601 !------------------------------------------------------------------------------
602 ! [1.0] Initialise:
603 !------------------------------------------------------------------------------
604
605 ds_inv = 1.0 / ds
606 ! do j = 1, dim2
607 ! do i = 1, dim1
608 ! coeff(i,j) = ( mapfac_m(i,j) * mapfac_m(i,j) ) * ds_inv
609 ! end do
610 ! end do
611 coeff(:,:) = ds_inv ! Calculate f.d. Del**2 Chi, rather than divergence.
612
613 !------------------------------------------------------------------------------
614 ! [2] Calculate divergence field:
615 !------------------------------------------------------------------------------
616
617 do j = 1, dim2
618 do i = 1, dim1+1
619 um(i,j) = u(i,j) / mapfac_u(i,j)
620 end do
621 end do
622
623 do j = 1, dim2+1
624 do i = 1, dim1
625 if (mapfac_v(i,j) > 0.000001) then
626 vm(i,j) = v(i,j) / mapfac_v(i,j)
627 else
628 vm(i,j)=0.0
629 end if
630 end do
631 end do
632
633 do j = 1, dim2
634 do i = 1, dim1
635 div(i,j) = coeff(i,j) * ( um(i+1,j) - um(i,j) + vm(i,j+1) - vm(i,j) )
636 end do
637 end do
638
639 end subroutine da_uv_to_div_c
640
641 subroutine da_uv_to_vor_c( dim1, dim2, ds, &
642 mapfac_m, mapfac_u, mapfac_v, &
643 u, v, vor )
644
645 !------------------------------------------------------------------------------
646 ! PURPOSE: Calculate vorticity on a co-ordinate surface, given an input
647 ! wind field on an Arakawa C-grid.
648 !
649 ! NOTE: Zero vorticity boundary conditions.
650 !
651 ! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker
652 ! d V d U
653 ! Vor = m^2 *[---(---) - ---(---) ]
654 ! dx m dy M
655 !------------------------------------------------------------------------------
656
657 implicit none
658
659 integer, intent(in):: dim1, dim2 ! Dimensions.
660 real, intent(in) :: ds ! Resolution.
661 real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts.
662 real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points.
663 real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points.
664 real, intent(in) :: u(1:dim1+1,1:dim2) ! v wind.
665 real, intent(in) :: v(1:dim1,1:dim2+1) ! v wind.
666 real, intent(out) :: vor(1:dim1+1,1:dim2+1) ! Vorticity.
667
668 integer :: i, j ! Loop counters.
669 real :: ds_inv ! 1/ds.
670 ! real :: mapfac_vor ! Map factor (vorticity pts)
671 real :: coeff(1:dim1,1:dim2) ! Coefficient.
672 real :: um(1:dim1+1,1:dim2) ! u-wind copy.
673 real :: vm(1:dim1,1:dim2+1) ! v-wind copy.
674
675 !------------------------------------------------------------------------------
676 ! [1.0] Initialise:
677 !------------------------------------------------------------------------------
678
679 vor(:,:) = 0.0
680
681 ds_inv = 1.0 / ds
682 ! do j = 1, dim2
683 ! do i = 1, dim1
684 ! mapfac_vor = 0.25 * ( mapfac_u(i,j-1) + mapfac_u(i,j) + &
685 ! mapfac_v(i-1,j) + mapfac_v(i,j) ) ! Average.
686 ! coeff(i,j) = ( mapfac_vor * mapfac_vor ) * ds_inv
687 ! end do
688 ! end do
689 coeff(:,:) = ds_inv ! Calculate f.d. Del**2 Chi, rather than divergence.
690
691 !------------------------------------------------------------------------------
692 ! [2] Calculate vorticity field:
693 !------------------------------------------------------------------------------
694
695 do j = 1, dim2
696 do i = 1, dim1+1
697 um(i,j) = u(i,j) / mapfac_u(i,j)
698 end do
699 end do
700
701 do j = 1, dim2+1
702 do i = 1, dim1
703 if (mapfac_v(i,j) > 0.000001) then
704 vm(i,j) = v(i,j) / mapfac_v(i,j)
705 else
706 vm(i,j) = 0.0
707 end if
708 end do
709 end do
710
711 do j = 2, dim2
712 do i = 2, dim1
713 vor(i,j) = coeff(i,j) * ( vm(i,j) - vm(i-1,j) - um(i,j) + um(i,j-1) )
714 end do
715 end do
716
717 ! Boundary values (extrapolation):
718 ! Note - not used in Del**2 calculation if overwritten with bcs there).
719 ! vor(1,1:dim2+1) = 2.0 * vor(2,1:dim2+1) - vor(3,1:dim2+1) ! West
720 ! vor(dim1+1,1:dim2+1) = 2.0 * vor(dim1,1:dim2+1) - vor(dim1-1,1:dim2+1) ! East
721 ! vor(1:dim1+1,1) = 2.0 * vor(1:dim1+1,2) - vor(1:dim1+1,3) ! South
722 ! vor(1:dim1+1,dim2+1) = 2.0 * vor(1:dim1+1,dim2) - vor(1:dim1+1,dim2-1) ! South
723
724 ! Boundary values (zero gradient):
725 ! Note - not used in Del**2 calculation if overwritten with bcs there).
726 vor(1,2:dim2) = vor(2,2:dim2) ! West
727 vor(dim1+1,2:dim2) = vor(dim1,2:dim2) ! East
728 vor(1:dim1+1,1) = vor(1:dim1+1,2) ! South
729 vor(1:dim1+1,dim2+1) = vor(1:dim1+1,dim2) ! South
730
731 end subroutine da_uv_to_vor_c
732
733 subroutine da_psichi_to_uv_c( dim1, dim2, ds, &
734 mapfac_u, mapfac_v, &
735 psi, chi, u, v )
736
737 !------------------------------------------------------------------------------
738 ! PURPOSE: Calculate u and v wind components on an Arakawa C-grid.
739 !
740 ! NOTE:
741 !
742 ! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker
743 !------------------------------------------------------------------------------
744
745 implicit none
746
747 integer, intent(in):: dim1, dim2 ! Dimensions.
748 real, intent(in) :: ds ! Resolution.
749 real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points.
750 real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points.
751 real, intent(in) :: psi(1:dim1+1,1:dim2+1) ! Streamfunction.
752 real, intent(in) :: chi(1:dim1,1:dim2) ! Velcoity potential.
753 real, intent(out) :: u(1:dim1+1,1:dim2) ! v wind.
754 real, intent(out) :: v(1:dim1,1:dim2+1) ! v wind.
755
756 integer :: i, j ! Loop counters.
757 integer :: its, ite, jts, jte ! WRF dims (dummies for now).
758 real :: ds_inv ! 1/ds.
759 real :: one_third ! 1/3.
760
761 ds_inv = 1.0 / ds
762 one_third = 1.0 / 3.0
763
764 ! u-wind component:
765 its = 2
766 ite = dim1
767 jts = 1
768 jte = dim2
769
770 do j = jts, jte
771 do i = its, ite
772 u(i,j) = mapfac_u(i,j) * ds_inv * &
773 ( -psi(i,j+1) + psi(i,j) + chi(i,j) - chi(i-1,j) )
774 end do
775 end do
776
777 ! Remaining points on E/W boundaries (extrapolation):
778 u(1,jts:jte) = 2.0 * u(2,jts:jte) - u(3,jts:jte)
779 u(dim1+1,jts:jte) = 2.0 * u(dim1,jts:jte) - u(dim1-1,jts:jte)
780
781 ! v-wind component:
782 its = 1
783 ite = dim1
784 jts = 2
785 jte = dim2
786
787 do j = jts, jte
788 do i = its, ite
789 v(i,j) = mapfac_v(i,j) * ds_inv * &
790 ( psi(i+1,j) - psi(i,j) + chi(i,j) - chi(i,j-1) )
791 end do
792 end do
793
794 ! Remaining points on S/N boundaries (extrapolation):
795 v(its:ite,1) = 2.0 * v(its:ite,2) - v(its:ite,3)
796 v(its:ite,dim2+1) = 2.0 * v(its:ite,dim2) - v(its:ite,dim2-1)
797
798 end subroutine da_psichi_to_uv_c
799
800 !------------------------------------------------------------------------------
801
802 subroutine da_del2a_to_a( dim1, dim2, n1, n2, fft_method, ifax1, ifax2, trigs1, trigs2, &
803 fft_coeffs, del2a, a )
804
805 implicit none
806
807 integer, intent(in):: dim1, dim2 ! Dimensions.
808 integer, intent(in):: n1, n2 ! Padded dimensions (n=dim-1+pad).
809 integer, intent(in):: fft_method ! 1=Cosine, 2=Sine transform.
810 integer, intent(in):: ifax1(1:num_fft_factors) ! FFT factors.
811 integer, intent(in):: ifax2(1:num_fft_factors) ! FFT factors.
812 real, intent(in) :: trigs1(1:3*n1) ! FFT trig functions.
813 real, intent(in) :: trigs2(1:3*n2) ! FFT trig functions.
814 real, intent(in) :: fft_coeffs(1:n1+1,1:n2+1) ! FFT coefficients.
815 real, intent(in) :: del2a(1:dim1,1:dim2) ! Del**2 a.
816 real, intent(out) :: a(1:dim1,1:dim2) ! Field a.
817
818 integer :: i, j ! Loop counters.
819 integer :: ij ! 1D array counter.
820 integer :: isign ! -1=Grid>spec, 1=Spec>Grid.
821 integer :: inc ! Stride between data points.
822 integer :: jump ! Increment between start of data vectors.
823 integer :: lot ! Number of data vectors.
824 integer :: n ! n+1 is the length of the data.
825 integer :: work_area ! Dimension of workspace.
826 real :: a2d(1:n1+1,1:n2+1) ! 2D data array.
827 real :: a1d(1:(n1+1)*(n2+1)) ! 1D data array.
828
829 work_area = ( n1 + 1 ) * ( n2 + 1 )
830
831 ! Fill 2D array structure
832 do j = 1, dim2
833 do i = 1, dim1
834 a2d(i,j) = del2a(i,j)
835 end do
836
837 ! Fill pad zone (and force b.c.s to satisfy solution type):
838 if ( fft_method == 1 ) then ! Cosine transform.
839 a2d(1,j) = a2d(2,j)
840 do i = dim1, n1+1
841 a2d(i,j) = a2d(dim1-1,j)
842 end do
843 else if ( fft_method == 2 ) then ! Sine transform:
844 a2d(1,j) = 0.0
845 do i = dim1, n1+1
846 a2d(i,j) = 0.0
847 end do
848 end if
849 end do
850
851 if ( fft_method == 1 ) then ! Cosine transform.
852 do i = 1, n1+1
853 a2d(i,1) = a2d(i,2)
854 do j = dim2, n2+1
855 a2d(i,j) = a2d(i,dim2-1)
856 end do
857 end do
858 else if ( fft_method == 2 ) then ! Sine transform:
859 do i = 1, n1+1
860 a2d(i,1) = 0.0
861 do j = dim2, n2+1
862 a2d(i,j) = 0.0
863 end do
864 end do
865 end if
866
867 ! Transfer to data array:
868 do j = 1, n2+1
869 do i = 1, n1+1
870 ij = (j-1) * (n1+1) + i
871 a1d(ij) = a2d(i,j)
872 end do
873 end do
874
875 !------------------------------------------------------------------------------
876 ! Perform double fast sine/cosine transform to get spectral del2a:
877 !------------------------------------------------------------------------------
878
879 isign = -1 ! Grid to spectral
880
881 ! 1st dimension:
882 inc = 1 ! Stride between data points.
883 jump = n1+1! Increment between start of data vectors.
884 lot = n2+1 ! Number of data vectors.
885 n = n1 ! n+1 is the length of the data.
886 if ( fft_method == 1 ) then
887 call fft551( isign, inc, jump, lot, n, &
888 ifax1, trigs1, a1d, work_area )
889 else if ( fft_method == 2 ) then
890 call fft661( isign, inc, jump, lot, n, &
891 ifax1, trigs1, a1d, work_area )
892 end if
893
894 ! 2nd dimension:
895 inc = n1+1 ! Stride between data points.
896 jump = 1 ! Increment between start of data vectors.
897 lot = n1+1 ! Number of data vectors.
898 n = n2 ! n+1 is the length of the data.
899
900 if ( fft_method == 1 ) then
901 call fft551( isign, inc, jump, lot, n, &
902 ifax2, trigs2, a1d, work_area )
903 else if ( fft_method == 2 ) then
904 call fft661( isign, inc, jump, lot, n, &
905 ifax2, trigs2, a1d, work_area )
906 end if
907
908 !------------------------------------------------------------------------------
909 ! Perform conversion from del2a to a in spectral space:
910 !------------------------------------------------------------------------------
911
912 ! Note fft_coeffs(1,1)=0 so a(k=0,l=0) is also 0.
913 do j = 1, n2+1
914 do i = 1, n1+1
915 ij = (j-1) * (n1+1) + i
916 a1d(ij) = fft_coeffs(i,j) * a1d(ij)
917 end do
918 end do
919
920 !------------------------------------------------------------------------------
921 ! Perform double fast sine/cosine transform to get gridpoint a:
922 !------------------------------------------------------------------------------
923
924 isign = 1 ! Spectral to grid.
925
926 ! 1st dimension:
927 inc = 1 ! Stride between data points.
928 jump = n1+1! Increment between start of data vectors.
929 lot = n2+1 ! Number of data vectors.
930 n = n1 ! n+1 is the length of the data.
931
932 if ( fft_method == 1 ) then
933 call fft551( isign, inc, jump, lot, n, &
934 ifax1, trigs1, a1d, work_area )
935 else if ( fft_method == 2 ) then
936 call fft661( isign, inc, jump, lot, n, &
937 ifax1, trigs1, a1d, work_area )
938 end if
939
940 ! 2nd dimension:
941 inc = n1+1 ! Stride between data points.
942 jump = 1 ! Increment between start of data vectors.
943 lot = n1+1 ! Number of data vectors.
944 n = n2 ! n+1 is the length of the data.
945
946 if ( fft_method == 1 ) then
947 call fft551( isign, inc, jump, lot, n, &
948 ifax2, trigs2, a1d, work_area )
949 else if ( fft_method == 2 ) then
950 call fft661( isign, inc, jump, lot, n, &
951 ifax2, trigs2, a1d, work_area )
952 end if
953
954 ! Transfer grid-point chi to 2D-array (throwing away pad):
955 do j = 1, dim2
956 do i = 1, dim1
957 ij = (j-1) * (n1+1) + i
958 a(i,j) = a1d(ij)
959 end do
960 end do
961
962 end subroutine da_del2a_to_a
963
964 !---------------------------------------------------------------------------------------------
965 subroutine da_test_inverse( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
966 n1, n2, fft_method, ifax1, ifax2, trigs1, trigs2, fft_coeffs, &
967 n1s, n2s, ifax1s, ifax2s, trigs1s, trigs2s, fft_coeffss, &
968 u1, v1, psi, chi )
969
970 !------------------------------------------------------------------------------
971 ! PURPOSE: Test u, v -> psi, chi calculation by performing inverse test.
972 !
973 ! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker
974 !------------------------------------------------------------------------------
975
976 implicit none
977
978 integer, intent(in):: dim1, dim2 ! Dimensions.
979 real, intent(in) :: ds ! Resolution.
980 real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts.
981 real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points.
982 real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points.
983 integer, intent(in):: n1, n2 ! Padded dimensions (n=dim-1+pad).
984 integer, intent(in):: fft_method ! 1=Cosine, 2=Sine transform.
985 integer, intent(in):: ifax1(1:num_fft_factors) ! FFT factors.
986 integer, intent(in):: ifax2(1:num_fft_factors) ! FFT factors.
987 real, intent(in) :: trigs1(1:3*n1) ! FFT trig functions.
988 real, intent(in) :: trigs2(1:3*n2) ! FFT trig functions.
989 real, intent(in) :: fft_coeffs(1:n1+1,1:n2+1) ! FFT coefficients.
990 integer, intent(in):: n1s, n2s ! Padded dimensions (n=dim-1+pad).
991 integer, intent(in):: ifax1s(1:num_fft_factors) ! FFT factors.
992 integer, intent(in):: ifax2s(1:num_fft_factors) ! FFT factors.
993 real, intent(in) :: trigs1s(1:3*n1) ! FFT trig functions.
994 real, intent(in) :: trigs2s(1:3*n2) ! FFT trig functions.
995 real, intent(in) :: fft_coeffss(1:n1+1,1:n2+1) ! FFT coefficients.
996
997 real, intent(in) :: u1(1:dim1+1,1:dim2) ! u
998 real, intent(in) :: v1(1:dim1,1:dim2+1) ! v
999 real, intent(in) :: psi(1:dim1+1,1:dim2+1) ! Streamfunction.
1000 real, intent(in) :: chi(1:dim1,1:dim2) ! Velocity potential.
1001
1002 real :: div(1:dim1,1:dim2) ! Divergence.
1003 real :: vor(1:dim1+1,1:dim2+1) ! Vorticity.
1004
1005 real :: u2(1:dim1+1,1:dim2) ! u
1006 real :: v2(1:dim1,1:dim2+1) ! v
1007 real :: u3(1:dim1+1,1:dim2) ! u
1008 real :: v3(1:dim1,1:dim2+1) ! v
1009 real :: psi1(1:dim1+1,1:dim2+1) ! streamfunction
1010 real :: chi1(1:dim1,1:dim2) ! divergence
1011
1012 write(6,'(a,i4)')' Using FFT method (1=Cosine, 2=Sine): ', fft_method
1013
1014 call da_psichi_to_uv_c( dim1, dim2, ds, &
1015 mapfac_u, mapfac_v, psi, chi, u2, v2 )
1016
1017 write(6,'(a,1pe12.4)')' Inverse test 1: Ratio error/field u = ', &
1018 sqrt(sum( ( u1(:,:) - u2(:,:) )**2 ) / sum( u1(:,:)**2 ))
1019 write(6,'(a,1pe12.4)')' Inverse test 1: Ratio error/field v = ', &
1020 sqrt(sum( ( v1(:,:) - v2(:,:) )**2 ) / sum( v1(:,:)**2 ))
1021
1022 call da_uv_to_div_c( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
1023 u2, v2, div )
1024 call da_uv_to_vor_c( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
1025 u2, v2, vor )
1026
1027 call da_del2a_to_a( dim1, dim2, n1, n2, fft_method, ifax1, ifax2, trigs1, trigs2, &
1028 fft_coeffs, div, chi1 )
1029 call da_del2a_to_a( dim1s, dim2s, n1s, n2s, fft_method, ifax1s, ifax2s, trigs1s, trigs2s, &
1030 fft_coeffss, vor, psi1 )
1031
1032 write(6,'(a,1pe12.4)')' Inverse test 2: Ratio error/field psi = ', &
1033 sqrt(sum( ( psi(:,:) - psi1(:,:) )**2 ) / sum( psi(:,:)**2 ))
1034 write(6,'(a,1pe12.4)')' Inverse test 2: Ratio error/field chi = ', &
1035 sqrt(sum( ( chi(:,:) - chi1(:,:) )**2 ) / sum( chi(:,:)**2 ))
1036
1037 call da_psichi_to_uv_c( dim1, dim2, ds, &
1038 mapfac_u, mapfac_v, psi1, chi1, u3, v3 )
1039
1040 write(6,'(a,1pe12.4)')' Inverse test 3: Ratio error/field u = ', &
1041 sqrt(sum( ( u3(:,:) - u2(:,:) )**2 ) / sum( u2(:,:)**2 ))
1042 write(6,'(a,1pe12.4)')' Inverse test 3: Ratio error/field v = ', &
1043 sqrt(sum( ( v3(:,:) - v2(:,:) )**2 ) / sum( v2(:,:)**2 ))
1044
1045 end subroutine da_test_inverse
1046
1047 !---------------------------------------------------------------------------------------------
1048
1049 #if 0
1050 ! routine never used
1051 subroutine da_test_inverse2( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
1052 u1, v1, psi, chi )
1053
1054 !------------------------------------------------------------------------------
1055 ! PURPOSE: Test u, v -> psi, chi calculation by performing inverse test.
1056 !
1057 ! HISTORY: 08/09/2005 - Creation of F90 version. Dale Barker
1058 !------------------------------------------------------------------------------
1059
1060 implicit none
1061
1062 integer, intent(in):: dim1, dim2 ! Dimensions.
1063 real, intent(in) :: ds ! Resolution.
1064 real, intent(in) :: mapfac_m(1:dim1,1:dim2) ! Map factor - mass pts.
1065 real, intent(in) :: mapfac_u(1:dim1+1,1:dim2) ! Map factor - u points.
1066 real, intent(in) :: mapfac_v(1:dim1,1:dim2+1) ! Map factor - u points.
1067
1068 real, intent(in) :: u1(1:dim1+1,1:dim2) ! u
1069 real, intent(in) :: v1(1:dim1,1:dim2+1) ! v
1070 real, intent(in) :: psi(1:dim1+1,1:dim2+1) ! Streamfunction.
1071 real, intent(in) :: chi(1:dim1,1:dim2) ! Velocity potential.
1072
1073 real :: div(1:dim1,1:dim2) ! Divergence.
1074 real :: vor(1:dim1+1,1:dim2+1) ! Vorticity.
1075
1076 real :: u2(1:dim1+1,1:dim2) ! u
1077 real :: v2(1:dim1,1:dim2+1) ! v
1078 real :: u3(1:dim1+1,1:dim2) ! u
1079 real :: v3(1:dim1,1:dim2+1) ! v
1080 real :: psi1(1:dim1+1,1:dim2+1) ! streamfunction
1081 real :: chi1(1:dim1,1:dim2) ! divergence
1082
1083 write(6,'(a,i4)')' Using SOR method'
1084
1085 call da_psichi_to_uv_c( dim1, dim2, ds, &
1086 mapfac_u, mapfac_v, psi, chi, u2, v2 )
1087
1088 write(6,'(a,1pe12.4)')' Inverse test 1: Ratio error/field u = ', &
1089 sqrt(sum( ( u1(:,:) - u2(:,:) )**2 ) / sum( u1(:,:)**2 ))
1090 write(6,'(a,1pe12.4)')' Inverse test 1: Ratio error/field v = ', &
1091 sqrt(sum( ( v1(:,:) - v2(:,:) )**2 ) / sum( v1(:,:)**2 ))
1092
1093 call da_uv_to_div_c( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
1094 u2, v2, div )
1095 call da_uv_to_vor_c( dim1, dim2, ds, mapfac_m, mapfac_u, mapfac_v, &
1096 u2, v2, vor )
1097
1098 call da_sor( dim1, dim2, ds, div, chi1 )
1099 call da_sor( dim1s, dim2s, ds, vor, psi1 )
1100
1101 write(6,'(a,1pe12.4)')' Inverse test 2: Ratio error/field psi = ', &
1102 sqrt(sum( ( psi(:,:) - psi1(:,:) )**2 ) / sum( psi(:,:)**2 ))
1103 write(6,'(a,1pe12.4)')' Inverse test 2: Ratio error/field chi = ', &
1104 sqrt(sum( ( chi(:,:) - chi1(:,:) )**2 ) / sum( chi(:,:)**2 ))
1105
1106 call da_psichi_to_uv_c( dim1, dim2, ds, &
1107 mapfac_u, mapfac_v, psi1, chi1, u3, v3 )
1108
1109 write(6,'(a,1pe12.4)')' Inverse test 3: Ratio error/field u = ', &
1110 sqrt(sum( ( u3(:,:) - u2(:,:) )**2 ) / sum( u2(:,:)**2 ))
1111 write(6,'(a,1pe12.4)')' Inverse test 3: Ratio error/field v = ', &
1112 sqrt(sum( ( v3(:,:) - v2(:,:) )**2 ) / sum( v2(:,:)**2 ))
1113
1114 end subroutine da_test_inverse2
1115 #endif
1116
1117 !----------------------------------------------------------------------------------------
1118
1119 !DALE: Use fortran 90 free-form, and lower case throughout?
1120 !DALE: Provide short description of code here, including assumptions, ref., etc.
1121
1122 SUBROUTINE da_sor ( dim1, dim2, ds, del2a, a )
1123
1124 IMPLICIT NONE
1125
1126 REAL, PARAMETER :: SMALLRES = 1.0E-4 ! DALE: Provide comments here.
1127
1128 INTEGER, PARAMETER :: MM = 20000 ! DALE: And here, etc.
1129 REAL, PARAMETER :: ALPHA = 1.8
1130 REAL, PARAMETER :: ALPHAOV4 = alpha / 4.0
1131
1132 integer, intent(in):: dim1, dim2 ! Dimensions.
1133 real, intent(in) :: ds ! Grid resolution.
1134 real, intent(in) :: del2a(1:dim1,1:dim2) ! Del**2 a.
1135 real, intent(out) :: a(1:dim1,1:dim2) ! Field a.
1136
1137 integer :: i, j ! Loop counters.
1138 real :: rd(1:dim1,1:dim2) ! Solution residual.
1139
1140 integer :: ids, ide, jds, jde
1141 INTEGER :: ITER
1142 INTEGER :: MI
1143 REAL :: CHI ( 1:dim1,1:dim2 )
1144 REAL :: CHIMX ( 1:dim1 )
1145 REAL :: EPX
1146 REAL :: FAC
1147 REAL :: FF ( 1:dim1 , 1:dim2 )
1148 REAL :: RDMAX ( 1:dim1 )
1149
1150 LOGICAL :: converged = .FALSE.
1151
1152 ids = 1
1153 ide = dim1
1154 jds = 1
1155 jde = dim2
1156
1157 rd = 0.0
1158 chi = 0.0
1159
1160 fac = 2.0 * ds * ds
1161 DO i = ids, ide
1162 DO j = jds, jde
1163 ff(i,j) = del2a(i,j) * fac
1164 END DO
1165 END DO
1166
1167 iter_loop : DO iter = 1, mm
1168 mi = iter
1169 chimx = 0.0
1170
1171 ! FIX? Crayx1 compiler does not like these, so comment out
1172 ! for the present
1173 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1174 !!csd$ parallel do private(i,j)
1175 DO i = ids+1, ide-1
1176 DO j = jds+1, jde-1
1177 chimx(i) = MAX(ABS(chi(i,j)),chimx(i))
1178 END DO
1179 END DO
1180 !!csd$ end parallel do
1181
1182 epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha
1183
1184 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1185 !!csd$ parallel do private(i,j)
1186 DO j = jds+1, jde-1
1187 DO i = ids+1, ide-1
1188 rd(i,j) = chi(i+1,j+1) + chi(i-1,j+1) + chi(i+1,j-1) + chi(i-1,j-1) - 4.0 * chi(i,j) - ff(i,j)
1189 chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1190 END DO
1191 END DO
1192 !!csd$ end parallel do
1193
1194 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1195 !!csd$ parallel do private(i,j)
1196 DO j = jde-1, jds+1, -1
1197 DO i = ide-1, ids+1, -1
1198 rd(i,j) = chi(i+1,j+1) + chi(i-1,j+1) + chi(i+1,j-1) + chi(i-1,j-1) - 4.0 * chi(i,j) - ff(i,j)
1199 chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1200 END DO
1201 END DO
1202 !!csd$ end parallel do
1203
1204 rdmax = 0.0
1205
1206 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1207 !!csd$ parallel do private(i,j)
1208 DO i = ids+1, ide-1
1209 DO j = jds+1, jde-1
1210 rdmax(i) = MAX(ABS(rd(i,j)),rdmax(i))
1211 END DO
1212 END DO
1213 !!csd$ end parallel do
1214
1215 IF (MAXVAL(rdmax) .lt. epx) THEN
1216 converged = .TRUE.
1217 EXIT iter_loop
1218 ELSE
1219 ! print '("iter No.",i5," Max. Residule=",e12.5)', mi, maxval(rdmax)
1220 END IF
1221
1222 END DO iter_loop
1223
1224 IF (converged ) THEN
1225 PRINT '(A,i3,A,I5,A)','k=',k,' Relaxation converged in ',mi,' iterations.'
1226 ELSE
1227 PRINT '(A,i3,A,I5,A)','k=',k,' Relaxation did not converge in',mm,' iterations.'
1228 ! STOP 'no_converge'
1229 END IF
1230
1231 a(:,:) = chi(:,:)
1232
1233 END SUBROUTINE da_sor
1234
1235 #if 0
1236 ! routine never used
1237 SUBROUTINE da_relax (psi, vor, residual, ids, ide, jds, jde, ds)
1238
1239 IMPLICIT NONE
1240
1241 INTEGER, intent(in) :: IDS, IDE
1242 INTEGER, intent(in) :: JDS, JDE
1243 REAL, intent(out) :: PSI ( ids:ide , jds:jde )
1244 REAL, intent(in) :: VOR ( ids:ide , jds:jde )
1245 REAL, intent(out) :: residual ( ids:ide , jds:jde )
1246 REAL, intent(in) :: DS
1247
1248 REAL, PARAMETER :: SMALLRES = 5.0E-7
1249 ! REAL, PARAMETER :: SMALLRES = 1.0E-6
1250 ! REAL, PARAMETER :: SMALLRES = 1.0E-4 ! DALE: Provide comments here.
1251
1252 INTEGER, PARAMETER :: MM = 20000 ! DALE: And here, etc.
1253 ! REAL, PARAMETER :: ALPHA = 1.8
1254 REAL, PARAMETER :: ALPHA = 1.2
1255 REAL, PARAMETER :: ALPHAOV4 = alpha / 4.0
1256
1257 INTEGER :: I
1258 INTEGER :: ITER
1259 INTEGER :: J
1260 INTEGER :: MI
1261
1262 REAL (kind=8) :: CHI ( ids:ide , jds:jde )
1263 REAL (kind=8) :: CHIMX ( ids:ide )
1264 REAL (kind=8) :: EPX
1265 REAL (kind=8) :: FAC
1266 REAL (kind=8) :: FF ( ids:ide , jds:jde )
1267 REAL (kind=8) :: RD ( ids:ide , jds:jde )
1268 REAL (kind=8) :: RDMAX ( ids:ide )
1269
1270 LOGICAL :: converged = .FALSE.
1271
1272 rd = 0.0
1273 chi = 0.0
1274
1275 fac = 2.0 * ds * ds
1276 DO i = ids, ide
1277 DO j = jds, jde
1278 ff(i,j) = vor(i,j) * fac
1279 END DO
1280 END DO
1281
1282 iter_loop : DO iter = 1, mm
1283 mi = iter
1284 chimx = 0.0
1285
1286
1287 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1288 !!csd$ parallel do private(i,j)
1289 DO i = ids+1, ide-1
1290 DO j = jds+1, jde-1
1291 chimx(i) = MAX(ABS(chi(i,j)),chimx(i))
1292 END DO
1293 END DO
1294 !!csd$ end parallel do
1295
1296 epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha
1297
1298 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1299 !!csd$ parallel do private(i,j)
1300 DO j = jds+1, jde-1
1301 DO i = ids+1, ide-1
1302 rd(i,j) = chi(i+1,j+1) + chi(i-1,j+1) + chi(i+1,j-1) + chi(i-1,j-1) - 4.0 * chi(i,j) - ff(i,j)
1303 chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1304 END DO
1305 END DO
1306 !!csd$ end parallel do
1307
1308 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1309 !!csd$ parallel do private(i,j)
1310 DO j = jde-1, jds+1, -1
1311 DO i = ide-1, ids+1, -1
1312 rd(i,j) = chi(i+1,j+1) + chi(i-1,j+1) + chi(i+1,j-1) + chi(i-1,j-1) - 4.0 * chi(i,j) - ff(i,j)
1313 chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1314 END DO
1315 END DO
1316 !!csd$ end parallel do
1317
1318 rdmax = 0.0
1319
1320 !!$OMP PARALLEL DO DEFAULT ( SHARED ) PRIVATE ( i , j )
1321 !!csd$ parallel do private(i,j)
1322 DO i = ids+1, ide-1
1323 DO j = jds+1, jde-1
1324 rdmax(i) = MAX(ABS(rd(i,j)),rdmax(i))
1325 END DO
1326 END DO
1327 !!csd$ end parallel do
1328
1329 IF (MAXVAL(rdmax) .lt. epx) THEN
1330 converged = .TRUE.
1331 EXIT iter_loop
1332 ELSE
1333 ! print '("iter No.",i5," Max. Residule=",e12.5)', mi, maxval(rdmax)
1334 END IF
1335
1336 END DO iter_loop
1337
1338 IF (converged ) THEN
1339 PRINT '(A,i3,A,I5,A)','k=',k,' Relaxation converged in ',mi,' iterations.'
1340 ELSE
1341 PRINT '(A,i3,A,I5,A)','k=',k,' Relaxation did not converge in',mm,' iterations.'
1342 ! STOP 'no_converge'
1343 END IF
1344
1345 psi(:,:) = chi(:,:)
1346 residual(:,:) = rd(:,:)
1347
1348 END SUBROUTINE da_relax
1349 #endif
1350
1351 #ifdef crayx1
1352
1353 subroutine getarg(i, harg)
1354 implicit none
1355 character(len=*) :: harg
1356 integer :: ierr, ilen, i
1357
1358 call pxfgetarg(i, harg, ilen, ierr)
1359 return
1360 end subroutine getarg
1361 #endif
1362
1363 end program gen_be_stage0_wrf
1364